1 #include <unistd.h>
2
3 #include "aligner.h"
4 #include "config.h"
5
printCluster(std::pair<double,RNAProfileAlignment * > cluster,unsigned int clusterNr,double minPairProb,const Options & options)6 void printCluster(std::pair<double, RNAProfileAlignment*> cluster,
7 unsigned int clusterNr, double minPairProb, const Options & options) {
8 std::cout << "RNA Structure Cluster Nr: " << clusterNr << std::endl;
9 std::cout << "Score: " << cluster.first << std::endl;
10 std::cout << "Members: " << cluster.second->getNumStructures() << std::endl << std::endl;
11 if (options.has(Options::FastaOutput)) {
12 cluster.second->printFastaAli(false);
13 std::cout << std::endl;
14 } else {
15
16 if (not options.has(Options::ShowOnlyScore)){
17 cluster.second->printSeqAli();
18
19 #ifdef HAVE_LIBRNA
20 // print alignment
21 cluster.second->printStrAli();
22 std::cout << std::endl;
23 #endif
24 }
25 }
26
27 // save profile
28 if (options.has(Options::SaveProfile)) {
29 std::string filename;
30 filename = options.generateFilename(Options::SaveProfile,".pro", "rna.pro",clusterNr);
31 cluster.second->save(filename);
32 }
33
34
35 if (not options.has(Options::ShowOnlyScore)){
36 // print consensus structure
37 std::cout << "Consensus sequence/structure:" << std::endl;
38 cluster.second->printConsensus(minPairProb);
39 }
40 std::cout << std::endl;
41
42 #ifdef HAVE_LIBG2
43 #ifdef HAVE_LIBRNA
44 // generate squiggle plots
45 if (options.has(Options::MakeSquigglePlot)) {
46 RNAProfileAlignment::SquigglePlotOptions sqOptions;
47 std::string filename;
48
49 // plot showing full base information
50 filename = options.generateFilename(Options::MakeSquigglePlot,".ps", "rnaprofile.ps",clusterNr);
51 sqOptions.greyColors = options.has(Options::SquiggleGreyColors);
52 sqOptions.minPairProb = minPairProb;
53 sqOptions.mostLikelySequence = false;
54 cluster.second->squigglePlot(filename,sqOptions);
55
56 // plot showing consensus sequence
57 filename = options.generateFilename(Options::MakeSquigglePlot,"_cons.ps", "rnaprofile_cs.ps",clusterNr);
58 sqOptions.mostLikelySequence = true;
59 cluster.second->squigglePlot(filename,sqOptions);
60 }
61 #endif
62 #endif
63 }
64
65
66
alignMultiple(std::vector<RNAProfileAlignment * > & inputListMult,Score & score,const Options & options,bool anchoring,RNAFuncs::AddXmlInfos & xmlInfos)67 void alignMultiple(std::vector<RNAProfileAlignment*> &inputListMult, Score &score,const Options &options,bool anchoring,RNAFuncs::AddXmlInfos &xmlInfos) {
68
69 std::vector<std::pair<double,RNAProfileAlignment*> > resultList;
70 progressiveAlign(inputListMult, resultList, score, options, anchoring);
71
72 double minPairProb;
73 options.get(Options::ConsensusMinPairProb, minPairProb, 0.5);
74
75 std::cout << std::endl;
76 std::cout << std::endl;
77 std::cout << "*** Results ***" << std::endl << std::endl;
78 std::cout << "Minimum basepair probability for consensus structure (-cmin): " << minPairProb << std::endl << std::endl;
79
80 unsigned int clusterNr = 1;
81 std::vector<std::pair<double,RNAProfileAlignment*> >::const_iterator cluster;
82 for (cluster=resultList.begin(); cluster!=resultList.end(); cluster++) {
83 printCluster(*cluster, clusterNr, minPairProb, options);
84 clusterNr++;
85 }
86 #ifdef HAVE_LIBXMLPLUSPLUS
87 #ifdef HAVE_LIBXML2
88 // generate xml
89 std::string outputFile;
90 if (options.has(Options::XmlOutputFile)) {
91 options.get(Options::XmlOutputFile,outputFile,std::string(""));
92 } else {
93 outputFile = "result.xml";
94 }
95
96 if (options.has(Options::GenerateXML)) {
97 RNAFuncs::printMAliXML(resultList,options,minPairProb,xmlInfos,outputFile);
98 }
99 #endif
100 #endif
101
102 /* TODO das geht alle noch nicht, geht aber in pairwise
103 if(!options.has(Options::ShowOnlyScore))
104 {
105 std::cout << "global optimal score: ";
106 }
107 std::cout << optScore << std::endl;
108
109 vector<RNAProfileAlignment*>::iterator it=inputMapProfile.begin();
110 RNAProfileAlignment *f=it;
111 if(!options.has(Options::ShowOnlyScore))
112 {
113 // generate dot file
114 makeDotFileAli(*f,options);
115
116 f->print();
117 std::cout << std::endl;
118 f->printConsensus();
119 }
120 */
121
122 // save profile alignment to binary file
123 /* if(options.has(Options::SaveMultipleAliFile))
124 {
125 std::string filename;
126
127 filename=generateFilename(options,Options::SaveMultipleAliFile,".mta", "unknown.dot");
128 std::ofstream s(filename.c_str());
129 f1->save(s);
130 }
131 */
132
133 /*
134 // generate squiggle plots
135 if(options.has(Options::MakeSquigglePlot))
136 {
137 RNAProfileAlignment::SquigglePlotOptions sqOptions;
138 std::string filename;
139
140 // plot showing full base information
141 filename=options.generateFilename(Options::MakeSquigglePlot,".ps", "rnaprofile.ps");
142 sqOptions.greyColors=options.has(Options::SquiggleGreyColors);
143 sqOptions.mostLikelySequence=false;
144 f->squigglePlot(filename,sqOptions);
145
146 // plot showing consensus sequence
147 filename=options.generateFilename(Options::MakeSquigglePlot,"_cons.ps", "rnaprofile_cons.ps");
148 sqOptions.greyColors=options.has(Options::SquiggleGreyColors);
149 sqOptions.mostLikelySequence=true;
150 f->squigglePlot(filename,sqOptions);
151 }
152 */
153
154
155 }
156
157 //void computeSpaceTimeInfo(const RNAForest* f1, const RNAForest* f2, const Score &score, const Options &options) {
158 //
159 // tms tmsStart, tmsEnd;
160 // Algebra<double,RNA_Alphabet> *algGlobClassic = new DoubleDist_Algebra(score);
161 //
162 // if (!options.has(Options::ShowOnlyScore)) {
163 // std::cout << "F1_NUMNODES" << ";";
164 // std::cout << "F2_NUMNODES" << ";";
165 // std::cout << "F1_DEGREE" << ";";
166 // std::cout << "F2_DEGREE" << ";";
167 // std::cout << "F1_LEAVES" << ";";
168 // std::cout << "F2_LEAVES" << ";";
169 // std::cout << "F1_DEPTH" << ";";
170 // std::cout << "F2_DEPTH" << ";";
171 // std::cout << "F1_NUMCSFS" << ";";
172 // std::cout << "F2_NUMCSFS" << ";";
173 // std::cout << "TABLE_SIZE_4D" << ";";
174 // std::cout << "TABLE_SIZE_2D" << ";";
175 // std::cout << "GLOBALI_TIME" << ";";
176 // std::cout << "GLOBALI_TIME_SPEEDUP" << ";";
177 // std::cout << "LOCALALI_TIME" << ";";
178 // std::cout << "LOCALALI_TIME_SPEEDUP" << ";";
179 // std::cout << std::endl;
180 // }
181 //
182 // std::cout << f1->size() << "\t";
183 // std::cout << f2->size() << "\t";
184 // std::cout << f1->maxDegree() << "\t";
185 // std::cout << f2->maxDegree() << "\t";
186 // std::cout << f1->numLeaves() << "\t";
187 // std::cout << f2->numLeaves() << "\t";
188 // std::cout << f1->maxDepth() << "\t";
189 // std::cout << f2->maxDepth() << "\t";
190 // std::cout << f1->getNumCSFs() << "\t";
191 // std::cout << f2->getNumCSFs() << "\t";
192 // std::cout << f1->size()*f2->size()*f1->maxDegree()*f1->maxDegree() << "\t";
193 // std::cout << f1->getNumCSFs()*f2->getNumCSFs() << "\t";
194 //
195 // bool topdown = false;
196 // bool anchored = options.has(Options::Anchoring);
197 // // global alignment
198 // {
199 // times(&tmsStart);
200 // AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,topdown,anchored,false,false);
201 // times(&tmsEnd);
202 // std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t";
203 // }
204 // // global alignment speedup
205 // {
206 // times(&tmsStart);
207 // AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,topdown,anchored,false,true);
208 // times(&tmsEnd);
209 // std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t";
210 // }
211 // // local alignment
212 // {
213 // times(&tmsStart);
214 // AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,topdown,anchored,true,false);
215 // times(&tmsEnd);
216 // std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t";
217 // }
218 // // local alignment speedup
219 // {
220 // times(&tmsStart);
221 // AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,topdown,anchored,true,true);
222 // times(&tmsEnd);
223 // std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t";
224 // }
225 // std::cout << std::endl;
226 //}
227
computeSpaceTimeInfo(const RNAForest * f1,const RNAForest * f2,const Score & score,const Options & options)228 void computeSpaceTimeInfo(const RNAForest* f1, const RNAForest* f2, const Score &score, const Options &options) {
229
230 // tms tmsStart, tmsEnd;
231 // Algebra<double,RNA_Alphabet> *algGlobClassic = new DoubleDist_Algebra(score);
232 if (!options.has(Options::ShowOnlyScore)) {
233 std::cout << "F1_NUMNODES" << ";";
234 std::cout << "F2_NUMNODES" << ";";
235 std::cout << "F1_DEGREE" << ";";
236 std::cout << "F2_DEGREE" << ";";
237 std::cout << "F1_LEAVES" << ";";
238 std::cout << "F2_LEAVES" << ";";
239 std::cout << "F1_DEPTH" << ";";
240 std::cout << "F2_DEPTH" << ";";
241 std::cout << "F1_NUMCSFS" << ";";
242 std::cout << "F2_NUMCSFS" << ";";
243 std::cout << "TABLE_SIZE_4D" << ";";
244 std::cout << "TABLE_SIZE_2D" << ";";
245
246 for (int topdown=0; topdown<=1; topdown++) {
247 for (int anchored=0; anchored<=1; anchored++) {
248 for (int local=0; local<=1; local++) {
249 for (int affine=0; affine <=1; affine++) {
250 if ((anchored && !options.has(Options::Anchoring))||(anchored && !topdown))
251 continue;
252 if (affine)
253 std::cout << "AFFINE_";
254 else
255 std::cout << "LINEAR_";
256 if (topdown)
257 std::cout << "TOPDOWN_";
258 else
259 std::cout << "BOTTOM_UP_";
260 if (anchored)
261 std::cout << "ANCHORED_";
262 else
263 std::cout << "NONANCHORED_";
264 if (local)
265 std::cout << "LOCAL" << ";";
266 else
267 std::cout << "GLOBAL" << ";";
268 }
269 }
270 }
271 }
272
273 std::cout << std::endl;
274 }
275 tms tmsStart, tmsEnd;
276 Algebra<double,RNA_Alphabet> *alg = new DoubleDist_Algebra(score);
277 AlgebraAffine<double,RNA_Alphabet> *alg_affine = new AffineDoubleDist_Algebra(score);
278
279 std::cout << f1->size() << "\t";
280 std::cout << f2->size() << "\t";
281 std::cout << f1->maxDegree() << "\t";
282 std::cout << f2->maxDegree() << "\t";
283 std::cout << f1->numLeaves() << "\t";
284 std::cout << f2->numLeaves() << "\t";
285 std::cout << f1->maxDepth() << "\t";
286 std::cout << f2->maxDepth() << "\t";
287 std::cout << f1->getNumCSFs() << "\t";
288 std::cout << f2->getNumCSFs() << "\t";
289 std::cout << f1->size()*f2->size()*f1->maxDegree()*f1->maxDegree() << "\t";
290 std::cout << f1->getNumCSFs()*f2->getNumCSFs() << "\t";
291
292 bool printBT = true; // never enabled
293 bool speedup = true; // always enabled
294
295 for (int topdown=0; topdown<=1; topdown++) {
296 for (int anchored=0; anchored<=1; anchored++) {
297 for (int local=0; local<=1; local++) {
298 for (int affine=0; affine <=1; affine++) {
299
300 // TODO man muesste eig abpruefen, ob der input fuer's anchoring ausreicht, nicht, ob es option ist
301 if ((anchored && !options.has(Options::Anchoring))||(anchored && !topdown))
302 continue;
303
304 std::cout << "\t" << std::flush;
305 //std::cout << std::endl;
306 //if (affine)
307 // std::cout << "affine ";
308 //else
309 // std::cout << "linear ";
310 //if (topdown)
311 // std::cout << "topdown ";
312 //else
313 // std::cout << "bottomup ";
314 //if (anchored)
315 // std::cout << "anchored ";
316 //else
317 // std::cout << "nonanchored ";
318 //if (local)
319 // std::cout << "local ";
320 //else
321 // std::cout << "global ";
322 //std::cout << std::endl;
323
324 if (affine) {
325 times(&tmsStart);
326 AlignmentAffine<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg_affine,topdown,anchored,local,!printBT,speedup);
327 times(&tmsEnd);
328 std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t" << std::flush;
329 }
330 else {
331 times(&tmsStart);
332 AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg,topdown,anchored,local,!printBT,speedup);
333 times(&tmsEnd);
334 std::cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t" << std::flush;
335 }
336
337 }
338 }
339 }
340 }
341 std::cout << std::endl;
342 std::cout << std::endl;
343 std::cout << std::endl;
344 }
345
alignPairwise(std::vector<RNAForest * > & inputListPW,Score & score,const Options & options,bool anchored,RNAFuncs::AddXmlInfos & xmlInfos)346 void alignPairwise(std::vector<RNAForest*> &inputListPW, Score &score, const Options &options, bool anchored, RNAFuncs::AddXmlInfos &xmlInfos) {
347
348 #ifdef HAVE_LIBG2
349 RNAFuncs::SquigglePlotOptions sqOptions;
350 // generate squiggle plot
351 if (options.has(Options::MakeSquigglePlot)) {
352 // get sq options
353 sqOptions.hideBaseNumbers = options.has(Options::SquiggleHideBaseNumbers);
354 options.get(Options::SquiggleBaseNumberInterval,sqOptions.baseNumInterval,(unsigned int)20);
355 sqOptions.greyColors = options.has(Options::SquiggleGreyColors);
356 options.get(Options::SquiggleScaleFactor,sqOptions.scale,1.0);
357 sqOptions.generateFIG = options.has(Options::SquiggleGenerateFIG);
358 #ifdef HAVE_LIBGD
359 sqOptions.generatePNG = options.has(Options::SquiggleGeneratePNG);
360 sqOptions.generateJPG = options.has(Options::SquiggleGenerateJPG);
361 #endif
362
363 }
364 #endif
365
366 RNAForest *f1 = inputListPW.front();
367 RNAForest *f2 = inputListPW.back();
368
369 if (options.has(Options::SpaceTimeInfo)) {
370 computeSpaceTimeInfo(f1, f2, score, options);
371 exit(EXIT_SUCCESS);
372 }
373
374 // for global we compute less entries in the matrix
375 bool topdown = options.has(Options::Topdown);
376 //bool anchored = options.has(Options::Anchoring);
377 bool local = options.has(Options::LocalSimilarity);
378 bool printBT = options.has(Options::Backtrace);
379
380 RNA_Algebra<double,RNA_Alphabet> *alg = NULL;
381 RNA_AlgebraAffine<double,RNA_Alphabet> * alg_affine = NULL;
382 Alignment<double,RNA_Alphabet,RNA_AlphaPair> * ali = NULL;
383
384 if (options.has(Options::Affine)) {
385 if (options.has(Options::CalculateDistance))
386 alg_affine = new AffineDoubleDistRNA_Algebra(score);
387 else if (options.has(Options::RIBOSUMScore))
388 alg_affine = new AffineRIBOSUM8560(score);
389 else
390 alg_affine = new AffineDoubleSimiRNA_Algebra(score);
391 ali = new AlignmentAffine<double,RNA_Alphabet,RNA_AlphaPair>(f1,f2,*alg_affine,topdown,anchored,local,printBT);
392 }
393 else {
394 if (options.has(Options::CalculateDistance))
395 alg = new DoubleDistRNA_Algebra(score);
396 else if (options.has(Options::RIBOSUMScore))
397 alg = new RIBOSUM8560(score);
398 else
399 alg = new DoubleSimiRNA_Algebra(score);
400
401 ali = new AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair>(f1,f2,*alg,topdown,anchored,local,printBT);
402 }
403
404 if (options.has(Options::Tables)) { // TODO mit stringstream zusammenbauen
405 if (options.has(Options::Topdown)) {
406 if (options.has(Options::Affine))
407 ali->printToFile(std::string("tables_affine_topdown.txt"));
408 else
409 ali->printToFile(std::string("tables_linear_topdown.txt"));
410 }
411 else {
412 if (options.has(Options::Affine))
413 ali->printToFile(std::string("tables_affine.txt"));
414 else
415 ali->printToFile(std::string("tables_linear.txt"));
416 }
417 }
418
419 RNA_Alignment fali;
420 fali.setStructureNames(f1->getName(),f2->getName());
421
422 #if defined(HAVE_LIBRNA) && defined(HAVE_LIBXMLPLUSPLUS) && defined(HAVE_LIBXML2)
423 double optScore = getAndPrintOptimum(ali, options);
424 #else
425 getAndPrintOptimum(ali, options);
426 #endif
427
428 std::vector<std::pair<unsigned int,unsigned int> > xsubopts;
429 std::vector<std::pair<unsigned int,unsigned int> > ysubopts;
430 if (!options.has(Options::ShowOnlyScore)) {
431
432 // read options for local simi
433 int suboptPercent;
434 options.get(Options::LocalSubopts,suboptPercent,100);
435 // for local simi
436 unsigned int xbasepos,ybasepos,xlen,ylen;
437 if (options.has(Options::SmallInLarge)) {
438 ali->getOptSILAlignment(fali,ybasepos);
439 std::cout << "starting at position: " << ybasepos << std::endl << std::endl;
440 }
441 else if (options.has(Options::LocalSimilarity)) {
442 ali->resetOptLocalAlignment(suboptPercent);
443 ali->getOptLocalAlignment(fali,xbasepos,ybasepos);
444 std::cout << "starting at positions: " << xbasepos << "," << ybasepos << std::endl << std::endl;
445 xmlInfos.xbasepos = xbasepos;
446 xmlInfos.ybasepos = ybasepos;
447 }
448 else
449 ali->getOptGlobalAlignment(fali);
450
451 // generate dot file
452 makeDotFileAli(fali,options);
453
454 // print alignment
455 //printRNAAlignment(RNA_Alignment fali, Options options);
456 std::string seq1,seq2,str1,str2;
457 fali.getSequenceAlignments(seq1,seq2);
458 fali.getStructureAlignment(str1,true);
459 fali.getStructureAlignment(str2,false);
460
461 if (options.has(Options::FastaOutput)) {
462 std::cout << fali.getStructureNameX() << std::endl;
463 std::cout << seq1 << std::endl;
464 std::cout << str1 << std::endl;
465 std::cout << fali.getStructureNameY() << std::endl;
466 std::cout << seq2 << std::endl;
467 std::cout << str2 << std::endl;
468 std::cout << std::endl;
469 }
470 else {
471 RNAFuncs::printAli(fali.getStructureNameX(),fali.getStructureNameY(),seq1,seq2,str1,str2);
472 }
473
474 xlen = seq1.size();
475 ylen = seq2.size();
476 if (options.has(Options::LocalSimilarity)) {
477 xsubopts.push_back(std::pair<unsigned int,unsigned int>(xbasepos,xlen));
478 ysubopts.push_back(std::pair<unsigned int,unsigned int>(ybasepos,ylen));
479 }
480
481 // for squiggle plot and subopts squiggle plot
482 unsigned int count = 1;
483
484 #ifdef HAVE_LIBG2
485 #ifdef HAVE_LIBRNA
486 char s[8];
487 if (options.has(Options::MakeSquigglePlot)) {
488 sprintf(s,"%d",count);
489 fali.squigglePlot(s,sqOptions);
490 }
491 #endif
492 #endif
493
494 // suboptimal alignments
495 if (options.has(Options::LocalSubopts)) {
496 while (ali->nextLocalSuboptimum()) {
497 count++;
498 std::cout << "local optimal score: ";
499 std::cout << ali->getLocalOptimum() << std::endl;
500 RNA_Alignment fali;
501 unsigned int xbasepos = 10, ybasepos = 10;
502 ali->getOptLocalAlignment(fali,xbasepos,ybasepos);
503 std::cout << "starting at positions: " << xbasepos << "," << ybasepos << std::endl << std::endl;
504
505 // ab hier und das pushback herausnehmen - fasta option abfragen (print subopt)
506 // print alignment
507 fali.getSequenceAlignments(seq1,seq2);
508 fali.getStructureAlignment(str1,true);
509 fali.getStructureAlignment(str2,false);
510 RNAFuncs::printAli(fali.getStructureNameX(),fali.getStructureNameY(),seq1,seq2,str1,str2);
511 xlen = seq1.size();
512 ylen = seq2.size();
513 xsubopts.push_back(std::pair<unsigned int,unsigned int>(xbasepos,xlen));
514 ysubopts.push_back(std::pair<unsigned int,unsigned int>(ybasepos,ylen));
515
516 #ifdef HAVE_LIBG2
517 #ifdef HAVE_LIBRNA
518 if (options.has(Options::MakeSquigglePlot)) {
519 sprintf(s,"%d",count);
520 fali.squigglePlot(s,sqOptions);
521 }
522 #endif
523 #endif
524 }
525 }
526 }
527
528 // generateXML(Options & options, RNA_Alignment & fali, seq1, seq2, str1, str2, optScore, xmlInfos)
529 // NOTE der folgende Teil druckt bei den subopts nur das optimale
530 #ifdef HAVE_LIBRNA
531 #ifdef HAVE_LIBXMLPLUSPLUS
532 #ifdef HAVE_LIBXML2
533 // generate xml
534 std::string outputFile;
535 if (options.has(Options::XmlOutputFile)) {
536 options.get(Options::XmlOutputFile,outputFile,std::string(""));
537 } else {
538 outputFile = "result.xml";
539 }
540 if (options.has(Options::GenerateXML)) {
541 RNAFuncs::printPAliXML(fali.getStructureNameX(),fali.getStructureNameY(),seq1,seq2,str1,str2,optScore,options,xmlInfos,outputFile);
542 }
543 #endif
544 #endif
545 #endif
546
547 // der plot2d bekommt alle opt und subopt. alis zum ausdruck
548 #ifdef HAVE_LIBG2
549 #ifdef HAVE_LIBRNA
550 // generate squiggle plot
551 if (options.has(Options::MakeSquigglePlot)) {
552 f1->plot2d("x_str", xsubopts, sqOptions);
553 f2->plot2d("y_str", ysubopts, sqOptions);
554 }
555 #endif
556 #endif
557
558 // clear input vector
559 //std::vector<RNAForest*>::const_iterator it;
560 //for (it = inputListPW.begin(); it!=inputListPW.end(); it++)
561 // delete *it;
562 inputListPW.clear();
563 if (options.has(Options::Affine))
564 delete alg_affine;
565 else
566 delete alg;
567 }
568
569
570
editPairwise(std::vector<RNAForestSZ * > & inputListSZ,Score & score,Options & options,bool anchored)571 void editPairwise(std::vector<RNAForestSZ*> &inputListSZ,Score &score,Options &options, bool anchored) {
572 // timeb t1,t2;
573 SZAlgebra<double,RNA_Alphabet> *alg;
574
575 // if(options.has(Options::CalculateDistance))
576 //alg=new DoubleDistSZAlgebra(score);
577 alg=new IntDistSZAlgebra(score);
578 // else
579 // alg=new DoubleSimiSZAlgebra(score);
580
581 RNAForestSZ *f1=inputListSZ.front();
582 RNAForestSZ *f2=inputListSZ.back();
583
584 // ftime(&t1);
585 Mapping<double,RNA_Alphabet> mapping(f1,f2,*alg);
586 // ftime(&t2);
587
588 std::cout << "Global optimum: " << mapping.getGlobalOptimum() << std::endl;
589 // std::cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << std::endl;
590 }
591
592
alignPairwiseSimple(std::vector<RNAForest * > & inputListPW,Score & score,Options & options,bool anchored)593 void alignPairwiseSimple(std::vector<RNAForest*> &inputListPW,Score &score,Options &options, bool anchored) {
594 // timeb t1,t2;
595 Algebra<double,RNA_Alphabet> *alg;
596
597 // if(options.has(Options::CalculateDistance))
598 alg = new DoubleDist_Algebra(score);
599 // else
600 // alg=new ScoreAlgebraSimple(score);
601
602 RNAForest *f1 = inputListPW.front();
603 RNAForest *f2 = inputListPW.back();
604
605 // ftime(&t1);
606 // no local, no print backtrace, but speedup
607 AlignmentLinear<double,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg,false,false,false,true);
608 // ftime(&t2);
609
610 std::cout << "Global optimum: " << ali.getGlobalOptimum() << std::endl;
611 // std::cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << std::endl;
612 }
613
getAndPrintOptimum(Alignment<double,RNA_Alphabet,RNA_AlphaPair> * ali,const Options & options)614 double getAndPrintOptimum(Alignment<double,RNA_Alphabet,RNA_AlphaPair>* ali, const Options & options) {
615 if (!options.has(Options::ShowOnlyScore)) {
616 if (options.has(Options::SmallInLarge))
617 std::cout << "small-in-large optimal score: ";
618 else if (options.has(Options::LocalSimilarity))
619 std::cout << "local optimal score: ";
620 else
621 std::cout << "global optimal score: ";
622 }
623
624 if (options.has(Options::SmallInLarge)) {
625 std::cout << ali->getSILOptimum() << std::endl;
626 return ali->getSILOptimum();
627 }
628 else if (options.has(Options::LocalSimilarity)) {
629 std::cout << ali->getLocalOptimum() << std::endl;
630 return ali->getLocalOptimum();
631 }
632 else {
633 std::cout << ali->getGlobalOptimum() << std::endl;
634 if (options.has(Options::RelativeScore))
635 std::cout << ali->getGlobalOptimumRelative() << std::endl;
636 return ali->getGlobalOptimum();
637 }
638 }
639
640
641