1 #include "Util.h"
2
3 extern short *pair_table;
4 //extern void *space(unsigned int size);
5 /*@exits@*/
6 //extern void nrerror(const char message[]);
7
8 extern "C" {
9 #include "utils.h"
10 }
11
Cout(std::string s)12 void Cout(std::string s){std::cout<<s;}
13
Str(double x)14 std::string Str(double x)
15 {
16 std::stringstream ss;
17 ss << x;
18 return ss.str();
19 // return Str((float)x);
20 /*
21 std::stringstream ss;
22 char cBuf[1000];
23 sprintf(cBuf,"%12.3f",x);
24 ss << cBuf;
25
26 return ss.str();*/
27 }
28
29
Str(float x)30 std::string Str(float x)
31 {
32 std::stringstream ss;
33 ss << x;
34 return ss.str();
35 }
36
Str(int x)37 std::string Str(int x)
38 {
39 std::stringstream ss;
40 ss << x;
41 return ss.str();
42 }
43
Str(unsigned int x)44 std::string Str(unsigned int x)
45 {
46 std::stringstream ss;
47 ss << x;
48 return ss.str();
49 }
50
PrintPairTable()51 std::string PrintPairTable(){
52 std::string s=std::string();
53 for(int i=0;i<=pair_table[0];i++){
54 s+=Str(pair_table[i])+" ";
55 }
56 s+="\n";
57 return s;
58 }
59
60
PrintBasePair(std::pair<int,int> bp)61 std::string PrintBasePair(std::pair<int,int> bp){
62 return "("+Str(bp.first)+","+Str(bp.second)+")";
63 }
64
65
PrintBasePairList(std::vector<std::pair<int,int>> list)66 std::string PrintBasePairList(std::vector<std::pair<int,int> > list){
67 std::string sText="";
68 for(size_t i=0;i<list.size();i++) sText+=PrintBasePair(list[i]);
69 return sText;
70 }
71
72
73 std::string
PairTableToStructure(std::vector<int> pTbl)74 PairTableToStructure(std::vector<int> pTbl)
75 {
76 std::string structure(pTbl.size(), '.');
77
78 int i=0;
79 for (std::vector<int>::const_iterator j = pTbl.begin();
80 j != pTbl.end();
81 j++, i++) {
82
83 // we processed already the i < *j case
84 if (i > *j) continue;
85
86 // unpaired position
87 if (*j == -1) {
88 structure[*j] = '.';
89 continue;
90 }
91
92 structure[i] = '(';
93 structure[*j] = ')';
94 }
95
96 return (structure);
97 }
98
99 std::vector<std::pair<int,int> >
PairTableToBasePairList(short int * pair_table)100 PairTableToBasePairList(short int * pair_table)
101 {
102 std::vector<std::pair<int,int> > pl;
103
104 //std::string structure(pair_table[0], '.');
105
106 int i=1;
107 for (int j = 1;j<= pair_table[0];j++, i++) {
108
109 int val=pair_table[j] ;
110 // we processed already the i < *j case
111 if (i > val) continue;
112
113 // unpaired position
114 if (val == 0) {
115 // structure[*j] = '.';
116 continue;
117 }
118 pl.push_back(std::make_pair(i,val));
119 //structure[i] = '(';
120 //structure[*j] = ')';
121 }
122
123 return pl;//(structure);
124 }
125
126
127 std::string
BasePairListToStructure(int length,std::vector<std::pair<int,int>> pl)128 BasePairListToStructure(int length, std::vector<std::pair<int,int> > pl)
129 {
130 std::string structure(length, '.');
131
132 for (std::vector<std::pair<int,int> >::const_iterator p = pl.begin();
133 p != pl.end();
134 p++) {
135
136 structure[(*p).first] = '(';
137 structure[(*p).second] = ')';
138 }
139
140 return (structure);
141 }
142
143 std::string
BasePairListToStructure1(int length,const std::vector<std::pair<int,int>> & pl)144 BasePairListToStructure1(int length, const std::vector<std::pair<int,int> > & pl)
145 //BasePairListToStructure1(int length, std::vector<std::pair<int,int> > pl)
146 {
147 std::string structure(length, '.');
148
149 //for(int i=0;i<pl.size();i++){
150 //structure[pl[i].first-1] = '(';
151 //structure[pl[i].second-1] = ')';
152 for (std::vector<std::pair<int,int> >::const_iterator p = pl.begin();
153 p != pl.end();
154 p++) {
155 structure[(*p).first-1] = '(';
156 structure[(*p).second-1] = ')';
157 }
158
159 return (structure);
160 }
161
162 /*
163 void
164 StringToTextFile(std::string sText, std::string filePath)
165 {
166 std::ofstream outfile(filePath.c_str());
167 outfile << sText;
168 }
169 */
170
171
IntroducesPseudoKnot(const std::vector<std::pair<int,int>> & node,const std::pair<int,int> & p1)172 bool IntroducesPseudoKnot(const std::vector<std::pair<int,int> >& node, const std::pair<int,int>& p1){
173 for(size_t i=0;i<node.size();i++){
174 std::pair<int,int> p2=node[i];
175 if(p1.first==p2.first || p1.first==p2.second || p1.second==p2.first || p1.second==p2.second) return true;
176 if(p1.first<=p2.first && p1.second>= p2.first && p1.second<=p2.second) return true;
177 if(p1.first>=p2.first && p1.first<=p2.second && p1.second>=p2.second) return true;
178 }
179 return false;
180 }
181
ConformationHasPair(const std::vector<std::pair<int,int>> & node,const std::pair<int,int> & p)182 bool ConformationHasPair(const std::vector<std::pair<int,int> >& node,const std::pair<int,int> & p){
183 for(size_t i=0;i<node.size();i++){
184 if(node[i].first==p.first && node[i].second==p.second) {
185 return true;
186 }
187 }
188 return false;
189 }
190
191
Conflict(const std::vector<std::pair<int,int>> & node,const std::pair<int,int> & p1)192 bool Conflict(const std::vector<std::pair<int,int> >& node, const std::pair<int,int>& p1){
193 if(ConformationHasPair(node,p1)) return true;
194 if(IntroducesPseudoKnot(node,p1)) return true;
195 return false;
196 }
197
Conflict(const std::vector<std::pair<int,int>> & node1,const std::vector<std::pair<int,int>> & node2)198 bool Conflict(const std::vector<std::pair<int,int> >& node1, const std::vector<std::pair<int,int> >& node2){
199 for(size_t i=0;i<node2.size();i++){
200 if(Conflict(node1,node2[i])) return true;
201 }
202 return false;
203
204 }
205
206 //std::vector<std::vector<std::pair<int,int> > >
207 void
ConformationToStacks(std::vector<std::vector<std::pair<int,int>>> & stacks,std::vector<std::pair<int,int>> node,int stacksize)208 ConformationToStacks(std::vector<std::vector<std::pair<int,int> > > & stacks, std::vector<std::pair<int,int> > node,int stacksize)
209 {
210
211 stacks.clear();
212
213
214 // return if node is empty
215 if (node.empty()) return;
216
217 sort(node.begin(), node.end());
218
219 //std::cout << "nodes" << std::endl
220 // << PrintBasePairList(node) << std::endl;
221
222 std::vector<std::pair<int,int> > stack;
223 copy(node.begin(), node.begin()+1, std::inserter(stack, stack.begin()));
224
225 //std::cout << "stack" << std::endl
226 // << PrintBasePairList(stack) << std::endl;
227
228 if (node.size()==1 && stacksize==1) {
229 stacks.push_back(stack);
230 return;
231 }
232 for(size_t i=1;i<node.size();i++){
233 bool added=false;
234 if(node[i].first==stack.back().first+1 &&
235 node[i].second==stack.back().second-1) {
236 added=true;
237 stack.push_back(std::pair<int,int>(stack.back().first+1,
238 stack.back().second-1));
239 }
240 else if(node[i].first==stack.back().first+2 &&
241 node[i].second==stack.back().second-1) {
242 added=true;
243 stack.push_back(std::pair<int,int>(stack.back().first+2,
244 stack.back().second-1));
245 }
246 else if(node[i].first==stack.back().first+1 &&
247 node[i].second==stack.back().second-2) {
248 added=true;
249 stack.push_back(std::pair<int,int>(stack.back().first+1,
250 stack.back().second-2));
251 }
252 if(stack.size()>=(size_t)stacksize && (!added || i==node.size()-1) ) {
253 stacks.push_back(stack);
254 stack.clear();
255 stack.push_back(node[i]);
256 }
257 if(!added) {
258 stack.clear();
259 stack.push_back(node[i]);
260 }
261 }
262 }
263
264
265 #ifdef WITH_DMALLOC
266 #define MG_space(S) calloc(1,(S))
267 #endif
268
269
MakePairTable(const char * structure)270 void MakePairTable(const char *structure)
271 {
272 /* returns array representation of structure.
273 table[i] is 0 if unpaired or j if (i.j) pair. */
274 short i,j,hx;
275 short length;
276 short *stack;
277 // short *table;
278
279 length = (short) strlen(structure);
280 stack = (short *) space(sizeof(short)*(length+1));
281 //table = (short *) space(sizeof(short)*(length+2));
282 pair_table[0] = length;
283
284 for (hx=0, i=1; i<=length; i++) {
285 switch (structure[i-1]) {
286 case '(':
287 stack[hx++]=i;
288 break;
289 case ')':
290 j = stack[--hx];
291 if (hx<0) {
292 // fprintf(stderr, "%s\n", structure);
293 std::cout<<structure<<" unbalanced brackets in Michael's MakePairTable"<<std::endl;
294 // nrerror("unbalanced brackets in Michael's MakePairTable");
295 }
296 pair_table[i]=j;
297 pair_table[j]=i;
298 break;
299 default: /* unpaired base, usually '.' */
300 pair_table[i]= 0;
301 break;
302 }
303 }
304 if (hx!=0) {
305 // fprintf(stderr, "%s\n", structure);
306 std::cout<<structure<<" unbalanced brackets in Michael's MakePairTable"<<std::endl;
307 //nrerror("unbalanced brackets in Michael's MakePairTable");
308 }
309 free(stack);
310 //return(table);
311 }
312
CreateRGBString(double r,double g,double b)313 std::string CreateRGBString(double r,double g,double b) {
314 return "["+Str(r)+" "+Str(g)+" "+Str(b)+"]";
315 }
316
317
318
319 /**
320 Code to create a .ps dotplot template. Based on PS_dot.c
321 */
322
PSFrontPlot(std::string sequence,std::vector<std::pair<int,int>> extrema)323 std::string PSFrontPlot(std::string sequence,std::vector<std::pair<int,int> > extrema){
324
325 std::string sText;
326 sText+=" %!PS-Adobe-3.0 EPSF-3.0\n";
327 sText+=" %%Title: Maximum Matching Color Plot\n";
328 sText+=" %%Creator: maxmatch.cpp\n";
329 sText+=" %%BoundingBox: 66 211 518 662\n";
330 sText+=" %%DocumentFonts: Helvetica\n";
331 sText+=" %%Pages: 1\n";
332 sText+=" %%EndComments\n";
333 sText+=" \n";
334 sText+=" %Options: \n";
335 sText+=" % \n";
336 sText+=" %Colored max matching matrixc.\n";
337 sText+=" % i j sqrt(p(i,j)) ubox\n";
338 sText+=" \n";
339 sText+=" %%BeginProlog\n";
340 sText+=" /DPdict 100 dict def\n";
341 sText+=" DPdict begin\n";
342 sText+=" /logscale false def\n";
343 sText+=" \n";
344 sText+=" \n";
345 sText+=" \n";
346 sText+=" /mylbox { % x y size [rgb] => -\n";
347 sText+=" exch 4 2 roll\n";
348 sText+=" len exch sub 1 add mybox\n";
349 sText+=" } bind def\n";
350 sText+=" \n";
351 sText+=" /mybox { % [rgb] size x y box - draws box centered on x,y\n";
352 sText+=" 2 index 0.5 mul add % x += 0.5\n";
353 sText+=" exch 2 index 0.5 mul add exch % x += 0.5\n";
354 sText+=" newpath\n";
355 sText+=" moveto\n";
356 sText+=" dup neg 0 rlineto\n";
357 sText+=" dup neg 0 exch rlineto\n";
358 sText+=" 0 rlineto\n";
359 sText+=" closepath\n";
360 sText+=" gsave\n";
361 sText+=" aload pop\n";
362 sText+=" setrgbcolor\n";
363 sText+=" fill\n";
364 sText+=" grestore\n";
365 sText+=" } bind def\n";
366 sText+=" \n";
367 sText+=" \n";
368 sText+=" \n";
369 sText+=" \n";
370 sText+=" /drawseq {\n";
371 sText+=" % print sequence along all 4 sides\n";
372 sText+=" [ \n";
373 // sText+=" [ [0.7 -0.3 0 ]\n"; remove bottom
374 sText+=" [0.7 0.7 len add 0]\n";
375 // sText+=" [-0.3 len sub -0.4 -90]\n"; remove left
376 sText+=" [-0.3 len sub 0.7 len add -90]\n";
377 sText+=" ] {\n";
378 sText+=" gsave\n";
379 sText+=" aload pop rotate translate\n";
380 sText+=" 0 1 len 1 sub {\n";
381 sText+=" dup 0 moveto\n";
382 sText+=" sequence exch 1 getinterval\n";
383 sText+=" show\n";
384 sText+=" } for\n";
385 sText+=" grestore\n";
386 sText+=" } forall\n";
387 sText+=" } bind def\n";
388 sText+=" \n";
389 sText+=" \n";
390 sText+=" /drawgrid{\n";
391 sText+=" 0.01 setlinewidth\n";
392 sText+=" len log 0.9 sub cvi 10 exch exp % grid spacing\n";
393 sText+=" dup 1 gt {\n";
394 sText+=" dup dup 20 div dup 2 array astore exch 40 div setdash\n";
395 sText+=" } { [0.3 0.7] 0.1 setdash } ifelse\n";
396 sText+=" 0 exch len {\n";
397 sText+=" dup dup\n";
398 sText+=" 0 moveto\n";
399 sText+=" len lineto \n";
400 sText+=" dup\n";
401 sText+=" len exch sub 0 exch moveto\n";
402 sText+=" len exch len exch sub lineto\n";
403 sText+=" stroke\n";
404 sText+=" } for\n";
405 sText+=" [] 0 setdash\n";
406 sText+=" 0.04 setlinewidth \n";
407 sText+=" currentdict /cutpoint known {\n";
408 sText+=" cutpoint 1 sub\n";
409 sText+=" dup dup -1 moveto len 1 add lineto\n";
410 sText+=" len exch sub dup\n";
411 sText+=" -1 exch moveto len 1 add exch lineto\n";
412 sText+=" stroke\n";
413 sText+=" } if\n";
414 sText+=" 0.5 neg dup translate\n";
415 sText+=" } bind def\n";
416 sText+=" \n";
417 sText+=" end\n";
418 sText+=" %%EndProlog\n";
419 sText+=" DPdict begin\n";
420 sText+=" %delete next line to get rid of title\n";
421 sText+=" %270 665 moveto /Helvetica findfont 14 scalefont setfont (dot.ps) show\n";
422 sText+=" %set up the matrix\n";
423 sText+=" /sequence { ("+sequence+") } def\n";
424 sText+=" /len { sequence length } bind def\n";
425 sText+=" 72 216 translate\n";
426 sText+=" 72 6 mul len 1 add div dup scale\n";
427 sText+=" /Helvetica findfont 0.95 scalefont setfont\n";
428 sText+=" \n";
429 sText+=" \n";
430 // sText+=" drawseq\n"; draw sequence after everything else has been drawn so it is not rendered invisible by colored squares
431 sText+=" 0.5 dup translate\n";
432 sText+=" % draw diagonal\n";
433 sText+=" 0.04 setlinewidth\n";
434 sText+=" 0 len moveto len 0 lineto stroke \n";
435
436 sText+=" %dynamic content\n";
437 sText+=" %drawgrid\n";
438 sText+=" %data starts here\n";
439 sText+=" % x y size rgb-array\n";
440
441 //Add the actual entries and their colours for the upper half of the matrix
442 for(std::vector<std::pair<int,int> >::iterator it=extrema.begin();it!=extrema.end();it++){
443 for(int i=it->first; i<=it->second; i++) {
444 for (int j=i; j<=it->second; j++) {
445
446 /*
447 double color_depth=ratio(i,j)/max_ratio;
448 //round to two decimals
449 color_depth=.01*(int)floor(color_depth*100+.5);
450 */
451 // string rgb=CreateRGBString(1.0,1.0,1.0);
452 // string rgb=CreateRGBString(0.5,0.5,0.5);
453 // string rgb=CreateRGBString(0.6,0.6,0.6);
454 std::string rgb=CreateRGBString(0.7,0.7,0.7);
455 //string rgb=CreateRGBString(0.8,0.8,0.8);
456 sText+=Str(j)+"\t"+Str(i)+"\t"+Str(1)+"\t"+rgb+"\t"+"mylbox\n";
457
458 }
459 }
460 }
461 sText+=" drawseq\n";
462 //draw a colour scale in the lower half of the matrix
463 // for(int i=0;i<10;i++) {
464 // string rgb=CreateRGBString(.1*i);
465 // sText+=Str(i+1)+"\t"+Str(l-1)+"\t"+Str(1)+"\t"+rgb+"\t"+"mylbox\n";
466 //}
467
468 sText+=" showpage\n";
469 return sText;
470
471 }
472