annotate deseq-hts_1.0/mex/interval_overlap.cpp @ 2:8ac6832cb42c draft

Deleted selected files
author vipints
date Wed, 27 Jun 2012 15:37:03 -0400
parents 94a108763d9e
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
1 /*
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
2 * This program is free software; you can redistribute it and/or modify
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
3 * it under the terms of the GNU General Public License as published by
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
4 * the Free Software Foundation; either version 3 of the License, or
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
5 * (at your option) any later version.
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
6 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
7 * Written (W) 2010-2011 Jonas Behr
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
8 * Copyright (C) 2010-2011 Max Planck Society
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
9 */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
10
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
11
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
12 #include <stdio.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
13 #include <stdarg.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
14 #include <errno.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
15 #include <sys/types.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
16 #include <sys/stat.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
17 #include <fcntl.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
18 #include <ctype.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
19 #include <sys/stat.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
20 #include <stdlib.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
21 #include <unistd.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
22 #include <string.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
23 #include <sys/mman.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
24 #include <time.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
25 #include <math.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
26 #include <limits>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
27 #include <mex.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
28 #include <assert.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
29 #include <vector>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
30 using std::vector;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
31 #include <algorithm>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
32 using std::sort;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
33 using std::min;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
34 using std::max;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
35
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
36 typedef struct {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
37 int start;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
38 int stop;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
39 int idx;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
40 int set_id;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
41 } interval_t;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
42
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
43 bool compare (interval_t i, interval_t j)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
44 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
45 return (i.start<j.start);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
46 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
47
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
48 bool overlaps(interval_t a, interval_t b)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
49 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
50 int v = min(a.stop,b.stop) - max(a.start,b.start) + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
51 return (v >= 1);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
52 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
53 bool leftOf(interval_t a, interval_t b)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
54 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
55 return (a.stop < b.start);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
56 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
57
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
58 void scan(interval_t f, vector<interval_t>* Wf, interval_t g, vector<interval_t>* Wg, vector<int>* overlap)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
59 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
60 vector<interval_t>::iterator i;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
61 i=Wg->begin();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
62 while (i<Wg->end())
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
63 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
64 interval_t g2 = *i;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
65 if (leftOf(g2,f))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
66 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
67 Wg->erase(i);// inefficient if Wg is large
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
68 // this moves all elements, therefore i is not incremented
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
69 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
70 else if (overlaps(g2,f))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
71 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
72 if (g2.set_id==1)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
73 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
74 //printf("overlap: [%i | %i, %i] [%i | %i, %i]\n", g2.idx, g2.start, g2.stop, f.idx, f.start, f.stop);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
75 overlap->push_back(g2.idx);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
76 overlap->push_back(f.idx);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
77 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
78 else if (f.set_id==1)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
79 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
80 //printf("overlap: [%i | %i, %i] [%i | %i, %i]\n", f.idx, f.start, f.stop, g2.idx, g2.start, g2.stop);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
81 overlap->push_back(f.idx);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
82 overlap->push_back(g2.idx);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
83 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
84 i++;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
85 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
86 else
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
87 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
88 printf("never happens??\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
89 i++;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
90 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
91 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
92 if (!leftOf(f, g))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
93 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
94 Wf->push_back(f);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
95 //printf("push: [%i, %i] size:%i\n", f.start, f.stop, Wf->size());
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
96 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
97 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
98
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
99 /*
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
100 * prhs[0] first interval set starts
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
101 * prhs[1] first interval set stops
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
102 * prhs[2] second interval set starts
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
103 * prhs[3] second interval set stops
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
104 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
105 * return:
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
106 * plhs[0] one based index in first interval set overlapping with a interval in the second set
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
107 * plhs[1] corresponding index in the second set
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
108 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
109 */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
110 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
111 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
112 if (nrhs!=4)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
113 mexErrMsgTxt("Expected 4 arguments: starts1, stops1, starts2, stops2 \n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
114 if (nlhs!=2)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
115 mexErrMsgTxt("Expected 2 output arguments \n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
116
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
117 int num_intervals1 = mxGetNumberOfElements(prhs[0]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
118 assert(num_intervals1 == mxGetNumberOfElements(prhs[1]));
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
119 int num_intervals2 = mxGetNumberOfElements(prhs[2]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
120 assert(num_intervals2 == mxGetNumberOfElements(prhs[3]));
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
121
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
122 //printf("num_intervals1: %i\n", num_intervals1);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
123 //printf("num_intervals2: %i\n", num_intervals2);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
124
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
125 double* starts1 = mxGetPr(prhs[0]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
126 double* stops1 = mxGetPr(prhs[1]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
127 double* starts2 = mxGetPr(prhs[2]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
128 double* stops2 = mxGetPr(prhs[3]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
129
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
130 vector<interval_t> intervals1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
131 for (int i=0; i<num_intervals1; i++)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
132 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
133 interval_t interval;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
134 interval.start = starts1[i];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
135 interval.stop = stops1[i];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
136 interval.set_id = 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
137 interval.idx = i+1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
138 intervals1.push_back(interval);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
139 //printf("int1: [%i, %i] \n",intervals1[i].start, intervals1[i].stop);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
140 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
141 interval_t i;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
142 i.start = std::numeric_limits<int>::max();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
143 i.stop = std::numeric_limits<int>::max();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
144 i.set_id = std::numeric_limits<int>::max();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
145 i.idx = std::numeric_limits<int>::max();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
146 intervals1.push_back(i);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
147
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
148 //printf("num_intervals1: %i\n", intervals1.size());
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
149 vector<interval_t> intervals2;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
150 for (int i=0; i<num_intervals2; i++)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
151 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
152 interval_t interval;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
153 interval.start = starts2[i];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
154 interval.stop = stops2[i];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
155 interval.set_id = 2;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
156 interval.idx = i+1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
157 intervals2.push_back(interval);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
158 //printf("int2: [%i, %i] \n",intervals2[i].start, intervals2[i].stop);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
159 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
160 intervals2.push_back(i);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
161 //printf("num_intervals2: %i\n", intervals2.size());
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
162
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
163 sort(intervals1.begin(), intervals1.end(), compare);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
164 sort(intervals2.begin(), intervals2.end(), compare);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
165
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
166
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
167 vector<int> overlap;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
168 vector<interval_t> Wx;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
169 vector<interval_t> Wy;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
170 vector<interval_t>::iterator x = intervals1.begin();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
171 vector<interval_t>::iterator y = intervals2.begin();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
172 while (x<intervals1.end() && y<intervals2.end())
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
173 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
174 //vector<interval_t>::iterator x;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
175 //vector<interval_t>::iterator y;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
176 //if (it1>intervals1.end())
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
177 // x = inf_interval();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
178 //else
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
179 // x = it1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
180 //if (it2>intervals2.end())
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
181 // y = inf_interval();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
182 //else
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
183 // y=it2;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
184
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
185 if (x->start <= y->start)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
186 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
187 scan(*x, &Wx, *y, &Wy, &overlap);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
188 x++;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
189 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
190 else
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
191 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
192 if (y<=intervals2.end())
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
193 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
194 scan(*y, &Wy, *x, &Wx, &overlap);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
195 y++;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
196 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
197 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
198 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
199
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
200 plhs[0] = mxCreateDoubleMatrix(1, overlap.size()/2, mxREAL);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
201 double* idx1 = mxGetPr(plhs[0]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
202
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
203 plhs[1] = mxCreateDoubleMatrix(1, overlap.size()/2, mxREAL);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
204 double* idx2 = mxGetPr(plhs[1]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
205
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
206 for (int i=0; i<overlap.size(); i+=2)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
207 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
208 //printf("ov: %i [%i, %i] \n", i, overlap[i], overlap[i+1]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
209 idx1[i/2] = (double) overlap[i];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
210 idx2[i/2] = (double) overlap[i+1];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
211 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
212 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
213
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
214
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
215
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
216
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
217