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