1 /***************************************************************************
2 * Copyright (C) 2014 by *
3 * Lam-Tung Nguyen <nltung@gmail.com> *
4 * *
5 * *
6 * This program is free software; you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation; either version 2 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program; if not, write to the *
18 * Free Software Foundation, Inc., *
19 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
20 ***************************************************************************/
21
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <cmath>
26 #include <assert.h>
27
28 #define GLOBAL_VARIABLES_DEFINITION
29
30 #if !defined WIN32 && !defined _WIN32 && !defined __WIN32__
31 #include <sys/resource.h>
32 #endif
33
34 #include "tree/phylotree.h"
35 #include "pllnni.h"
36 #include "alignment/alignment.h"
37
38 /* program options */
39 int nni0;
40 int nni5;
41 extern VerboseMode verbose_mode;
42 int NNI_MAX_NR_STEP = 10;
43
44 /* program options */
45 extern Params *globalParams;
46 extern Alignment *globalAlignment;
47
48 /**
49 * map from newick tree string to frequencies that a tree is revisited during tree search
50 */
51 StringIntMap pllTreeCounter;
52
53
54 /*
55 * ****************************************************************************
56 * pllUFBoot area
57 * ****************************************************************************
58 */
59
60 pllUFBootData * pllUFBootDataPtr = NULL;
61
62
compareDouble(const void * a,const void * b)63 int compareDouble(const void * a, const void * b) {
64 if (*(double*) a > *(double*) b)
65 return 1;
66 else if (*(double*) a < *(double*) b)
67 return -1;
68 else
69 return 0;
70 }
71
getNNIList(pllInstance * tr)72 pllNNIMove *getNNIList(pllInstance* tr) {
73 static pllNNIMove* nniList;
74 if (nniList == NULL) {
75 nniList = (pllNNIMove*) malloc(2 * (tr->mxtips - 3) * sizeof(pllNNIMove));
76 ASSERT(nniList != NULL);
77 }
78 return nniList;
79 }
80
getNonConflictNNIList(pllInstance * tr)81 pllNNIMove *getNonConflictNNIList(pllInstance* tr) {
82 static pllNNIMove* nonConfNNIList;
83 if (nonConfNNIList == NULL) {
84 nonConfNNIList = (pllNNIMove*) malloc((tr->mxtips - 3) * sizeof(pllNNIMove));
85 ASSERT(nonConfNNIList != NULL);
86 }
87 return nonConfNNIList;
88 }
89
pllDoRandomNNIs(pllInstance * tr,partitionList * pr,int numNNI)90 double pllDoRandomNNIs(pllInstance *tr, partitionList *pr, int numNNI) {
91 int numInBrans = tr->mxtips - 3;
92 int numNNIinStep = (int) numInBrans / 5;
93
94 // devided in multiple round, each round collect 1/5 of numNNI
95 int cnt1 = 0;
96 unordered_set<int> selectedNodes;
97 vector<nodeptr> selectedBrans;
98 vector<nodeptr> brans;
99 do {
100 int cnt2 = 0;
101 selectedNodes.clear();
102 selectedBrans.clear();
103 brans.clear();
104 pllGetAllInBran(tr, brans);
105 ASSERT(brans.size() == numInBrans);
106 while(cnt2 < numNNIinStep && cnt2 < numNNI) {
107 int branIndex = random_int(numInBrans);
108 if (selectedNodes.find(brans[branIndex]->number) == selectedNodes.end() &&
109 selectedNodes.find(brans[branIndex]->back->number) == selectedNodes.end()) {
110 selectedNodes.insert(brans[branIndex]->number);
111 selectedNodes.insert(brans[branIndex]->back->number);
112 selectedBrans.push_back(brans[branIndex]);
113 cnt2++;
114 }
115 }
116 for (vector<nodeptr>::iterator it = selectedBrans.begin(); it != selectedBrans.end(); ++it) {
117 int nniType = random_int(2);
118 doOneNNI(tr, pr, (*it), nniType, TOPO_ONLY);
119 }
120 cnt1 += selectedBrans.size();
121 if (numNNI - cnt1 < numNNIinStep) {
122 numNNIinStep = numNNI - cnt1;
123 }
124 } while (cnt1 < numNNI);
125 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
126 pllOptimizeBranchLengths(tr, pr, 1);
127 return tr->likelihood;
128 }
129
pllGetAllInBran(pllInstance * tr,vector<nodeptr> & branlist)130 void pllGetAllInBran(pllInstance *tr, vector<nodeptr> &branlist) {
131 nodeptr p = tr->start->back;
132 nodeptr q = p->next;
133 while (q != p) {
134 pllGetAllInBranForSubtree(tr, q->back, branlist);
135 q = q->next;
136 }
137 }
138
pllGetAllInBranForSubtree(pllInstance * tr,nodeptr p,vector<nodeptr> & branlist)139 void pllGetAllInBranForSubtree(pllInstance *tr, nodeptr p, vector<nodeptr> &branlist) {
140 if (!isTip(p->number, tr->mxtips) && !isTip(p->back->number, tr->mxtips)) {
141 branlist.push_back(p);
142 nodeptr q = p->next;
143 while (q != p) {
144 pllGetAllInBranForSubtree(tr, q->back, branlist);
145 q = q->next;
146 }
147 }
148 }
149
pllPerturbTree(pllInstance * tr,partitionList * pr,vector<pllNNIMove> & nnis)150 double pllPerturbTree(pllInstance *tr, partitionList *pr, vector<pllNNIMove> &nnis) {
151 //printf("Perturbing %d NNIs \n", numNNI);
152 for (vector<pllNNIMove>::iterator it = nnis.begin(); it != nnis.end(); it++) {
153 printf("Do pertubing NNI (%d - %d) with logl = %10.4f \n", (*it).p->number, (*it).p->back->number, (*it).likelihood);
154 doOneNNI(tr, pr, (*it).p, (*it).nniType, TOPO_ONLY);
155 updateBranchLengthForNNI(tr, pr, (*it));
156
157 }
158 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
159 pllOptimizeBranchLengths(tr, pr, 1);
160 return tr->likelihood;
161 }
162
quicksort_nni(pllNNIMove * arr,int left,int right)163 void quicksort_nni(pllNNIMove* arr,int left, int right) {
164 int i = left, j = right;
165 pllNNIMove tmp, pivot = arr[(left + right) / 2];
166
167 /* partition */
168 while (i <= j) {
169 while (arr[i].likelihood < pivot.likelihood)
170 i++;
171 while (pivot.likelihood < arr[j].likelihood)
172 j--;
173 if (i <= j) {
174 tmp = arr[i];
175 arr[i] = arr[j];
176 arr[j] = tmp;
177
178 i++;
179 j--;
180 }
181 };
182
183 /* recursion */
184 if (left < j)
185 quicksort_nni(arr, left, j);
186 if (i < right)
187 quicksort_nni(arr, i, right);
188 }
189
190 //TODO: Workaround for memory leak problem when calling setupTopol within doNNISearch
_setupTopol(pllInstance * tr)191 topol *_setupTopol(pllInstance* tr) {
192 static topol* tree;
193 if (tree == NULL)
194 tree = setupTopol(tr->mxtips);
195 return tree;
196 }
197
getAffectedBranches(pllInstance * tr,nodeptr p)198 vector<string> getAffectedBranches(pllInstance* tr, nodeptr p) {
199 vector<string> res;
200 res.push_back(getBranString(p));
201 nodeptr q = p->back;
202 nodeptr p_nei = p->next;
203 nodeptr q_nei = q->next;
204 while (p_nei != p) {
205 res.push_back(getBranString(p_nei));
206 if (!isTip(p_nei->back->number, tr->mxtips)) {
207 res.push_back(getBranString(p_nei->back->next));
208 res.push_back(getBranString(p_nei->back->next->next));
209 }
210 p_nei = p_nei->next;
211 }
212 while (q_nei != q) {
213 res.push_back(getBranString(q_nei));
214 if (!isTip(q_nei->back->number, tr->mxtips)) {
215 res.push_back(getBranString(q_nei->back->next));
216 res.push_back(getBranString(q_nei->back->next->next));
217 }
218 q_nei = q_nei->next;
219 }
220 return res;
221 }
222
getBranString(nodeptr p)223 string getBranString(nodeptr p) {
224 stringstream branString;
225 nodeptr q = p->back;
226 if (p->number < q->number) {
227 branString << p->number << "-" << q->number;
228 } else {
229 branString << q->number << "-" << p->number;
230 }
231 return branString.str();
232 }
233
getAffectedNodes(pllInstance * tr,nodeptr p)234 set<int> getAffectedNodes(pllInstance* tr, nodeptr p) {
235 set<int> nodeSet;
236 nodeptr q = p->back;
237 nodeptr p_nei = p->next;
238 nodeptr q_nei = q->next;
239 nodeSet.insert(p->number);
240 nodeSet.insert(q->number);
241 while (p_nei != p) {
242 nodeptr nei = p_nei->back;
243 if (isTip(nei->number, tr->mxtips)) {
244 nodeSet.insert(nei->number);
245 } else {
246 nodeSet.insert(nei->number);
247 nodeSet.insert(nei->next->back->number);
248 nodeSet.insert(nei->next->next->back->number);
249 }
250 p_nei = p_nei->next;
251 }
252 while (q_nei != q) {
253 nodeptr nei = q_nei->back;
254 if (isTip(nei->number, tr->mxtips)) {
255 nodeSet.insert(nei->number);
256 } else {
257 nodeSet.insert(nei->number);
258 nodeSet.insert(nei->next->back->number);
259 nodeSet.insert(nei->next->next->back->number);
260 }
261 q_nei = q_nei->next;
262 }
263 return nodeSet;
264 }
265
pllEvalAllNNIs(pllInstance * tr,partitionList * pr,SearchInfo & searchinfo)266 void pllEvalAllNNIs(pllInstance *tr, partitionList *pr, SearchInfo &searchinfo) {
267 /* DTH: mimic IQTREE::optimizeNNI 's first call to IQTREE::saveCurrentTree */
268 if((globalParams->online_bootstrap == PLL_TRUE) &&
269 (globalParams->gbo_replicates > 0)){
270 tr->fastScaling = PLL_FALSE;
271 pllEvaluateLikelihood(tr, pr, tr->start, PLL_FALSE, PLL_TRUE);
272 pllSaveCurrentTree(tr, pr, tr->start);
273 }
274
275 nodeptr p = tr->start->back;
276 nodeptr q = p->next;
277 while (q != p) {
278 evalNNIForSubtree(tr, pr, q->back, searchinfo);
279 q = q->next;
280 }
281 }
282
283 /*
284 void pllSaveQuartetForSubTree(pllInstance *tr, nodeptr p, SearchInfo &searchinfo) {
285 if (!isTip(p->number, tr->mxtips) && !isTip(p->back->number, tr->mxtips)) {
286 evalNNIForBran(tr, pr, p, searchinfo);
287 nodeptr q = p->next;
288 while (q != p) {
289 pllSaveQuartetForSubTree(tr, q->back, searchinfo);
290 q = q->next;
291 }
292 }
293 }
294
295 void pllSaveAllQuartet(pllInstance *tr, SearchInfo &searchinfo) {
296 nodeptr p = tr->start->back;
297 nodeptr q = p->next;
298 while (q != p) {
299 pllSaveQuartetForSubTree(tr, q->back, searchinfo);
300 }
301 }
302 */
303
pllDoNNISearch(pllInstance * tr,partitionList * pr,SearchInfo & searchinfo)304 double pllDoNNISearch(pllInstance* tr, partitionList *pr, SearchInfo &searchinfo) {
305 double initLH = tr->likelihood;
306 double finalLH = initLH;
307 vector<pllNNIMove> selectedNNIs;
308 unordered_set<int> selectedNodes;
309 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
310
311 /* data structure to store the initial tree topology + branch length */
312 topol* curTree = _setupTopol(tr);
313 saveTree(tr, curTree, numBranches);
314
315 /* evaluate NNIs */
316 pllEvalAllNNIs(tr, pr, searchinfo);
317
318 if (globalParams->speednni) {
319 searchinfo.aBranches.clear();
320 }
321
322 /* apply non-conflicting positive NNIs */
323 if (searchinfo.posNNIList.size() != 0) {
324 sort(searchinfo.posNNIList.begin(), searchinfo.posNNIList.end(), comparePLLNNIMove);
325 if (verbose_mode >= VB_DEBUG) {
326 cout << "curScore: " << searchinfo.curLogl << endl;
327 for (int i = 0; i < searchinfo.posNNIList.size(); i++) {
328 cout << "Logl of positive NNI " << i << " : " << searchinfo.posNNIList[i].likelihood << endl;
329 }
330 }
331 for (vector<pllNNIMove>::reverse_iterator rit = searchinfo.posNNIList.rbegin(); rit != searchinfo.posNNIList.rend(); ++rit) {
332 if (selectedNodes.find((*rit).p->number) == selectedNodes.end() && selectedNodes.find((*rit).p->back->number) == selectedNodes.end()) {
333 selectedNNIs.push_back((*rit));
334 selectedNodes.insert((*rit).p->number);
335 selectedNodes.insert((*rit).p->back->number);
336 } else {
337 continue;
338 }
339 }
340
341 /* Applying all independent NNI moves */
342 searchinfo.curNumAppliedNNIs = selectedNNIs.size();
343 for (vector<pllNNIMove>::iterator it = selectedNNIs.begin(); it != selectedNNIs.end(); it++) {
344 /* do the topological change */
345 doOneNNI(tr, pr, (*it).p, (*it).nniType, TOPO_ONLY);
346 if (globalParams->speednni) {
347 vector<string> aBranches = getAffectedBranches(tr, (*it).p);
348 searchinfo.aBranches.insert(aBranches.begin(), aBranches.end());
349 }
350 updateBranchLengthForNNI(tr, pr, (*it));
351 if (numBranches > 1 && !tr->useRecom) {
352 pllUpdatePartials(tr, pr, (*it).p, PLL_TRUE);
353 pllUpdatePartials(tr, pr, (*it).p->back, PLL_TRUE);
354 } else {
355 pllUpdatePartials(tr, pr, (*it).p, PLL_FALSE);
356 pllUpdatePartials(tr, pr, (*it).p->back, PLL_FALSE);
357 }
358 }
359 if (selectedNNIs.size() != 0) {
360 //pllEvaluateLikelihood(tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
361 pllOptimizeBranchLengths(tr, pr, 1);
362 if (globalParams->count_trees) {
363 countDistinctTrees(tr, pr);
364 }
365 int numNNI = selectedNNIs.size();
366 /* new tree likelihood should not be smaller the likelihood of the computed best NNI */
367 while (tr->likelihood < selectedNNIs.back().likelihood) {
368 if (numNNI == 1) {
369 printf("ERROR: new logl=%10.4f after applying only the best NNI < best NNI logl=%10.4f\n",
370 tr->likelihood, selectedNNIs[0].likelihood);
371 ASSERT(0);
372 } else {
373 cout << "Best logl: " << selectedNNIs.back().likelihood << " / " << "NNI step " << searchinfo.curNumNNISteps<< " / Applying " << numNNI << " NNIs give logl: " << tr->likelihood << " (worse than best)";
374 cout << " / Roll back tree ... " << endl;
375 //restoreTL(rl, tr, 0, pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
376 if (!restoreTree(curTree, tr, pr)) {
377 printf("ERROR: failed to roll back tree \n");
378 ASSERT(0);
379 }
380 // If tree log-likelihood decreases only apply the best NNI
381 numNNI = 1;
382 int count = numNNI;
383 for (vector<pllNNIMove>::reverse_iterator rit = selectedNNIs.rbegin(); rit != selectedNNIs.rend(); ++rit) {
384 doOneNNI(tr, pr, (*rit).p, (*rit).nniType, TOPO_ONLY);
385 updateBranchLengthForNNI(tr, pr, (*rit));
386 if (numBranches > 1 && !tr->useRecom) {
387 pllUpdatePartials(tr, pr, (*rit).p, PLL_TRUE);
388 pllUpdatePartials(tr, pr, (*rit).p->back, PLL_TRUE);
389 } else {
390 pllUpdatePartials(tr, pr, (*rit).p, PLL_FALSE);
391 pllUpdatePartials(tr, pr, (*rit).p->back, PLL_FALSE);
392 }
393 count--;
394 if (count == 0) {
395 break;
396 }
397 }
398 //pllEvaluateLikelihood(tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
399 pllOptimizeBranchLengths(tr, pr, 1);
400 //cout << "Number of NNIs reduced to " << numNNI << ": " << tr->likelihood << endl;
401
402 /* Only apply the best NNI after the tree has been rolled back */
403 searchinfo.curNumAppliedNNIs = numNNI;
404 }
405 }
406 if (tr->likelihood - initLH < 0.1) {
407 searchinfo.curNumAppliedNNIs = 0;
408 }
409 finalLH = tr->likelihood;
410 }
411 } else {
412 searchinfo.curNumAppliedNNIs = 0;
413 }
414 //freeTopol(curTree);
415 return finalLH;
416 }
417
418
updateBranchLengthForNNI(pllInstance * tr,partitionList * pr,pllNNIMove & nni)419 void updateBranchLengthForNNI(pllInstance* tr, partitionList *pr, pllNNIMove &nni) {
420 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
421 /* apply branch lengths */
422 for (int i = 0; i < numBranches; i++) {
423 nni.p->z[i] = nni.z0[i];
424 nni.p->back->z[i] = nni.z0[i];
425 nni.p->next->z[i] = nni.z1[i];
426 nni.p->next->back->z[i] = nni.z1[i];
427 nni.p->next->next->z[i] = nni.z2[i];
428 nni.p->next->next->back->z[i] = nni.z2[i];
429 nni.p->back->next->z[i] = nni.z3[i];
430 nni.p->back->next->back->z[i] = nni.z3[i];
431 nni.p->back->next->next->z[i] = nni.z4[i];
432 nni.p->back->next->next->back->z[i] = nni.z4[i];
433 }
434 /* update partial likelihood */
435 // if (numBranches > 1 && !tr->useRecom) {
436 // pllNewviewGeneric(tr, pr, nni.p, PLL_TRUE);
437 // pllNewviewGeneric(tr, pr, nni.p->back, PLL_TRUE);
438 // } else {
439 // pllNewviewGeneric(tr, pr, nni.p, PLL_FALSE);
440 // pllNewviewGeneric(tr, pr, nni.p->back, PLL_FALSE);
441 // }
442 }
443
pllOptimizeOneBranch(pllInstance * tr,partitionList * pr,nodeptr p)444 void pllOptimizeOneBranch(pllInstance *tr, partitionList *pr, nodeptr p) {
445 nodeptr q;
446 int i;
447 double z[PLL_NUM_BRANCHES], z0[PLL_NUM_BRANCHES];
448 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
449
450 #ifdef _DEBUG_UPDATE
451 double
452 startLH;
453
454 pllEvaluateLikelihood (tr, p);
455
456 startLH = tr->likelihood;
457 #endif
458
459 q = p->back;
460
461 for(i = 0; i < numBranches; i++)
462 z0[i] = q->z[i];
463
464 if(numBranches > 1)
465 makenewzGeneric(tr, pr, p, q, z0, PLL_ITERATIONS, z, PLL_TRUE);
466 else
467 makenewzGeneric(tr, pr, p, q, z0, PLL_ITERATIONS, z, PLL_FALSE);
468
469 for(i = 0; i < numBranches; i++)
470 {
471 if(!tr->partitionConverged[i])
472 {
473 if(PLL_ABS(z[i] - z0[i]) > PLL_DELTAZ)
474 {
475 tr->partitionSmoothed[i] = PLL_FALSE;
476 }
477
478 p->z[i] = q->z[i] = z[i];
479 }
480 }
481
482 #ifdef _DEBUG_UPDATE
483 pllEvaluateLikelihood (tr, p);
484
485 if(tr->likelihood <= startLH)
486 {
487 if(fabs(tr->likelihood - startLH) > 0.01)
488 {
489 printf("%f %f\n", startLH, tr->likelihood);
490 ASSERT(0);
491 }
492 }
493 #endif
494 }
495
doOneNNI(pllInstance * tr,partitionList * pr,nodeptr p,int swap,NNI_Type nni_type,SearchInfo * searchinfo)496 double doOneNNI(pllInstance *tr, partitionList *pr, nodeptr p, int swap, NNI_Type nni_type, SearchInfo *searchinfo) {
497 ASSERT(swap == 0 || swap == 1);
498 nodeptr q;
499 nodeptr tmp;
500 q = p->back;
501 ASSERT(!isTip(q->number, tr->mxtips));
502 ASSERT(!isTip(p->number, tr->mxtips));
503 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
504
505 if (swap == 1) {
506 tmp = p->next->back;
507 hookup(p->next, q->next->back, q->next->z, numBranches);
508 hookup(q->next, tmp, tmp->z, numBranches);
509 } else {
510 tmp = p->next->next->back;
511 hookup(p->next->next, q->next->back, q->next->z, numBranches);
512 hookup(q->next, tmp, tmp->z, numBranches);
513 }
514
515 if (nni_type == TOPO_ONLY) {
516 return 0.0;
517 }
518
519 if (numBranches > 1 && !tr->useRecom) {
520 pllUpdatePartials(tr, pr, p, PLL_TRUE);
521 pllUpdatePartials(tr, pr, q, PLL_TRUE);
522 } else {
523 pllUpdatePartials(tr, pr, p, PLL_FALSE);
524 pllUpdatePartials(tr, pr, q, PLL_FALSE);
525 }
526 // Optimize the central branch
527 pllOptimizeOneBranch(tr, pr, p);
528 if((globalParams->online_bootstrap == PLL_TRUE) && (globalParams->gbo_replicates > 0)){
529 tr->fastScaling = PLL_FALSE;
530 pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_TRUE); // DTH: modified the last arg
531 pllSaveCurrentTree(tr, pr, p);
532 }else{
533 pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_FALSE);
534 }
535 // if after optimizing the central branch we already obtain better logl
536 // then there is no need for optimizing other branches
537 if (tr->likelihood > searchinfo->curLogl) {
538 return tr->likelihood;
539 }
540 // Optimize 4 other branches
541 if (nni_type == NNI5) {
542 if (numBranches > 1 && !tr->useRecom) {
543 pllUpdatePartials(tr, pr, q, PLL_TRUE);
544 } else {
545 pllUpdatePartials(tr, pr, q, PLL_FALSE);
546 }
547 nodeptr r;
548 r = p->next;
549 if (numBranches > 1 && !tr->useRecom) {
550 pllUpdatePartials(tr, pr, r, PLL_TRUE);
551 } else {
552 pllUpdatePartials(tr, pr, r, PLL_FALSE);
553 }
554 pllOptimizeOneBranch(tr, pr, r);
555 // pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_FALSE);
556 // if (tr->likelihood > searchinfo->curLogl) {
557 // return tr->likelihood;
558 // }
559 r = p->next->next;
560 if (numBranches > 1 && !tr->useRecom)
561 pllUpdatePartials(tr, pr, r, PLL_TRUE);
562 else
563 pllUpdatePartials(tr, pr, r, PLL_FALSE);
564 pllOptimizeOneBranch(tr, pr, r);
565 // pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_FALSE);
566 // if (tr->likelihood > searchinfo->curLogl) {
567 // return tr->likelihood;
568 // }
569 if (numBranches > 1 && !tr->useRecom)
570 pllUpdatePartials(tr, pr, p, PLL_TRUE);
571 else
572 pllUpdatePartials(tr, pr, p, PLL_FALSE);
573 pllOptimizeOneBranch(tr, pr, p);
574 // pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_FALSE);
575 // if (tr->likelihood > searchinfo->curLogl) {
576 // return tr->likelihood;
577 // }
578 // optimize 2 branches at node q
579 r = q->next;
580 if (numBranches > 1 && !tr->useRecom)
581 pllUpdatePartials(tr, pr, r, PLL_TRUE);
582 else
583 pllUpdatePartials(tr, pr, r, PLL_FALSE);
584 pllOptimizeOneBranch(tr, pr, r);
585 // pllEvaluateLikelihood(tr, pr, p, PLL_FALSE, PLL_FALSE);
586 // if (tr->likelihood > searchinfo->curLogl) {
587 // return tr->likelihood;
588 // }
589 r = q->next->next;
590 if (numBranches > 1 && !tr->useRecom)
591 pllUpdatePartials(tr, pr, r, PLL_TRUE);
592 else
593 pllUpdatePartials(tr, pr, r, PLL_FALSE);
594 pllOptimizeOneBranch(tr, pr, r);
595 if((globalParams->online_bootstrap == PLL_TRUE) &&
596 (globalParams->gbo_replicates > 0)){
597 tr->fastScaling = PLL_FALSE;
598 pllEvaluateLikelihood(tr, pr, r, PLL_FALSE, PLL_TRUE); // DTH: modified the last arg
599 pllSaveCurrentTree(tr, pr, r);
600 }else{
601 pllEvaluateLikelihood(tr, pr, r, PLL_FALSE, PLL_FALSE);
602 }
603 }
604 return tr->likelihood;
605 }
606
estBestLoglImp(SearchInfo * searchinfo)607 double estBestLoglImp(SearchInfo* searchinfo) {
608 double res = 0.0;
609 int index = floor(searchinfo->deltaLogl.size() * 5 / 100);
610 set<double>::reverse_iterator ri;
611 for (ri = searchinfo->deltaLogl.rbegin(); ri != searchinfo->deltaLogl.rend(); ++ri) {
612 //cout << (*ri) << " ";
613 --index;
614 if (index == 0) {
615 res = (*ri);
616 break;
617 }
618 }
619 //cout << res << endl;
620 //cout << searchinfo->deltaLogl.size() << endl;
621 return res;
622 }
623
convertQuartet2String(nodeptr p)624 string convertQuartet2String(nodeptr p) {
625 nodeptr q = p->back;
626 int pNr = p->number;
627 int qNr = q->number;
628 int pNei1Nr = p->next->back->number;
629 int pNei2Nr = p->next->next->back->number;
630 int qNei1Nr = q->next->back->number;
631 int qNei2Nr = q->next->next->back->number;
632 stringstream middle;
633 stringstream left;
634 stringstream right;
635 stringstream res;
636 if (pNr < qNr) {
637 middle << "-" << pNr << "-" << qNr << "-";
638 } else {
639 middle << "-" << qNr << "-" << pNr << "-";
640 }
641 if (pNei1Nr < pNei2Nr) {
642 left << pNei1Nr << "-" << pNei2Nr;
643 } else {
644 left << pNei2Nr << "-" << pNei1Nr;
645 }
646 if (qNei1Nr < qNei2Nr) {
647 right << qNei1Nr << "-" << qNei2Nr;
648 } else {
649 right << qNei2Nr << "-" << qNei1Nr;
650 }
651 res << left.str() << middle.str() << right.str();
652 return res.str();
653 }
654
countDistinctTrees(pllInstance * pllInst,partitionList * pllPartitions)655 void countDistinctTrees(pllInstance* pllInst, partitionList *pllPartitions) {
656 pllTreeToNewick(pllInst->tree_string, pllInst, pllPartitions, pllInst->start->back, PLL_FALSE,
657 PLL_TRUE, PLL_FALSE, PLL_FALSE, PLL_FALSE, PLL_SUMMARIZE_LH, PLL_FALSE, PLL_FALSE);
658 PhyloTree mtree;
659 mtree.rooted = false;
660 mtree.aln = globalAlignment;
661 mtree.readTreeString(string(pllInst->tree_string));
662 // mtree.root = mtree.findNodeName(globalAlignment->getSeqName(0));
663 mtree.setRootNode(mtree.params->root);
664 ostringstream ostr;
665 mtree.printTree(ostr, WT_TAXON_ID | WT_SORT_TAXA);
666 string tree_str = ostr.str();
667 if (pllTreeCounter.find(tree_str) == pllTreeCounter.end()) {
668 // not found in hash_map
669 pllTreeCounter[tree_str] = 1;
670 } else {
671 // found in hash_map
672 pllTreeCounter[tree_str]++;
673 }
674 }
675
evalNNIForBran(pllInstance * tr,partitionList * pr,nodeptr p,SearchInfo & searchinfo)676 int evalNNIForBran(pllInstance* tr, partitionList *pr, nodeptr p, SearchInfo &searchinfo) {
677 nodeptr q = p->back;
678 ASSERT(!isTip(p->number, tr->mxtips));
679 ASSERT(!isTip(q->number, tr->mxtips));
680 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
681 int numPosNNI = 0;
682 int i;
683 pllNNIMove nni0; // dummy NNI to store backup information
684 nni0.p = p;
685 nni0.nniType = 0;
686 nni0.likelihood = searchinfo.curLogl;
687 for (i = 0; i < PLL_NUM_BRANCHES; i++) {
688 nni0.z0[i] = p->z[i];
689 nni0.z1[i] = p->next->z[i];
690 nni0.z2[i] = p->next->next->z[i];
691 nni0.z3[i] = q->next->z[i];
692 nni0.z4[i] = q->next->next->z[i];
693 }
694
695 pllNNIMove bestNNI;
696
697 /* do an NNI move of type 1 */
698 double lh1 = doOneNNI(tr, pr, p, 0, searchinfo.nni_type, &searchinfo);
699 if (globalParams->count_trees)
700 countDistinctTrees(tr, pr);
701 pllNNIMove nni1;
702 nni1.p = p;
703 nni1.nniType = 0;
704 // Store the optimized branch lengths
705 for (i = 0; i < PLL_NUM_BRANCHES; i++) {
706 nni1.z0[i] = p->z[i];
707 nni1.z1[i] = p->next->z[i];
708 nni1.z2[i] = p->next->next->z[i];
709 nni1.z3[i] = q->next->z[i];
710 nni1.z4[i] = q->next->next->z[i];
711 }
712 nni1.likelihood = lh1;
713 nni1.loglDelta = lh1 - nni0.likelihood;
714 nni1.negLoglDelta = -nni1.loglDelta;
715
716 /* Restore previous NNI move */
717 doOneNNI(tr, pr, p, 0, TOPO_ONLY);
718 /* Restore the old branch length */
719 for (i = 0; i < PLL_NUM_BRANCHES; i++) {
720 p->z[i] = nni0.z0[i];
721 q->z[i] = nni0.z0[i];
722 p->next->z[i] = nni0.z1[i];
723 p->next->back->z[i] = nni0.z1[i];
724 p->next->next->z[i] = nni0.z2[i];
725 p->next->next->back->z[i] = nni0.z2[i];
726 q->next->z[i] = nni0.z3[i];
727 q->next->back->z[i] = nni0.z3[i];
728 q->next->next->z[i] = nni0.z4[i];
729 q->next->next->back->z[i] = nni0.z4[i];
730 }
731
732 /* do an NNI move of type 2 */
733 double lh2 = doOneNNI(tr, pr, p, 1, searchinfo.nni_type, &searchinfo);
734 if (globalParams->count_trees)
735 countDistinctTrees(tr, pr);
736
737 // Create the nniMove struct to store this move
738 pllNNIMove nni2;
739 nni2.p = p;
740 nni2.nniType = 1;
741 // Store the optimized central branch length
742 for (i = 0; i < PLL_NUM_BRANCHES; i++) {
743 nni2.z0[i] = p->z[i];
744 nni2.z1[i] = p->next->z[i];
745 nni2.z2[i] = p->next->next->z[i];
746 nni2.z3[i] = q->next->z[i];
747 nni2.z4[i] = q->next->next->z[i];
748 }
749 nni2.likelihood = lh2;
750 nni2.loglDelta = lh2 - nni0.likelihood;
751 nni2.negLoglDelta = -nni2.loglDelta;
752
753 if (nni2.likelihood > nni1.likelihood) {
754 bestNNI = nni2;
755 } else {
756 bestNNI = nni1;
757 }
758
759 if (bestNNI.likelihood > searchinfo.curLogl + 1e-6) {
760 numPosNNI++;
761 searchinfo.posNNIList.push_back(bestNNI);
762 }
763
764 /* Restore previous NNI move */
765 doOneNNI(tr, pr, p, 1, TOPO_ONLY);
766 /* Restore the old branch length */
767 for (i = 0; i < PLL_NUM_BRANCHES; i++) {
768 p->z[i] = nni0.z0[i];
769 q->z[i] = nni0.z0[i];
770 p->next->z[i] = nni0.z1[i];
771 p->next->back->z[i] = nni0.z1[i];
772 p->next->next->z[i] = nni0.z2[i];
773 p->next->next->back->z[i] = nni0.z2[i];
774 q->next->z[i] = nni0.z3[i];
775 q->next->back->z[i] = nni0.z3[i];
776 q->next->next->z[i] = nni0.z4[i];
777 q->next->next->back->z[i] = nni0.z4[i];
778 }
779
780 // Re-compute the partial likelihood vectors
781 // TODO: One could save these instead of recomputation
782 if (numBranches > 1 && !tr->useRecom) {
783 pllUpdatePartials(tr, pr, p, PLL_TRUE);
784 pllUpdatePartials(tr, pr, p->back, PLL_TRUE);
785 } else {
786 pllUpdatePartials(tr, pr, p, PLL_FALSE);
787 pllUpdatePartials(tr, pr, p->back, PLL_FALSE);
788 }
789
790 return numPosNNI;
791 }
792
793 //void recomputePartial(pllInstance tr, partitionList pr, nodeptr p) {
794 // if (numBranches > 1 && !tr->useRecom) {
795 // pllUpdatePartials(tr, pr, p, PLL_TRUE);
796 // pllUpdatePartials(tr, pr, p->back, PLL_TRUE);
797 // } else {
798 // pllUpdatePartials(tr, pr, p, PLL_FALSE);
799 // pllUpdatePartials(tr, pr, p->back, PLL_FALSE);
800 // }
801 //}
802
isAffectedBranch(nodeptr p,SearchInfo & searchinfo)803 bool isAffectedBranch(nodeptr p, SearchInfo &searchinfo) {
804 string branString = getBranString(p);
805 if (searchinfo.aBranches.find(branString) != searchinfo.aBranches.end()) {
806 return true;
807 } else {
808 return false;
809 }
810 }
811
evalNNIForSubtree(pllInstance * tr,partitionList * pr,nodeptr p,SearchInfo & searchinfo)812 void evalNNIForSubtree(pllInstance* tr, partitionList *pr, nodeptr p, SearchInfo &searchinfo) {
813 if (!isTip(p->number, tr->mxtips) && !isTip(p->back->number, tr->mxtips)) {
814 if (globalParams->speednni && searchinfo.curNumNNISteps != 1) {
815 if (isAffectedBranch(p, searchinfo)) {
816 evalNNIForBran(tr, pr, p, searchinfo);
817 }
818 } else {
819 evalNNIForBran(tr, pr, p, searchinfo);
820 }
821 nodeptr q = p->next;
822 while (q != p) {
823 evalNNIForSubtree(tr, pr, q->back, searchinfo);
824 q = q->next;
825 }
826 }
827 }
828
829 /**
830 * DTH:
831 * The PLL version of saveCurrentTree function
832 * @param tr: the tree (a pointer to a pllInstance)
833 * @param pr: pointer to a partitionList (this one keeps tons of tree info)
834 * @param p: root?
835 */
pllSaveCurrentTree(pllInstance * tr,partitionList * pr,nodeptr p)836 void pllSaveCurrentTree(pllInstance* tr, partitionList *pr, nodeptr p){
837 double cur_logl = tr->likelihood;
838
839 struct pllHashItem * item_ptr = (struct pllHashItem *) malloc(sizeof(struct pllHashItem));
840 item_ptr->data = (int *) malloc(sizeof(int));
841 item_ptr->next = NULL;
842 item_ptr->str = NULL;
843
844 unsigned int tree_index = -1;
845 char * tree_str = NULL;
846 pllTree2StringREC(tr->tree_string, tr, pr, tr->start->back, PLL_FALSE,
847 PLL_FALSE, PLL_FALSE, PLL_FALSE, PLL_TRUE, PLL_SUMMARIZE_LH, PLL_FALSE, PLL_FALSE);
848 tree_str = (char *) malloc (strlen(tr->tree_string) + 1);
849 strcpy(tree_str, tr->tree_string);
850
851 pllBoolean is_stored = PLL_FALSE;
852
853 // if(globalParam->store_candidate_trees){
854 // is_stored = pllHashSearch(pllUFBootDataPtr->treels, tree_str, &(item_ptr->data));
855 // }
856
857 if(is_stored){ // if found the tree_str
858 pllUFBootDataPtr->duplication_counter++;
859 tree_index = *((int *)item_ptr->data);
860 if (cur_logl <= pllUFBootDataPtr->treels_logl[tree_index] + 1e-4) {
861 if (cur_logl < pllUFBootDataPtr->treels_logl[tree_index] - 5.0)
862 if (verbose_mode >= VB_MED)
863 printf("Current lh %f is much worse than expected %f\n",
864 cur_logl, pllUFBootDataPtr->treels_logl[tree_index]);
865 /* free(tree_str);
866 free(item_ptr->data);
867 free(item_ptr);*/
868 return;
869 }
870 if (verbose_mode >= VB_MAX)
871 printf("Updated logl %f to %f\n", pllUFBootDataPtr->treels_logl[tree_index], cur_logl);
872 pllUFBootDataPtr->treels_logl[tree_index] = cur_logl;
873
874 // if (pllUFBootDataPtr->save_all_br_lens) {
875 // pllTree2StringREC(tr->tree_string, tr, pr, tr->start->back, PLL_TRUE,
876 // PLL_FALSE, PLL_FALSE, PLL_FALSE, PLL_TRUE, PLL_SUMMARIZE_LENGTH, PLL_FALSE, PLL_FALSE);
877 // char * tree_str_br_lens = (char *) malloc (strlen(tr->tree_string) + 1);
878 // strcpy(tree_str_br_lens, tr->tree_string);
879 // pllUFBootDataPtr->treels_newick[tree_index] = tree_str_br_lens;
880 // }
881 if (pllUFBootDataPtr->boot_samples == NULL) {
882 (pllUFBootDataPtr->treels_ptnlh)[tree_index] =
883 (double *) malloc(pllUFBootDataPtr->n_patterns * sizeof(double));
884 pllComputePatternLikelihood(tr, (pllUFBootDataPtr->treels_ptnlh)[tree_index], &cur_logl);
885 /* free(tree_str);
886 free(item_ptr->data);
887 free(item_ptr);*/
888 return;
889 }
890 if (verbose_mode >= VB_MAX)
891 printf("Update treels_logl[%d] := %f\n", tree_index, cur_logl);
892
893 } else {
894 if (pllUFBootDataPtr->logl_cutoff != 0.0 && cur_logl <= pllUFBootDataPtr->logl_cutoff + 1e-4){
895 /* free(tree_str);
896 free(item_ptr->data);
897 free(item_ptr);*/
898 return;
899 }
900
901 if(pllUFBootDataPtr->treels_size == pllUFBootDataPtr->candidate_trees_count)
902 pllResizeUFBootData();
903
904 tree_index = pllUFBootDataPtr->candidate_trees_count;
905 pllUFBootDataPtr->candidate_trees_count++;
906 // if (globalParam->store_candidate_trees){
907 // *((int *)item_ptr->data) = tree_index;
908 // item_ptr->str = tree_str;
909 // pllHashAdd(pllUFBootDataPtr->treels, pllHashString(tree_str, pllUFBootDataPtr->treels->size), tree_str, item_ptr->data);
910 // }
911 pllUFBootDataPtr->treels_logl[tree_index] = cur_logl;
912
913 if (verbose_mode >= VB_MAX)
914 printf("Add treels_logl[%d] := %f\n", tree_index, cur_logl);
915 }
916
917 //if (write_intermediate_trees)
918 // printTree(out_treels, WT_NEWLINE | WT_BR_LEN);
919
920 double *pattern_lh = (double *) malloc(pllUFBootDataPtr->n_patterns * sizeof(double));
921 if(!pattern_lh) outError("Not enough dynamic memory!");
922 pllComputePatternLikelihood(tr, pattern_lh, &cur_logl);
923
924 if (pllUFBootDataPtr->boot_samples == NULL) {
925 // for runGuidedBootstrap
926 pllUFBootDataPtr->treels_ptnlh[tree_index] = pattern_lh;
927 } else {
928 // online bootstrap
929 int nptn = pllUFBootDataPtr->n_patterns;
930 int updated = 0;
931 int nsamples = globalParams->gbo_replicates;
932 for (int sample = 0; sample < nsamples; sample++) {
933 double rell = 0.0;
934 for (int ptn = 0; ptn < nptn; ptn++)
935 rell += pattern_lh[ptn] * pllUFBootDataPtr->boot_samples[sample][ptn];
936
937 if (rell > pllUFBootDataPtr->boot_logl[sample] + globalParams->ufboot_epsilon ||
938 (rell > pllUFBootDataPtr->boot_logl[sample] - globalParams->ufboot_epsilon &&
939 random_double() <= 1.0/(pllUFBootDataPtr->boot_counts[sample]+1))) {
940 // if (!globalParam->store_candidate_trees)
941 {
942 is_stored = pllHashSearch(pllUFBootDataPtr->treels, tree_str, &(item_ptr->data));
943 if(is_stored)
944 tree_index = *((int *)item_ptr->data);
945 else{
946 *((int *)item_ptr->data) = tree_index = pllUFBootDataPtr->candidate_trees_count - 1;
947 item_ptr->str = tree_str;
948 pllHashAdd(pllUFBootDataPtr->treels, pllHashString(tree_str, pllUFBootDataPtr->treels->size), tree_str, item_ptr->data);
949 // pllHashAdd(pllUFBootDataPtr->treels, tree_str, item_ptr->data);
950 }
951 }
952 if (rell <= pllUFBootDataPtr->boot_logl[sample] +
953 globalParams->ufboot_epsilon) {
954 pllUFBootDataPtr->boot_counts[sample]++;
955 } else {
956 pllUFBootDataPtr->boot_counts[sample] = 1;
957 }
958 if(rell > pllUFBootDataPtr->boot_logl[sample])
959 pllUFBootDataPtr->boot_logl[sample] = rell;
960 pllUFBootDataPtr->boot_trees[sample] = tree_index;
961 updated++;
962 }
963 }
964 /* if (updated && verbose_mode >= VB_MAX)
965 printf("%d boot trees updated\n", updated);*/
966 }
967 // if (pllUFBootDataPtr->save_all_br_lens) {
968 // pllTree2StringREC(tr->tree_string, tr, pr, tr->start->back, PLL_TRUE,
969 // PLL_FALSE, PLL_FALSE, PLL_FALSE, PLL_TRUE, PLL_SUMMARIZE_LH, PLL_FALSE, PLL_FALSE);
970 // char * s = (char *) malloc (strlen(tr->tree_string) + 1);
971 // strcpy(s, tr->tree_string);
972 // pllUFBootDataPtr->treels_newick[tree_index] = s;
973 // }
974
975 // if(!globalParam->store_candidate_trees){
976 // free(tree_str);
977 // free(item_ptr->data);
978 // free(item_ptr);
979 // }
980 if (pllUFBootDataPtr->boot_samples){
981 free(pattern_lh);
982 pllUFBootDataPtr->treels_ptnlh[tree_index] = NULL;
983 }
984
985 // printf("Done freeing: max = %d, count = %d, size = %d\n",
986 // pllUFBootDataPtr->max_candidate_trees,
987 // pllUFBootDataPtr->candidate_trees_count,
988 // pllUFBootDataPtr->treels_size);
989 }
990
991 /**
992 * DTH:
993 * Extract the array of site log likelihood to be kept in ptnlh
994 * And update *cur_log
995 * @param tr: the tree (pointer to an pllInstance)
996 * @param ptnlh: to-be-kept array of site log likelihood
997 * @param cur_logl: pointer to current tree log likelihood
998 */
pllComputePatternLikelihood(pllInstance * tr,double * ptnlh,double * cur_logl)999 void pllComputePatternLikelihood(pllInstance* tr, double * ptnlh, double * cur_logl){
1000 int i;
1001 double tree_logl = 0;
1002 for(i = 0; i < pllUFBootDataPtr->n_patterns; i++){
1003 ptnlh[i] = tr->lhs[i];
1004 tree_logl += tr->lhs[i] * tr->aliaswgt[i];
1005 }
1006 *cur_logl = tree_logl;
1007 }
1008
1009 /**
1010 * DTH:
1011 * Resize some of the arrays in UFBootData if they're full
1012 * Along with update treels_size (to track the size of these arrays)
1013 */
pllResizeUFBootData()1014 void pllResizeUFBootData(){
1015 int count = pllUFBootDataPtr->candidate_trees_count;
1016 pllUFBootDataPtr->treels_size = 2 * count;
1017
1018 double * rtreels_logl =
1019 (double *) malloc(2 * count * (sizeof(double)));
1020 if(!rtreels_logl) outError("Not enough dynamic memory!");
1021 // memset(rtreels_logl, 0, 2 * count * sizeof(double));
1022 memcpy(rtreels_logl, pllUFBootDataPtr->treels_logl, count * sizeof(double));
1023 free(pllUFBootDataPtr->treels_logl);
1024 pllUFBootDataPtr->treels_logl = rtreels_logl;
1025
1026 // char ** rtreels_newick =
1027 // (char **) malloc(2 * count * (sizeof(char *)));
1028 // if(!rtreels_newick) outError("Not enough dynamic memory!");
1029 // memset(rtreels_newick, 0, 2 * count * sizeof(char *));
1030 // memcpy(rtreels_newick, pllUFBootDataPtr->treels_newick, count * sizeof(char *));
1031 // free(pllUFBootDataPtr->treels_newick);
1032 // pllUFBootDataPtr->treels_newick = rtreels_newick;
1033
1034 double ** rtreels_ptnlh =
1035 (double **) malloc(2 * count * (sizeof(double *)));
1036 if(!rtreels_ptnlh) outError("Not enough dynamic memory!");
1037 memset(rtreels_ptnlh, 0, 2 * count * sizeof(double *));
1038 memcpy(rtreels_ptnlh, pllUFBootDataPtr->treels_ptnlh, count * sizeof(double *));
1039 free(pllUFBootDataPtr->treels_ptnlh);
1040 pllUFBootDataPtr->treels_ptnlh = rtreels_ptnlh;
1041 }
1042
1043
1044 /**
1045 * DTH:
1046 * (Based on function Tree2StringREC of PLL)
1047 * Print out the tree topology with IQTree taxa ID (starts at 0) instead of PLL taxa ID (starts at 1)
1048 * @param All are the same as in PLL's
1049 * 2014.4.23: REPLACE getBranchLength(tr, pr, perGene, p) BY pllGetBranchLength(tr, p, pr->numberOfPartitions)
1050 * BECAUSE OF LIB NEW DECLARATION: pllGetBranchLength (pllInstance *tr, nodeptr p, int partition_id);
1051 */
pllTree2StringREC(char * treestr,pllInstance * tr,partitionList * pr,nodeptr p,pllBoolean printBranchLengths,pllBoolean printNames,pllBoolean printLikelihood,pllBoolean rellTree,pllBoolean finalPrint,int perGene,pllBoolean branchLabelSupport,pllBoolean printSHSupport)1052 char *pllTree2StringREC(char *treestr, pllInstance *tr, partitionList *pr, nodeptr p, pllBoolean printBranchLengths, pllBoolean printNames,
1053 pllBoolean printLikelihood, pllBoolean rellTree, pllBoolean finalPrint, int perGene, pllBoolean branchLabelSupport, pllBoolean printSHSupport)
1054 {
1055 char * result = treestr; // DTH: added this var to be able to remove the '\n' at the end
1056 char *nameptr;
1057
1058 if(isTip(p->number, tr->mxtips)){
1059 if(printNames){
1060 nameptr = tr->nameList[p->number];
1061 sprintf(treestr, "%s", nameptr);
1062 }else
1063 sprintf(treestr, "%d", p->number - 1);
1064
1065 while (*treestr) treestr++;
1066 }else{
1067 *treestr++ = '(';
1068 treestr = pllTree2StringREC(treestr, tr, pr, p->next->back, printBranchLengths, printNames, printLikelihood, rellTree,
1069 finalPrint, perGene, branchLabelSupport, printSHSupport);
1070 *treestr++ = ',';
1071 treestr = pllTree2StringREC(treestr, tr, pr, p->next->next->back, printBranchLengths, printNames, printLikelihood, rellTree,
1072 finalPrint, perGene, branchLabelSupport, printSHSupport);
1073 if(p == tr->start->back){
1074 *treestr++ = ',';
1075 treestr = pllTree2StringREC(treestr, tr, pr, p->back, printBranchLengths, printNames, printLikelihood, rellTree,
1076 finalPrint, perGene, branchLabelSupport, printSHSupport);
1077 }
1078 *treestr++ = ')';
1079 }
1080
1081 if(p == tr->start->back){
1082 if(printBranchLengths && !rellTree)
1083 sprintf(treestr, ":0.0;\n");
1084 else
1085 sprintf(treestr, ";\n");
1086 }else{
1087 if(rellTree || branchLabelSupport || printSHSupport){
1088 if(( !isTip(p->number, tr->mxtips)) &&
1089 ( !isTip(p->back->number, tr->mxtips)))
1090 {
1091 ASSERT(p->bInf != (branchInfo *)NULL);
1092 if(rellTree)
1093 sprintf(treestr, "%d:%8.20f", p->bInf->support, p->z[0]);
1094 if(branchLabelSupport)
1095 sprintf(treestr, ":%8.20f[%d]", p->z[0], p->bInf->support);
1096 if(printSHSupport)
1097 sprintf(treestr, ":%8.20f[%d]", pllGetBranchLength(tr, p, pr->numberOfPartitions), p->bInf->support);
1098 }else{
1099 if(rellTree || branchLabelSupport)
1100 sprintf(treestr, ":%8.20f", p->z[0]);
1101 if(printSHSupport)
1102 sprintf(treestr, ":%8.20f", pllGetBranchLength(tr, p, pr->numberOfPartitions));
1103 }
1104 }else{
1105 if(printBranchLengths)
1106 sprintf(treestr, ":%8.20f", pllGetBranchLength(tr, p, pr->numberOfPartitions));
1107 else
1108 sprintf(treestr, "%s", "\0");
1109 }
1110 }
1111
1112 if(result[strlen(result) - 1] == '\n') result[strlen(result) - 1] = '\0';
1113 while (*treestr) treestr++;
1114 return treestr;
1115 }
1116
1117
1118