1 /**
2 * PLL (version 1.0.0) a software library for phylogenetic inference
3 * Copyright (C) 2013 Tomas Flouri and Alexandros Stamatakis
4 *
5 * Derived from
6 * RAxML-HPC, a program for sequential and parallel estimation of phylogenetic
7 * trees by Alexandros Stamatakis
8 *
9 * This program is free software: you can redistribute it and/or modify it
10 * under the terms of the GNU General Public License as published by the Free
11 * Software Foundation, either version 3 of the License, or (at your option)
12 * any later version.
13 *
14 * This program is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
17 * more details.
18 *
19 * You should have received a copy of the GNU General Public License along with
20 * this program. If not, see <http://www.gnu.org/licenses/>.
21 *
22 * For any other enquiries send an Email to Tomas Flouri
23 * Tomas.Flouri@h-its.org
24 *
25 * When publishing work that uses PLL please cite PLL
26 *
27 * @file searchAlgo.c
28 * @brief Collection of routines for performing likelihood computation and branch optimization.
29 *
30 * Detailed description to appear soon.
31 */
32 #include "mem_alloc.h"
33
34 #ifndef WIN32
35 #include <sys/times.h>
36 #include <sys/types.h>
37 #include <sys/time.h>
38 #include <unistd.h>
39 #endif
40
41 #include <math.h>
42 #include <time.h>
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <ctype.h>
46 #include <string.h>
47 #include <assert.h>
48 #include <errno.h>
49
50 #include "pll.h"
51 #include "pllInternal.h"
52
53 typedef struct bInf {
54 double likelihood;
55 nodeptr node;
56 } bestInfo;
57
58 typedef struct iL {
59 bestInfo *list;
60 int n;
61 int valid;
62 } infoList;
63
64 double treeOptimizeRapid(pllInstance *tr, partitionList *pr, int mintrav, int maxtrav, bestlist *bt, infoList *iList);
65 nniMove getBestNNIForBran(pllInstance* tr, partitionList *pr, nodeptr p, double curLH);
66 void evalNNIForSubtree(pllInstance* tr, partitionList *pr, nodeptr p, nniMove* nniList, int* cnt, int* cnt_nni, double curLH);
67
68
69 static int cmp_nni(const void* nni1, const void* nni2);
70 static void pllTraverseUpdate (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q, int mintrav, int maxtrav, pllRearrangeList * bestList);
71 static int pllStoreRearrangement (pllRearrangeList * bestList, pllRearrangeInfo * rearr);
72 static int pllTestInsertBIG (pllInstance * tr, partitionList * pr, nodeptr p, nodeptr q, pllRearrangeList * bestList);
73 static int pllTestSPR (pllInstance * tr, partitionList * pr, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList);
74 static void pllCreateSprInfoRollback (pllInstance * tr, pllRearrangeInfo * rearr, int numBranches);
75 static void pllCreateNniInfoRollback (pllInstance * tr, pllRearrangeInfo * rearr);
76 static void pllCreateRollbackInfo (pllInstance * tr, pllRearrangeInfo * rearr, int numBranches);
77 static void pllRollbackNNI (pllInstance * tr, partitionList * pr, pllRollbackInfo * ri);
78 static void pllRollbackSPR (partitionList * pr, pllRollbackInfo * ri);
79
80 extern partitionLengths pLengths[PLL_MAX_MODEL];
81
initrav(pllInstance * tr,partitionList * pr,nodeptr p)82 pllBoolean initrav (pllInstance *tr, partitionList *pr, nodeptr p)
83 {
84 nodeptr q;
85
86 if (!isTip(p->number, tr->mxtips))
87 {
88 q = p->next;
89
90 do
91 {
92 if (! initrav(tr, pr, q->back)) return PLL_FALSE;
93 q = q->next;
94 }
95 while (q != p);
96
97 pllUpdatePartials(tr, pr, p, PLL_FALSE);
98 }
99
100 return PLL_TRUE;
101 }
102
103
104 /** @brief Optimize the length of a specific branch
105
106 Optimize the length of the branch connecting \a p and \a p->back
107 for each partition (\a tr->numBranches) in library instance \a tr.
108
109 @param tr
110 The library instance
111
112 @param pr
113 Partition list
114
115 @param p
116 Endpoints of branch to be optimized
117 */
update(pllInstance * tr,partitionList * pr,nodeptr p)118 void update(pllInstance *tr, partitionList *pr, nodeptr p)
119 {
120 nodeptr q;
121 int i;
122 double z[PLL_NUM_BRANCHES], z0[PLL_NUM_BRANCHES];
123 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
124
125 #ifdef _DEBUG_UPDATE
126 double
127 startLH;
128
129 pllEvaluateLikelihood (tr, p);
130
131 startLH = tr->likelihood;
132 #endif
133
134 q = p->back;
135
136 for(i = 0; i < numBranches; i++)
137 z0[i] = q->z[i];
138
139 if(numBranches > 1)
140 makenewzGeneric(tr, pr, p, q, z0, PLL_NEWZPERCYCLE, z, PLL_TRUE);
141 else
142 makenewzGeneric(tr, pr, p, q, z0, PLL_NEWZPERCYCLE, z, PLL_FALSE);
143
144 for(i = 0; i < numBranches; i++)
145 {
146 if(!tr->partitionConverged[i])
147 {
148 if(PLL_ABS(z[i] - z0[i]) > PLL_DELTAZ)
149 {
150 tr->partitionSmoothed[i] = PLL_FALSE;
151 }
152
153 p->z[i] = q->z[i] = z[i];
154 }
155 }
156
157 #ifdef _DEBUG_UPDATE
158 pllEvaluateLikelihood (tr, p);
159
160 if(tr->likelihood <= startLH)
161 {
162 if(fabs(tr->likelihood - startLH) > 0.01)
163 {
164 printf("%f %f\n", startLH, tr->likelihood);
165 assert(0);
166 }
167 }
168 #endif
169 }
170
171 /** @brief Branch length optimization of subtree
172
173 Optimize the length of branch connected by \a p and \a p->back, and the
174 lengths of all branches in the subtrees rooted at \a p->next and \a p->next->next
175
176 @param tr
177 The library instance
178
179 @param pr
180 Partition list
181
182 @param p
183 Endpoint of branches to be optimized
184 */
smooth(pllInstance * tr,partitionList * pr,nodeptr p)185 void smooth (pllInstance *tr, partitionList *pr, nodeptr p)
186 {
187 nodeptr q;
188 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
189
190 update(tr, pr, p); /* Adjust branch */
191
192 if (! isTip(p->number, tr->mxtips))
193 { /* Adjust descendants */
194 q = p->next;
195 while (q != p)
196 {
197 smooth(tr, pr, q->back);
198 q = q->next;
199 }
200
201 if(numBranches > 1 && !tr->useRecom)
202 pllUpdatePartials(tr, pr,p, PLL_TRUE);
203 else
204 pllUpdatePartials(tr, pr,p, PLL_FALSE);
205 }
206 }
207
208 /** @brief Check whether the branches in all partitions have been optimized
209
210 Check if all branches in all partitions have reached the threshold for
211 optimization. If at least one branch can be optimized further return \b PLL_FALSE.
212
213 @param tr
214 The library instance
215
216 @return
217 If at least one branch can be further optimized return \b PLL_FALSE,
218 otherwise \b PLL_TRUE.
219
220 */
allSmoothed(pllInstance * tr,int numBranches)221 static pllBoolean allSmoothed(pllInstance *tr, int numBranches)
222 {
223 int i;
224 pllBoolean result = PLL_TRUE;
225
226 for(i = 0; i < numBranches; i++)
227 {
228 if(tr->partitionSmoothed[i] == PLL_FALSE)
229 result = PLL_FALSE;
230 else
231 tr->partitionConverged[i] = PLL_TRUE;
232 }
233
234 return result;
235 }
236
237
238 /** @brief Optimize all branch lenghts of a tree
239
240 Perform \a maxtimes rounds of branch length optimization by running smooth()
241 on all neighbour nodes of node \a tr->start.
242
243 @param tr
244 The library instance
245
246 @param maxtimes
247 Number of optimization rounds to perform
248 */
249 /* do maxtimes rounds of branch length optimization */
smoothTree(pllInstance * tr,partitionList * pr,int maxtimes)250 void smoothTree (pllInstance *tr, partitionList *pr, int maxtimes)
251 {
252 nodeptr p, q;
253 int i, count = 0;
254 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
255
256 p = tr->start;
257 for(i = 0; i < numBranches; i++)
258 tr->partitionConverged[i] = PLL_FALSE;
259
260 while (--maxtimes >= 0)
261 {
262 for(i = 0; i < numBranches; i++)
263 tr->partitionSmoothed[i] = PLL_TRUE;
264
265 smooth(tr, pr, p->back);
266 if (!isTip(p->number, tr->mxtips))
267 {
268 q = p->next;
269 while (q != p)
270 {
271 smooth(tr, pr, q->back);
272 q = q->next;
273 }
274 }
275 count++;
276
277 if (allSmoothed(tr, numBranches)) break;
278 }
279
280 for(i = 0; i < numBranches; i++)
281 tr->partitionConverged[i] = PLL_FALSE;
282 }
283
284
285 /** @brief Optimize the branch length of edges around a specific node
286
287 Optimize \a maxtimes the branch length of all (3) edges around a given node
288 \a p of the tree of library instance \a tr.
289
290 @param tr
291 The library instance
292
293 @param p
294 The node around which to optimize the edges
295
296 @param maxtimes
297 Number of optimization rounds to perform
298 */
localSmooth(pllInstance * tr,partitionList * pr,nodeptr p,int maxtimes)299 void localSmooth (pllInstance *tr, partitionList *pr, nodeptr p, int maxtimes)
300 {
301 nodeptr q;
302 int i;
303 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
304 if (isTip(p->number, tr->mxtips)) return;
305
306 for(i = 0; i < PLL_NUM_BRANCHES; i++)
307 tr->partitionConverged[i] = PLL_FALSE;
308
309 while (--maxtimes >= 0)
310 {
311 for(i = 0; i < PLL_NUM_BRANCHES; i++)
312 tr->partitionSmoothed[i] = PLL_TRUE;
313
314 q = p;
315 do
316 {
317 update(tr, pr, q);
318 q = q->next;
319 }
320 while (q != p);
321
322 if (allSmoothed(tr, numBranches))
323 break;
324 }
325
326 for(i = 0; i < PLL_NUM_BRANCHES; i++)
327 {
328 tr->partitionSmoothed[i] = PLL_FALSE;
329 tr->partitionConverged[i] = PLL_FALSE;
330 }
331 }
332
333
334
335
336 /** @brief Reset an \a infoList
337
338 Resets an \a infoList by setting elements \a node and \a likelihood
339 of each element of the \a bestInfo list structure to \b NULL and
340 \a PLL_UNLIKELY, respectively.
341
342 @param iList
343 The given \a infoList.
344 */
resetInfoList(infoList * iList)345 static void resetInfoList(infoList *iList)
346 {
347 int
348 i;
349
350 iList->valid = 0;
351
352 for(i = 0; i < iList->n; i++)
353 {
354 iList->list[i].node = (nodeptr)NULL;
355 iList->list[i].likelihood = PLL_UNLIKELY;
356 }
357 }
358
359 /** @brief Initialize an \a infoList
360
361 Initialize an \a infoList by creating a \a bestInfo list structure
362 of \a n elements and setting the attributes \a node and \a likelihood
363 of each element of the \a bestInfo list structure to \b NULL and
364 \a PLL_UNLIKELY, respectively.
365
366 @param iList
367 The given \a infoList.
368
369 @param n
370 Number of elements to be created in the \a bestInfo list.
371 */
initInfoList(infoList * iList,int n)372 static void initInfoList(infoList *iList, int n)
373 {
374 int
375 i;
376
377 iList->n = n;
378 iList->valid = 0;
379 iList->list = (bestInfo *)rax_malloc(sizeof(bestInfo) * (size_t)n);
380
381 for(i = 0; i < n; i++)
382 {
383 iList->list[i].node = (nodeptr)NULL;
384 iList->list[i].likelihood = PLL_UNLIKELY;
385 }
386 }
387
388 /** @brief Deallocate the contents of an \a infoList
389
390 Deallocate the contents of a given \a infoList by freeing
391 the memory used by its \a bestInfo list structure.
392
393 @param iList
394 The \a infoList to be used.
395 */
freeInfoList(infoList * iList)396 static void freeInfoList(infoList *iList)
397 {
398 rax_free(iList->list);
399 }
400
401
402 /** @brief Insert a record in an \a infoList
403
404 Insert the pair \a likelihood and \node into list \a iList
405 \b only if there already exists a pair in \a iList
406 whose \a likelihood attribute is smaller than the given \a
407 likelihood. The insertion is done by replacing the smallest
408 likelihood pair with the new pair.
409
410 @param node
411 The given node
412
413 @param likelihood
414 The given likelihood
415
416 @param iList
417 The given \a infoList where the record will possibly be appended.
418 */
insertInfoList(nodeptr node,double likelihood,infoList * iList)419 static void insertInfoList(nodeptr node, double likelihood, infoList *iList)
420 {
421 int
422 i,
423 min = 0;
424
425 double
426 min_l = iList->list[0].likelihood;
427
428 for(i = 1; i < iList->n; i++)
429 {
430 if(iList->list[i].likelihood < min_l)
431 {
432 min = i;
433 min_l = iList->list[i].likelihood;
434 }
435 }
436
437 if(likelihood > min_l)
438 {
439 iList->list[min].likelihood = likelihood;
440 iList->list[min].node = node;
441 if(iList->valid < iList->n)
442 iList->valid += 1;
443 }
444 }
445
446
447 /** @brief Optimize branch lengths of region
448
449 Optimize the branch lenghts of only a specific region. The branch optimization starts
450 at a node \a p and is carried out in all nodes with distance upto \a region edges from
451 \a p.
452
453 @param tr
454 The library instance.
455
456 @param p
457 Node to start branch optimization from.
458
459 @param region
460 The allowed node distance from \p for which to still perform branch optimization.
461 */
smoothRegion(pllInstance * tr,partitionList * pr,nodeptr p,int region)462 void smoothRegion (pllInstance *tr, partitionList *pr, nodeptr p, int region)
463 {
464 nodeptr q;
465
466 update(tr, pr, p); /* Adjust branch */
467
468 if (region > 0)
469 {
470 if (!isTip(p->number, tr->mxtips))
471 {
472 q = p->next;
473 while (q != p)
474 {
475 smoothRegion(tr, pr, q->back, --region);
476 q = q->next;
477 }
478
479 pllUpdatePartials(tr, pr,p, PLL_FALSE);
480 }
481 }
482 }
483
484 /** @brief Wrapper function for optimizing the branch length of a region \a maxtimes times
485
486 Optimize the branch lengths of a specific region \a maxtimes times. The branch optimization
487 starts at a given node \a p and is carried out in all nodes with distance upto \a region
488 from \a p.
489
490 @param tr
491 The library instance.
492
493 @param p
494 Node to start branch optimization from.
495
496 @param maxtimes
497 Number of times to perform branch optimization.
498
499 @param region
500 The allwed node distance from \p for which to still perform branch optimization.
501
502 @todo
503 In the previous version (before the model-sep merge) the loops were controlled by tr->numBranches,
504 and now they are controlled by a constant PLL_NUM_BRANCHES. What is right?
505 */
regionalSmooth(pllInstance * tr,partitionList * pr,nodeptr p,int maxtimes,int region)506 void regionalSmooth (pllInstance *tr, partitionList *pr, nodeptr p, int maxtimes, int region)
507 {
508 nodeptr q;
509 int i;
510 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
511
512 if (isTip(p->number, tr->mxtips)) return; /* Should be an error */
513
514 for(i = 0; i < PLL_NUM_BRANCHES; i++)
515 tr->partitionConverged[i] = PLL_FALSE;
516
517 while (--maxtimes >= 0)
518 {
519 for(i = 0; i < PLL_NUM_BRANCHES; i++)
520 tr->partitionSmoothed[i] = PLL_TRUE;
521
522 q = p;
523 do
524 {
525 smoothRegion(tr, pr, q, region);
526 q = q->next;
527 }
528 while (q != p);
529
530 if (allSmoothed(tr, numBranches))
531 break;
532 }
533
534 for(i = 0; i < PLL_NUM_BRANCHES; i++) {
535 tr->partitionSmoothed[i] = PLL_FALSE;
536 tr->partitionConverged[i] = PLL_FALSE;
537 }
538 }
539
540
541
542
543 /** @brief Split the tree into two components and optimize new branch length
544
545 Split the tree into two components. The disconnection point is node \a p.
546 First, a branch length is computed for the newly created branch between nodes
547 \a p->next->back and \a p->next->next->back and then the two nodes are
548 connected (hookup). Disconnection is done by setting \a p->next->next->back
549 and \a p->next->back to \b NULL.
550
551 @param tr
552 The library instance
553
554 @param p
555 The node at which the tree should be decomposed into two components.
556
557 @param numBranches
558 Number of branches per partition
559
560 @return
561 Node from the disconnected component
562
563 @todo
564 Why do we return this node?
565
566 @image html removeBIG.png "The diagram shows in blue color the new edge that is created and in red the edges that are removed"
567 */
removeNodeBIG(pllInstance * tr,partitionList * pr,nodeptr p,int numBranches)568 nodeptr removeNodeBIG (pllInstance *tr, partitionList *pr, nodeptr p, int numBranches)
569 {
570 // double zqr[numBranches], result[numBranches];
571 double* zqr = rax_malloc(numBranches*sizeof(double)), *result = rax_malloc(numBranches*sizeof(double));
572 nodeptr q, r;
573 int i;
574
575 q = p->next->back;
576 r = p->next->next->back;
577
578 for(i = 0; i < numBranches; i++)
579 zqr[i] = q->z[i] * r->z[i];
580
581 makenewzGeneric(tr, pr, q, r, zqr, PLL_ITERATIONS, result, PLL_FALSE);
582
583 for(i = 0; i < numBranches; i++)
584 tr->zqr[i] = result[i];
585
586 hookup(q, r, result, numBranches);
587
588 p->next->next->back = p->next->back = (node *) NULL;
589
590 rax_free(result);
591 rax_free(zqr);
592 return q;
593 }
594
595 /** @brief Split the tree into two components and recompute likelihood
596
597 Split the tree into two component. The disconnection point is node \a p.
598 Set the branch length of the new node between \a p->next->back and
599 \a p->next->next->back to \a tr->currentZQR and then decompose the tree
600 into two components by setting \a p->next->back and \a p->next->next->back
601 to \b NULL.
602
603 @param tr
604 The library instance
605
606 @param p
607 The node at which the tree should be decomposed into two components.
608
609 @return q
610 the node after \a p
611
612 @todo
613 Why do we return this node? Why do we set to tr->currentZQR and not compute
614 new optimized length? What is tr->currentZQR?
615 */
removeNodeRestoreBIG(pllInstance * tr,partitionList * pr,nodeptr p)616 nodeptr removeNodeRestoreBIG (pllInstance *tr, partitionList *pr, nodeptr p)
617 {
618 nodeptr q, r;
619
620 q = p->next->back;
621 r = p->next->next->back;
622
623 pllUpdatePartials(tr, pr,q, PLL_FALSE);
624 pllUpdatePartials(tr, pr,r, PLL_FALSE);
625
626 hookup(q, r, tr->currentZQR, pr->perGeneBranchLengths?pr->numberOfPartitions:1);
627
628 p->next->next->back = p->next->back = (node *) NULL;
629
630 return q;
631 }
632
633 /** @brief Connect two disconnected tree components
634
635 Connect two disconnected components by specifying an internal edge from one
636 component and a leaf from the other component. The internal edge \a e is the
637 edge between \a q and \a q->back. The leaf is specified by \a p.
638 Edge \a e is removed and two new edges are created. The first one is an edge
639 between \a p->next and \a q, and the second one is between \a p->next->next
640 and \a q->back. The new likelihood vector for node \a p is computed.
641
642 @note The function makes use of the \a thoroughInsertion flag
643
644 @todo
645 What is tr->lzi ? What is thorough insertion? Why do we optimize branch lengths
646 that will be removed? Add explanation
647
648 @image html pll.png "The diagram shows in blue colors the new edges that are created and in red the edge that is removed"
649 */
insertBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q)650 pllBoolean insertBIG (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q)
651 {
652 nodeptr r, s;
653 int i;
654 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
655
656 r = q->back;
657 s = p->back;
658
659 for(i = 0; i < numBranches; i++)
660 tr->lzi[i] = q->z[i];
661
662 if(tr->thoroughInsertion)
663 {
664 double * zqr = rax_malloc(numBranches*sizeof(double)),
665 *zqs = rax_malloc(numBranches*sizeof(double)),
666 *zrs = rax_malloc(numBranches*sizeof(double));
667 double lzqr, lzqs, lzrs, lzsum, lzq, lzr, lzs, lzmax;
668 double *defaultArray=rax_malloc(numBranches*sizeof(double));
669 double *e1 = rax_malloc(numBranches*sizeof(double)),
670 *e2 = rax_malloc(numBranches*sizeof(double)),
671 *e3 = rax_malloc(numBranches*sizeof(double));
672 double *qz;
673
674 qz = q->z;
675
676 for(i = 0; i < numBranches; i++)
677 defaultArray[i] = PLL_DEFAULTZ;
678
679 makenewzGeneric(tr, pr, q, r, qz, PLL_ITERATIONS, zqr, PLL_FALSE);
680 /* the branch lengths values will be estimated using q, r and s
681 * q-s are not connected, but both q and s have a valid LH vector , so we can call makenewzGeneric to get a value for
682 * lzsum, which is then use to generate reasonable starting values e1, e2, e3 for the new branches we create after the insertion
683 */
684
685 makenewzGeneric(tr, pr, q, s, defaultArray, PLL_ITERATIONS, zqs, PLL_FALSE);
686 makenewzGeneric(tr, pr, r, s, defaultArray, PLL_ITERATIONS, zrs, PLL_FALSE);
687
688
689 for(i = 0; i < numBranches; i++)
690 {
691 lzqr = (zqr[i] > PLL_ZMIN) ? log(zqr[i]) : log(PLL_ZMIN);
692 lzqs = (zqs[i] > PLL_ZMIN) ? log(zqs[i]) : log(PLL_ZMIN);
693 lzrs = (zrs[i] > PLL_ZMIN) ? log(zrs[i]) : log(PLL_ZMIN);
694 lzsum = 0.5 * (lzqr + lzqs + lzrs);
695
696 lzq = lzsum - lzrs;
697 lzr = lzsum - lzqs;
698 lzs = lzsum - lzqr;
699 lzmax = log(PLL_ZMAX);
700
701 if (lzq > lzmax) {lzq = lzmax; lzr = lzqr; lzs = lzqs;}
702 else if (lzr > lzmax) {lzr = lzmax; lzq = lzqr; lzs = lzrs;}
703 else if (lzs > lzmax) {lzs = lzmax; lzq = lzqs; lzr = lzrs;}
704
705 e1[i] = exp(lzq);
706 e2[i] = exp(lzr);
707 e3[i] = exp(lzs);
708 }
709 hookup(p->next, q, e1, numBranches);
710 hookup(p->next->next, r, e2, numBranches);
711 hookup(p, s, e3, numBranches);
712 rax_free(e3);
713 rax_free(e2);
714 rax_free(e1);
715 rax_free(defaultArray);
716 rax_free(zrs);
717 rax_free(zqs);
718 rax_free(zqr);
719
720 }
721 else
722 {
723 double *z = rax_malloc(numBranches*sizeof(double));
724
725 for(i = 0; i < numBranches; i++)
726 {
727 z[i] = sqrt(q->z[i]);
728
729 if(z[i] < PLL_ZMIN)
730 z[i] = PLL_ZMIN;
731 if(z[i] > PLL_ZMAX)
732 z[i] = PLL_ZMAX;
733 }
734
735 hookup(p->next, q, z, numBranches);
736 hookup(p->next->next, r, z, numBranches);
737 rax_free(z);
738 }
739
740 pllUpdatePartials(tr, pr,p, PLL_FALSE);
741
742 if(tr->thoroughInsertion)
743 {
744 localSmooth(tr, pr, p, PLL_MAX_LOCAL_SMOOTHING_ITERATIONS);
745 for(i = 0; i < numBranches; i++)
746 {
747 tr->lzq[i] = p->next->z[i];
748 tr->lzr[i] = p->next->next->z[i];
749 tr->lzs[i] = p->z[i];
750 }
751 }
752
753 return PLL_TRUE;
754 }
755
756 /** @brief Connect two disconnected tree components without optimizing branch lengths
757
758 Connect two disconnected components by specifying an internal edge from one
759 component and a leaf from the other component. The internal edge \a e is the
760 edge between \a q and \a q->back. The leaf is specified by \a p.
761 Edge \a e is removed and two new edges are created. The first one is an edge
762 between \a p->next and \a q, and the second one is between \a p->next->next
763 and \a q->back. The new likelihood vector for node \a p is computed.
764
765 @note The function makes use of the \a thoroughInsertion flag
766
767 @todo
768 What is the difference between this and insertBIG?
769 */
insertRestoreBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q)770 pllBoolean insertRestoreBIG (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q)
771 {
772 nodeptr r, s;
773
774 r = q->back;
775 s = p->back;
776
777 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
778
779 if(tr->thoroughInsertion)
780 {
781 hookup(p->next, q, tr->currentLZQ, numBranches);
782 hookup(p->next->next, r, tr->currentLZR, numBranches);
783 hookup(p, s, tr->currentLZS, numBranches);
784 }
785 else
786 {
787 double z[PLL_NUM_BRANCHES];
788 int i;
789
790 for(i = 0; i < numBranches; i++)
791 {
792 double zz;
793 zz = sqrt(q->z[i]);
794 if(zz < PLL_ZMIN)
795 zz = PLL_ZMIN;
796 if(zz > PLL_ZMAX)
797 zz = PLL_ZMAX;
798 z[i] = zz;
799 }
800
801 hookup(p->next, q, z, numBranches);
802 hookup(p->next->next, r, z, numBranches);
803 }
804
805 pllUpdatePartials(tr, pr,p, PLL_FALSE);
806
807 return PLL_TRUE;
808 }
809
810
restoreTopologyOnly(pllInstance * tr,bestlist * bt,int numBranches)811 static void restoreTopologyOnly(pllInstance *tr, bestlist *bt, int numBranches)
812 {
813 nodeptr p = tr->removeNode;
814 nodeptr q = tr->insertNode;
815 double qz[PLL_NUM_BRANCHES], pz[PLL_NUM_BRANCHES], p1z[PLL_NUM_BRANCHES], p2z[PLL_NUM_BRANCHES];
816 nodeptr p1, p2, r, s;
817 double currentLH = tr->likelihood;
818 int i;
819
820 p1 = p->next->back;
821 p2 = p->next->next->back;
822
823 //memcpy(p1z, p1->z, numBranches*sizeof(double));
824 //memcpy(p2z, p2->z, numBranches*sizeof(double));
825 //memcpy(qz, q->z, numBranches*sizeof(double));
826 //memcpy(pz, p->z, numBranches*sizeof(double));
827 for(i = 0; i < numBranches; i++)
828 {
829 p1z[i] = p1->z[i];
830 p2z[i] = p2->z[i];
831 }
832
833 hookup(p1, p2, tr->currentZQR, numBranches);
834
835 p->next->next->back = p->next->back = (node *) NULL;
836 for(i = 0; i < numBranches; i++)
837 {
838 qz[i] = q->z[i];
839 pz[i] = p->z[i];
840 }
841
842 r = q->back;
843 s = p->back;
844
845 if(tr->thoroughInsertion)
846 {
847 hookup(p->next, q, tr->currentLZQ, numBranches);
848 hookup(p->next->next, r, tr->currentLZR, numBranches);
849 hookup(p, s, tr->currentLZS, numBranches);
850 }
851 else
852 {
853 double z[PLL_NUM_BRANCHES];
854 for(i = 0; i < numBranches; i++)
855 {
856 z[i] = sqrt(q->z[i]);
857 if(z[i] < PLL_ZMIN)
858 z[i] = PLL_ZMIN;
859 if(z[i] > PLL_ZMAX)
860 z[i] = PLL_ZMAX;
861 }
862 hookup(p->next, q, z, numBranches);
863 hookup(p->next->next, r, z, numBranches);
864 }
865
866 tr->likelihood = tr->bestOfNode;
867
868 saveBestTree(bt, tr, numBranches);
869
870 tr->likelihood = currentLH;
871
872 hookup(q, r, qz, numBranches);
873
874 p->next->next->back = p->next->back = (nodeptr) NULL;
875
876 if(tr->thoroughInsertion)
877 hookup(p, s, pz, numBranches);
878
879 hookup(p->next, p1, p1z, numBranches);
880 hookup(p->next->next, p2, p2z, numBranches);
881 }
882
883 /** @brief Test the
884 */
testInsertBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q)885 pllBoolean testInsertBIG (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q)
886 {
887
888 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
889
890 double qz[PLL_NUM_BRANCHES], pz[PLL_NUM_BRANCHES];
891 nodeptr r;
892 double startLH = tr->endLH;
893 int i;
894
895 r = q->back;
896 for(i = 0; i < numBranches; i++)
897 {
898 qz[i] = q->z[i];
899 pz[i] = p->z[i];
900 }
901
902 if (! insertBIG(tr, pr, p, q)) return PLL_FALSE;
903
904 pllEvaluateLikelihood (tr, pr, p->next->next, PLL_FALSE, PLL_FALSE);
905
906 if(tr->likelihood > tr->bestOfNode)
907 {
908 tr->bestOfNode = tr->likelihood;
909 tr->insertNode = q;
910 tr->removeNode = p;
911 for(i = 0; i < numBranches; i++)
912 {
913 tr->currentZQR[i] = tr->zqr[i];
914 tr->currentLZR[i] = tr->lzr[i];
915 tr->currentLZQ[i] = tr->lzq[i];
916 tr->currentLZS[i] = tr->lzs[i];
917 }
918 }
919
920 if(tr->likelihood > tr->endLH)
921 {
922 tr->insertNode = q;
923 tr->removeNode = p;
924 for(i = 0; i < numBranches; i++)
925 tr->currentZQR[i] = tr->zqr[i];
926 tr->endLH = tr->likelihood;
927 }
928
929 /* reset the topology so that it is the same as it was before calling insertBIG */
930 hookup(q, r, qz, numBranches);
931
932 p->next->next->back = p->next->back = (nodeptr) NULL;
933
934 if(tr->thoroughInsertion)
935 {
936 nodeptr s = p->back;
937 hookup(p, s, pz, numBranches);
938 }
939
940 if((tr->doCutoff) && (tr->likelihood < startLH))
941 {
942 tr->lhAVG += (startLH - tr->likelihood);
943 tr->lhDEC++;
944 if((startLH - tr->likelihood) >= tr->lhCutoff)
945 return PLL_FALSE;
946 else
947 return PLL_TRUE;
948 }
949 else
950 return PLL_TRUE;
951 }
952
953
954 /** @brief Recursively traverse tree and test insertion
955
956 Recursively traverses the tree structure starting from node \a q and
957 tests the insertion of the component specified by leaf \a p at the edge
958 between \a q and \a q->back.
959
960 @param tr
961 PLL instance
962
963 @param pr
964 List of partitions
965 @param p
966 Leaf node of one tree component
967
968 @param q
969 Endpoint node of the edge to test the insertion
970
971 @param mintrav
972 Minimum radius around \a q to test the insertion
973
974 @param maxtrav
975 Maximum radius around \a q to test the insertion\
976 */
addTraverseBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q,int mintrav,int maxtrav)977 void addTraverseBIG(pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q, int mintrav, int maxtrav)
978 {
979 if (--mintrav <= 0)
980 {
981 if (! testInsertBIG(tr, pr, p, q)) return;
982
983 }
984
985 if ((!isTip(q->number, tr->mxtips)) && (--maxtrav > 0))
986 {
987 addTraverseBIG(tr, pr, p, q->next->back, mintrav, maxtrav);
988 addTraverseBIG(tr, pr, p, q->next->next->back, mintrav, maxtrav);
989 }
990 }
991
992
993
994
995 /** @brief Compute the best SPR movement
996
997 Compute all SPR moves starting from \a p in the space defined by \a mintrav and
998 \a maxtrav and store the best in the \a tr structure.
999
1000 @param tr
1001 PLL instancve
1002
1003 @param pr
1004 List of partitions
1005
1006 @param p
1007 Node from which to start the SPR moves testing
1008
1009 @param mintrav
1010 Minimum distance from \a p where to start testing SPRs
1011
1012 @param maxtrav
1013 Maximum distance from \a p where to test SPRs
1014
1015 @return
1016 0,1 or \b PLL_BADREAR
1017
1018 @todo
1019 fix the return value
1020 */
rearrangeBIG(pllInstance * tr,partitionList * pr,nodeptr p,int mintrav,int maxtrav)1021 int rearrangeBIG(pllInstance *tr, partitionList *pr, nodeptr p, int mintrav, int maxtrav)
1022 {
1023 double p1z[PLL_NUM_BRANCHES], p2z[PLL_NUM_BRANCHES], q1z[PLL_NUM_BRANCHES], q2z[PLL_NUM_BRANCHES];
1024 nodeptr p1, p2, q, q1, q2;
1025 int mintrav2, i;
1026 pllBoolean doP = PLL_TRUE, doQ = PLL_TRUE;
1027 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
1028
1029 if (maxtrav < 1 || mintrav > maxtrav) return (0);
1030 q = p->back;
1031
1032
1033
1034
1035 if (!isTip(p->number, tr->mxtips) && doP)
1036 {
1037 p1 = p->next->back;
1038 p2 = p->next->next->back;
1039
1040
1041 if(!isTip(p1->number, tr->mxtips) || !isTip(p2->number, tr->mxtips))
1042 {
1043 for(i = 0; i < numBranches; i++)
1044 {
1045 p1z[i] = p1->z[i];
1046 p2z[i] = p2->z[i];
1047 }
1048
1049 if (! removeNodeBIG(tr, pr, p, numBranches)) return PLL_BADREAR;
1050
1051 if (!isTip(p1->number, tr->mxtips))
1052 {
1053 addTraverseBIG(tr, pr, p, p1->next->back,
1054 mintrav, maxtrav);
1055
1056 addTraverseBIG(tr, pr, p, p1->next->next->back,
1057 mintrav, maxtrav);
1058 }
1059
1060 if (!isTip(p2->number, tr->mxtips))
1061 {
1062 addTraverseBIG(tr, pr, p, p2->next->back,
1063 mintrav, maxtrav);
1064 addTraverseBIG(tr, pr, p, p2->next->next->back,
1065 mintrav, maxtrav);
1066 }
1067
1068 hookup(p->next, p1, p1z, numBranches);
1069 hookup(p->next->next, p2, p2z, numBranches);
1070 pllUpdatePartials(tr, pr,p, PLL_FALSE);
1071 }
1072 }
1073
1074 if (!isTip(q->number, tr->mxtips) && maxtrav > 0 && doQ)
1075 {
1076 q1 = q->next->back;
1077 q2 = q->next->next->back;
1078
1079 /*if (((!q1->tip) && (!q1->next->back->tip || !q1->next->next->back->tip)) ||
1080 ((!q2->tip) && (!q2->next->back->tip || !q2->next->next->back->tip))) */
1081 if (
1082 (
1083 ! isTip(q1->number, tr->mxtips) &&
1084 (! isTip(q1->next->back->number, tr->mxtips) || ! isTip(q1->next->next->back->number, tr->mxtips))
1085 )
1086 ||
1087 (
1088 ! isTip(q2->number, tr->mxtips) &&
1089 (! isTip(q2->next->back->number, tr->mxtips) || ! isTip(q2->next->next->back->number, tr->mxtips))
1090 )
1091 )
1092 {
1093
1094 for(i = 0; i < numBranches; i++)
1095 {
1096 q1z[i] = q1->z[i];
1097 q2z[i] = q2->z[i];
1098 }
1099
1100 if (! removeNodeBIG(tr, pr, q, numBranches)) return PLL_BADREAR;
1101
1102 mintrav2 = mintrav > 2 ? mintrav : 2;
1103
1104 if (/*! q1->tip*/ !isTip(q1->number, tr->mxtips))
1105 {
1106 addTraverseBIG(tr, pr, q, q1->next->back,
1107 mintrav2 , maxtrav);
1108 addTraverseBIG(tr, pr, q, q1->next->next->back,
1109 mintrav2 , maxtrav);
1110 }
1111
1112 if (/*! q2->tip*/ ! isTip(q2->number, tr->mxtips))
1113 {
1114 addTraverseBIG(tr, pr, q, q2->next->back,
1115 mintrav2 , maxtrav);
1116 addTraverseBIG(tr, pr, q, q2->next->next->back,
1117 mintrav2 , maxtrav);
1118 }
1119
1120 hookup(q->next, q1, q1z, numBranches);
1121 hookup(q->next->next, q2, q2z, numBranches);
1122
1123 pllUpdatePartials(tr, pr,q, PLL_FALSE);
1124 }
1125 }
1126
1127 return 1;
1128 }
1129
1130
1131
1132
1133 /** @brief Perform an SPR move?
1134
1135 @param tr
1136 PLL instance
1137
1138 @param pr
1139 List of partitions
1140
1141 @param mintrav
1142
1143 @param maxtrav
1144
1145 @param adef
1146
1147 @param bt
1148
1149 @param iList
1150
1151 */
treeOptimizeRapid(pllInstance * tr,partitionList * pr,int mintrav,int maxtrav,bestlist * bt,infoList * iList)1152 double treeOptimizeRapid(pllInstance *tr, partitionList *pr, int mintrav, int maxtrav, bestlist *bt, infoList *iList)
1153 {
1154 int i, index,
1155 *perm = (int*)NULL;
1156
1157 nodeRectifier(tr);
1158
1159
1160
1161 if (maxtrav > tr->mxtips - 3)
1162 maxtrav = tr->mxtips - 3;
1163
1164
1165
1166 resetInfoList(iList);
1167
1168 resetBestTree(bt);
1169
1170 tr->startLH = tr->endLH = tr->likelihood;
1171
1172 if(tr->doCutoff)
1173 {
1174 if(tr->bigCutoff)
1175 {
1176 if(tr->itCount == 0)
1177 tr->lhCutoff = 0.5 * (tr->likelihood / -1000.0);
1178 else
1179 tr->lhCutoff = 0.5 * ((tr->lhAVG) / ((double)(tr->lhDEC)));
1180 }
1181 else
1182 {
1183 if(tr->itCount == 0)
1184 tr->lhCutoff = tr->likelihood / -1000.0;
1185 else
1186 tr->lhCutoff = (tr->lhAVG) / ((double)(tr->lhDEC));
1187 }
1188
1189 tr->itCount = tr->itCount + 1;
1190 tr->lhAVG = 0;
1191 tr->lhDEC = 0;
1192 }
1193
1194 /*
1195 printf("DoCutoff: %d\n", tr->doCutoff);
1196 printf("%d %f %f %f\n", tr->itCount, tr->lhAVG, tr->lhDEC, tr->lhCutoff);
1197
1198 printf("%d %d\n", mintrav, maxtrav);
1199 */
1200
1201 for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
1202 {
1203 tr->bestOfNode = PLL_UNLIKELY;
1204
1205 if(tr->permuteTreeoptimize)
1206 index = perm[i];
1207 else
1208 index = i;
1209
1210 if(rearrangeBIG(tr, pr, tr->nodep[index], mintrav, maxtrav))
1211 {
1212 if(tr->thoroughInsertion)
1213 {
1214 if(tr->endLH > tr->startLH)
1215 {
1216 /* commit the best SPR found by rearrangeBIG */
1217 restoreTreeFast(tr, pr);
1218 tr->startLH = tr->endLH = tr->likelihood;
1219 saveBestTree(bt, tr, pr->perGeneBranchLengths?pr->numberOfPartitions:1);
1220 }
1221 else
1222 {
1223 if(tr->bestOfNode != PLL_UNLIKELY)
1224 restoreTopologyOnly(tr, bt, pr->perGeneBranchLengths?pr->numberOfPartitions:1);
1225 }
1226 }
1227 else
1228 {
1229 insertInfoList(tr->nodep[index], tr->bestOfNode, iList);
1230 if(tr->endLH > tr->startLH)
1231 {
1232 restoreTreeFast(tr, pr);
1233 tr->startLH = tr->endLH = tr->likelihood;
1234 }
1235 }
1236 }
1237 }
1238
1239 if(!tr->thoroughInsertion)
1240 {
1241 tr->thoroughInsertion = PLL_TRUE;
1242
1243 for(i = 0; i < iList->valid; i++)
1244 {
1245 tr->bestOfNode = PLL_UNLIKELY;
1246
1247 if(rearrangeBIG(tr, pr, iList->list[i].node, mintrav, maxtrav))
1248 {
1249 if(tr->endLH > tr->startLH)
1250 {
1251 restoreTreeFast(tr, pr);
1252 tr->startLH = tr->endLH = tr->likelihood;
1253 saveBestTree(bt, tr, pr->perGeneBranchLengths?pr->numberOfPartitions:1);
1254 }
1255 else
1256 {
1257
1258 if(tr->bestOfNode != PLL_UNLIKELY)
1259 {
1260 restoreTopologyOnly(tr, bt, pr->perGeneBranchLengths?pr->numberOfPartitions:1);
1261 }
1262 }
1263 }
1264 }
1265
1266 tr->thoroughInsertion = PLL_FALSE;
1267 }
1268
1269 if(tr->permuteTreeoptimize)
1270 rax_free(perm);
1271
1272 return tr->startLH;
1273 }
1274
1275
1276
1277
testInsertRestoreBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q)1278 pllBoolean testInsertRestoreBIG (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q)
1279 {
1280 if(tr->thoroughInsertion)
1281 {
1282 if (! insertBIG(tr, pr, p, q)) return PLL_FALSE;
1283
1284 pllEvaluateLikelihood (tr, pr, p->next->next, PLL_FALSE, PLL_FALSE);
1285 }
1286 else
1287 {
1288 if (! insertRestoreBIG(tr, pr, p, q)) return PLL_FALSE;
1289
1290 {
1291 nodeptr x, y;
1292 x = p->next->next;
1293 y = p->back;
1294
1295 if(! isTip(x->number, tr->mxtips) && isTip(y->number, tr->mxtips))
1296 {
1297 while ((! x->x))
1298 {
1299 if (! (x->x))
1300 pllUpdatePartials(tr, pr,x, PLL_FALSE);
1301 }
1302 }
1303
1304 if(isTip(x->number, tr->mxtips) && !isTip(y->number, tr->mxtips))
1305 {
1306 while ((! y->x))
1307 {
1308 if (! (y->x))
1309 pllUpdatePartials(tr, pr,y, PLL_FALSE);
1310 }
1311 }
1312
1313 if(!isTip(x->number, tr->mxtips) && !isTip(y->number, tr->mxtips))
1314 {
1315 while ((! x->x) || (! y->x))
1316 {
1317 if (! (x->x))
1318 pllUpdatePartials(tr, pr,x, PLL_FALSE);
1319 if (! (y->x))
1320 pllUpdatePartials(tr, pr,y, PLL_FALSE);
1321 }
1322 }
1323
1324 }
1325
1326 tr->likelihood = tr->endLH;
1327 }
1328
1329 return PLL_TRUE;
1330 }
1331
restoreTreeFast(pllInstance * tr,partitionList * pr)1332 void restoreTreeFast(pllInstance *tr, partitionList *pr)
1333 {
1334 removeNodeRestoreBIG(tr, pr, tr->removeNode);
1335 testInsertRestoreBIG(tr, pr, tr->removeNode, tr->insertNode);
1336 }
1337
1338 /*
1339 static void myfwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream)
1340 {
1341 size_t
1342 bytes_written = fwrite(ptr, size, nmemb, stream);
1343
1344 assert(bytes_written == nmemb);
1345 }
1346
1347 static void myfread(void *ptr, size_t size, size_t nmemb, FILE *stream)
1348 {
1349 size_t
1350 bytes_read;
1351
1352 bytes_read = fread(ptr, size, nmemb, stream);
1353
1354 assert(bytes_read == nmemb);
1355 }
1356
1357 static void readTree(pllInstance *tr, partitionList *pr, FILE *f)
1358 {
1359 int
1360 nodeNumber,
1361 x = tr->mxtips + 3 * (tr->mxtips - 1);
1362
1363 nodeptr
1364 startAddress;
1365
1366 myfread(&nodeNumber, sizeof(int), 1, f);
1367
1368 tr->start = tr->nodep[nodeNumber];
1369
1370
1371 myfread(&startAddress, sizeof(nodeptr), 1, f);
1372
1373 myfread(tr->nodeBaseAddress, sizeof(node), x, f);
1374
1375 {
1376 int i;
1377
1378 size_t
1379 offset;
1380
1381 pllBoolean
1382 addIt;
1383
1384 if(startAddress > tr->nodeBaseAddress)
1385 {
1386 addIt = PLL_FALSE;
1387 offset = (size_t)startAddress - (size_t)tr->nodeBaseAddress;
1388 }
1389 else
1390 {
1391 addIt = PLL_TRUE;
1392 offset = (size_t)tr->nodeBaseAddress - (size_t)startAddress;
1393 }
1394
1395 for(i = 0; i < x; i++)
1396 {
1397 if(addIt)
1398 {
1399 tr->nodeBaseAddress[i].next = (nodeptr)((size_t)tr->nodeBaseAddress[i].next + offset);
1400 tr->nodeBaseAddress[i].back = (nodeptr)((size_t)tr->nodeBaseAddress[i].back + offset);
1401 }
1402 else
1403 {
1404
1405 tr->nodeBaseAddress[i].next = (nodeptr)((size_t)tr->nodeBaseAddress[i].next - offset);
1406 tr->nodeBaseAddress[i].back = (nodeptr)((size_t)tr->nodeBaseAddress[i].back - offset);
1407 }
1408 }
1409
1410 }
1411
1412 pllEvaluateLikelihood (tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
1413
1414 printBothOpen("RAxML Restart with likelihood: %1.50f\n", tr->likelihood);
1415 }
1416
1417 static void readCheckpoint(pllInstance *tr, partitionList *pr)
1418 {
1419 int
1420 restartErrors = 0,
1421 model;
1422
1423 FILE
1424 *f = myfopen(binaryCheckpointInputName, "r");
1425 */
1426 /* cdta */
1427 /*
1428 myfread(&(tr->ckp), sizeof(checkPointState), 1, f);
1429
1430
1431
1432 if(tr->ckp.searchConvergenceCriterion != tr->searchConvergenceCriterion)
1433 {
1434 printf("restart error, you are trying to re-start a run where the ML search criterion was turned %s\n", (tr->ckp.searchConvergenceCriterion)?"ON":"OFF");
1435 restartErrors++;
1436 }
1437
1438 if(tr->ckp.rateHetModel != tr->rateHetModel)
1439 {
1440 printf("restart error, you are trying to re-start a run with a different model of rate heterogeneity, the checkpoint was obtained under: %s\n", (tr->ckp.rateHetModel == PLL_GAMMA)?"GAMMA":"PSR");
1441 restartErrors++;
1442 }
1443
1444 if(tr->ckp.maxCategories != tr->maxCategories)
1445 {
1446 printf("restart error, you are trying to re-start a run with %d per-site rate categories, the checkpoint was obtained with: %d\n", tr->maxCategories, tr->ckp.maxCategories);
1447 restartErrors++;
1448 }
1449
1450 if(tr->ckp.NumberOfModels != pr->numberOfPartitions)
1451 {
1452 printf("restart error, you are trying to re-start a run with %d partitions, the checkpoint was obtained with: %d partitions\n", (int)pr->numberOfPartitions, tr->ckp.NumberOfModels);
1453 restartErrors++;
1454 }
1455
1456 if(tr->ckp.numBranches != pr->perGeneBranchLengths?pr->numberOfPartitions:1)
1457 {
1458 printf("restart error, you are trying to re-start a run where independent per-site branch length estimates were turned %s\n", (tr->ckp.numBranches > 1)?"ON":"OFF");
1459 restartErrors++;
1460 }
1461
1462 if(tr->ckp.originalCrunchedLength != tr->originalCrunchedLength)
1463 {
1464 printf("restart error, you are trying to re-start a run with %d site patterns, the checkpoint was obtained with: %d site patterns\n", tr->ckp.originalCrunchedLength, tr->originalCrunchedLength);
1465 restartErrors++;
1466 }
1467
1468 if(tr->ckp.mxtips != tr->mxtips)
1469 {
1470 printf("restart error, you are trying to re-start a run with %d taxa, the checkpoint was obtained with: %d taxa\n", tr->mxtips, tr->ckp.mxtips);
1471 restartErrors++;
1472 }
1473
1474 if(strcmp(tr->ckp.seq_file, seq_file) != 0)
1475 {
1476 printf("restart error, you are trying to re-start from alignemnt file %s, the checkpoint was obtained with file: %s\n", tr->ckp.seq_file, seq_file);
1477 restartErrors++;
1478 }
1479
1480 printf("REstart errors: %d\n", restartErrors);
1481
1482 if(restartErrors > 0)
1483 {
1484 printf("User induced errors with the restart from checkpoint, exiting ...\n");
1485
1486 if(restartErrors > 4)
1487 printf(" ... maybe you should do field work instead of trying to use a computer ...\n");
1488 if(restartErrors > 6)
1489 printf(" ... kala eisai telios ilithios;\n");
1490
1491 exit(-1);
1492 }
1493
1494 tr->ntips = tr->mxtips;
1495
1496 tr->startLH = tr->ckp.tr_startLH;
1497 tr->endLH = tr->ckp.tr_endLH;
1498 tr->likelihood = tr->ckp.tr_likelihood;
1499 tr->bestOfNode = tr->ckp.tr_bestOfNode;
1500
1501 tr->lhCutoff = tr->ckp.tr_lhCutoff;
1502 tr->lhAVG = tr->ckp.tr_lhAVG;
1503 tr->lhDEC = tr->ckp.tr_lhDEC;
1504 tr->itCount = tr->ckp.tr_itCount;
1505 tr->thoroughInsertion = tr->ckp.tr_thoroughInsertion;
1506
1507
1508
1509 accumulatedTime = tr->ckp.accumulatedTime;
1510 */
1511 /* printf("Accumulated time so far: %f\n", accumulatedTime); */
1512 /*
1513 tr->optimizeRateCategoryInvocations = tr->ckp.tr_optimizeRateCategoryInvocations;
1514
1515
1516 myfread(tr->tree0, sizeof(char), tr->treeStringLength, f);
1517 myfread(tr->tree1, sizeof(char), tr->treeStringLength, f);
1518
1519 if(tr->searchConvergenceCriterion)
1520 {
1521 int bCounter = 0;
1522
1523 if((tr->ckp.state == PLL_FAST_SPRS && tr->ckp.fastIterations > 0) ||
1524 (tr->ckp.state == PLL_SLOW_SPRS && tr->ckp.thoroughIterations > 0))
1525 {
1526
1527 #ifdef _DEBUG_CHECKPOINTING
1528 printf("parsing Tree 0\n");
1529 #endif
1530
1531 treeReadTopologyString(tr->tree0, tr);
1532
1533 bitVectorInitravSpecial(tr->bitVectors, tr->nodep[1]->back, tr->mxtips, tr->vLength, tr->h, 0, PLL_BIPARTITIONS_RF, (branchInfo *)NULL,
1534 &bCounter, 1, PLL_FALSE, PLL_FALSE, tr->threadID);
1535
1536 assert(bCounter == tr->mxtips - 3);
1537 }
1538
1539 bCounter = 0;
1540
1541 if((tr->ckp.state == PLL_FAST_SPRS && tr->ckp.fastIterations > 1) ||
1542 (tr->ckp.state == PLL_SLOW_SPRS && tr->ckp.thoroughIterations > 1))
1543 {
1544
1545 #ifdef _DEBUG_CHECKPOINTING
1546 printf("parsing Tree 1\n");
1547 #endif
1548
1549 treeReadTopologyString(tr->tree1, tr);
1550
1551 bitVectorInitravSpecial(tr->bitVectors, tr->nodep[1]->back, tr->mxtips, tr->vLength, tr->h, 1, PLL_BIPARTITIONS_RF, (branchInfo *)NULL,
1552 &bCounter, 1, PLL_FALSE, PLL_FALSE, tr->threadID);
1553
1554 assert(bCounter == tr->mxtips - 3);
1555 }
1556 }
1557
1558 myfread(tr->rateCategory, sizeof(int), tr->originalCrunchedLength, f);
1559 myfread(tr->patrat, sizeof(double), tr->originalCrunchedLength, f);
1560 myfread(tr->patratStored, sizeof(double), tr->originalCrunchedLength, f);
1561
1562 */
1563 /* need to read this as well in checkpoints, otherwise the branch lengths
1564 in the output tree files will be wrong, not the internal branch lengths though */
1565 /*
1566 //TODO: Same problem as writing the checkpoint
1567 //myfread(tr->fracchanges, sizeof(double), pr->numberOfPartitions, f);
1568 myfread(&(tr->fracchange), sizeof(double), 1, f);
1569 */
1570 /* pInfo */
1571 /*
1572 for(model = 0; model < pr->numberOfPartitions; model++)
1573 {
1574 int
1575 dataType = pr->partitionData[model]->dataType;
1576
1577 myfread(&(pr->partitionData[model]->numberOfCategories), sizeof(int), 1, f);
1578 myfread(pr->partitionData[model]->perSiteRates, sizeof(double), tr->maxCategories, f);
1579 myfread(pr->partitionData[model]->EIGN, sizeof(double), pLengths[dataType].eignLength, f);
1580 myfread(pr->partitionData[model]->EV, sizeof(double), pLengths[dataType].evLength, f);
1581 myfread(pr->partitionData[model]->EI, sizeof(double), pLengths[dataType].eiLength, f);
1582
1583 myfread(pr->partitionData[model]->frequencies, sizeof(double), pLengths[dataType].frequenciesLength, f);
1584 myfread(pr->partitionData[model]->tipVector, sizeof(double), pLengths[dataType].tipVectorLength, f);
1585 myfread(pr->partitionData[model]->substRates, sizeof(double), pLengths[dataType].substRatesLength, f);
1586 myfread(&(pr->partitionData[model]->alpha), sizeof(double), 1, f);
1587
1588 if(pr->partitionData[model]->protModels == PLL_LG4M || pr->partitionData[model]->protModels == PLL_LG4X)
1589 {
1590 int
1591 k;
1592
1593 for(k = 0; k < 4; k++)
1594 {
1595 myfread(pr->partitionData[model]->EIGN_LG4[k], sizeof(double), pLengths[dataType].eignLength, f);
1596 myfread(pr->partitionData[model]->EV_LG4[k], sizeof(double), pLengths[dataType].evLength, f);
1597 myfread(pr->partitionData[model]->EI_LG4[k], sizeof(double), pLengths[dataType].eiLength, f);
1598 myfread(pr->partitionData[model]->frequencies_LG4[k], sizeof(double), pLengths[dataType].frequenciesLength, f);
1599 myfread(pr->partitionData[model]->tipVector_LG4[k], sizeof(double), pLengths[dataType].tipVectorLength, f);
1600 myfread(pr->partitionData[model]->substRates_LG4[k], sizeof(double), pLengths[dataType].substRatesLength, f);
1601 }
1602 }
1603
1604 pllMakeGammaCats(pr->partitionData[model]->alpha, pr->partitionData[model]->gammaRates, 4, tr->useMedian);
1605 }
1606
1607 #if (defined(_FINE_GRAIN_MPI) || defined(_USE_PTHREADS))
1608 pllMasterBarrier (tr, pr, PLL_THREAD_COPY_INIT_MODEL);
1609 #endif
1610
1611 updatePerSiteRates(tr, pr, PLL_FALSE);
1612
1613 readTree(tr, pr, f);
1614
1615 fclose(f);
1616
1617 }
1618
1619 void restart(pllInstance *tr, partitionList *pr)
1620 {
1621 readCheckpoint(tr, pr);
1622
1623 switch(tr->ckp.state)
1624 {
1625 case PLL_REARR_SETTING:
1626 break;
1627 case PLL_FAST_SPRS:
1628 break;
1629 case PLL_SLOW_SPRS:
1630 break;
1631 default:
1632 assert(0);
1633 }
1634 }
1635 */
1636
1637 /* The number of maximum smoothing iterations is given explicitely */
1638 /** @brief Optimize branch lenghts and evaluate likelihood of topology
1639
1640 Optimize the branch lengths \a maxSmoothIterations times and evaluate
1641 the likelihood of tree. The resulting likelihood is placed in
1642 \a tr->likelihood
1643
1644 @param tr
1645 The PLL instance
1646
1647 @param pr
1648 List of partitions
1649
1650 @param maxSmoothIterations
1651 Number of times to optimize branch lengths
1652 */
1653 void
pllOptimizeBranchLengths(pllInstance * tr,partitionList * pr,int maxSmoothIterations)1654 pllOptimizeBranchLengths (pllInstance *tr, partitionList *pr, int maxSmoothIterations) /* Evaluate a user tree */
1655 {
1656 smoothTree(tr, pr, maxSmoothIterations); /* former (32 * smoothFactor) */
1657
1658 pllEvaluateLikelihood (tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
1659 }
1660
1661 /** @brief Perform an NNI move
1662
1663 Modify the tree topology of instance \a tr by performing an NNI (Neighbour Neighbor
1664 Interchange) move at node \a p. Let \a q be \a p->back. If \a swap is set to \b PLL_NNI_P_NEXT
1665 then the subtrees rooted at \a p->next->back and \a q->next->back will be swapped. Otherwise,
1666 if \a swap is set to \b PLL_NNI_P_NEXTNEXT then the subtrees rooted at \a p->next->next->back and
1667 \a q->next->back are swapped. For clarity, see the illustration.
1668
1669 @param tr
1670 PLL instance
1671
1672 @param p
1673 Node to use as origin for performing NNI
1674
1675 @param swap
1676 Which node to use for the NNI move. \b PLL_NNI_P_NEXT uses node p->next while \b PLL_NNI_P_NEXTNEXT uses p->next->next
1677
1678 @return
1679 In case of success \b PLL_TRUE, otherwise \b PLL_FALSE
1680
1681 @todo
1682 Started error checking here. Instead of checking the errors in the specified way, implement a variadic
1683 function where we pass the results of each check and the error code we want to assign if there is at
1684 least one negative result
1685
1686 @image html nni.png "In case \a swap is set to \b PLL_NNI_P_NEXT then the dashed red edge between \a p and \a r is removed and the blue edges are created. If \a swap is set to \b PLL_INIT_P_NEXTNEXT then the dashed red edge between \a p and \a s is removed and the green edges are created. In both cases the black dashed edge is removed"
1687 */
pllTopologyPerformNNI(pllInstance * tr,nodeptr p,int swap)1688 int pllTopologyPerformNNI(pllInstance * tr, nodeptr p, int swap)
1689 {
1690 nodeptr q, r;
1691
1692 q = p->back;
1693 if (isTip(q->number, tr->mxtips))
1694 {
1695 errno = PLL_NNI_Q_TIP;
1696 return (PLL_FALSE);
1697 }
1698 if (isTip(p->number, tr->mxtips))
1699 {
1700 errno = PLL_NNI_P_TIP;
1701 return (PLL_FALSE);
1702 }
1703 assert(!isTip(q->number, tr->mxtips));
1704 assert(!isTip(p->number, tr->mxtips));
1705
1706
1707 if(swap == PLL_NNI_P_NEXT)
1708 {
1709 r = p->next->back;
1710 hookupFull(p->next, q->next->back, q->next->z);
1711 hookupFull(q->next, r, p->next->z);
1712 }
1713 else
1714 {
1715 r = p->next->next->back;
1716 hookupFull(p->next->next, q->next->back, q->next->z);
1717 hookupFull(q->next, r, p->next->next->z);
1718 }
1719
1720 return PLL_TRUE;
1721 }
1722
1723 /** @brief Compares 2 NNI moves */
cmp_nni(const void * nni1,const void * nni2)1724 static int cmp_nni(const void* nni1, const void* nni2) {
1725 nniMove* myNNI1 = (nniMove*) nni1;
1726 nniMove* myNNI2 = (nniMove*) nni2;
1727 return (int) (1000000.f * myNNI1->deltaLH - 1000000.f * myNNI2->deltaLH);
1728 }
1729
1730 /** @brief Gets the best NNI move for a branch
1731
1732 @param tr
1733 PLL instance
1734
1735 @param pr
1736 List of partitions
1737
1738 @param p
1739 Node to use as origin for performing NNI
1740
1741 @param curLH
1742 The current likelihood
1743
1744 @return
1745 The best NNI move
1746
1747 */
getBestNNIForBran(pllInstance * tr,partitionList * pr,nodeptr p,double curLH)1748 nniMove getBestNNIForBran(pllInstance* tr, partitionList *pr, nodeptr p,
1749 double curLH) {
1750 nodeptr q = p->back;
1751 assert( ! isTip(p->number, tr->mxtips));
1752 assert( ! isTip(q->number, tr->mxtips));
1753 #ifdef _DEBUG_NNI
1754 pllTreeToNewick(tr->tree_string, tr, tr->start->back, TRUE, FALSE, 0, 0, 0, SUMMARIZE_LH, 0,0);
1755 fprintf(stderr, "%s\n", tr->tree_string);
1756 #endif
1757
1758 /* Backup the current branch length */
1759 double z0[PLL_NUM_BRANCHES];
1760 int i;
1761 for (i = 0; i < pr->numberOfPartitions; i++) {
1762 z0[i] = p->z[i];
1763 }
1764 #ifdef _DEBUG_NNI
1765 double lhOld = tr->likelihood;
1766 printf("lhOld: %f \n", lhOld);
1767 #endif
1768 double lh0 = curLH;
1769
1770
1771 #ifdef _DEBUG_NNI
1772 printf("lh0: %f \n", lh0);
1773 #endif
1774 nniMove nni0; // nni0 means no NNI move is done
1775 nni0.p = p;
1776 nni0.nniType = 0;
1777 nni0.deltaLH = 0;
1778 for (i = 0; i < pr->numberOfPartitions; i++) {
1779 nni0.z[i] = p->z[i];
1780 }
1781
1782 /* Save the scaling factor */
1783 // Now try to do an NNI move of type 1
1784 pllTopologyPerformNNI(tr, p, PLL_NNI_P_NEXT);
1785 double lh1 = tr->likelihood;
1786 /* Update branch lengths */
1787 pllUpdatePartials(tr, pr, p, PLL_FALSE);
1788 pllUpdatePartials(tr, pr, q, PLL_FALSE);
1789 update(tr, pr, p);
1790 pllEvaluateLikelihood (tr, pr, p, PLL_FALSE, PLL_FALSE);
1791
1792 nniMove nni1;
1793 nni1.p = p;
1794 nni1.nniType = 1;
1795 // Store the optimized und unoptimized central branch length
1796 for (i = 0; i < pr->numberOfPartitions; i++) {
1797 nni1.z[i] = p->z[i];
1798 nni1.z0[i] = z0[i];
1799 }
1800 nni1.likelihood = lh1;
1801 nni1.deltaLH = lh1 - lh0;
1802 #ifdef _DEBUG_NNI
1803 printf("Delta likelihood of the 1.NNI move: %f\n", nni1.deltaLH);
1804 #endif
1805
1806 /* Restore previous NNI move */
1807 pllTopologyPerformNNI(tr, p, PLL_NNI_P_NEXT);
1808 /* Restore the old branch length */
1809 for (i = 0; i < pr->numberOfPartitions; i++) {
1810 p->z[i] = z0[i];
1811 p->back->z[i] = z0[i];
1812 }
1813
1814 #ifdef _DEBUG_NNI
1815 printf("Restore topology\n");
1816 pllTreeToNewick(tr->tree_string, tr, tr->start->back, TRUE, FALSE, 0, 0, 0, SUMMARIZE_LH, 0,0);
1817 fprintf(stderr, "%s\n", tr->tree_string);
1818 pllEvaluateLikelihood (tr, tr->start, TRUE);
1819 printf("Likelihood after restoring from NNI 1: %f\n", tr->likelihood);
1820 #endif
1821 /* Try to do an NNI move of type 2 */
1822 pllTopologyPerformNNI(tr, p, 2);
1823 double lh2 = tr->likelihood;
1824 /* Update branch lengths */
1825 pllUpdatePartials(tr, pr, p, PLL_FALSE);
1826 pllUpdatePartials(tr, pr, q, PLL_FALSE);
1827 update(tr, pr, p);
1828 pllEvaluateLikelihood (tr, pr, p, PLL_FALSE, PLL_FALSE);
1829
1830 // Create the nniMove struct to store this move
1831 nniMove nni2;
1832 nni2.p = p;
1833 nni2.nniType = 2;
1834
1835 // Store the optimized and unoptimized central branch length
1836 for (i = 0; i < pr->numberOfPartitions; i++) {
1837 nni2.z[i] = p->z[i];
1838 nni2.z0[i] = z0[i];
1839 }
1840 nni2.likelihood = lh2;
1841 nni2.deltaLH = lh2 - lh0;
1842 #ifdef _DEBUG_NNI
1843 printf("Delta likelihood of the 2.NNI move: %f\n", nni2.deltaLH);
1844 #endif
1845
1846 /* Restore previous NNI move */
1847 pllTopologyPerformNNI(tr, p, 2);
1848 pllUpdatePartials(tr, pr, p, PLL_FALSE);
1849 pllUpdatePartials(tr, pr, p->back, PLL_FALSE);
1850 /* Restore the old branch length */
1851 for (i = 0; i < pr->numberOfPartitions; i++) {
1852 p->z[i] = z0[i];
1853 p->back->z[i] = z0[i];
1854 }
1855 if (nni1.deltaLH > 0 && nni1.deltaLH >= nni2.deltaLH) {
1856 return nni1;
1857 } else if (nni1.deltaLH > 0 && nni1.deltaLH < nni2.deltaLH) {
1858 return nni2;
1859 } else if (nni1.deltaLH < 0 && nni2.deltaLH > 0) {
1860 return nni2;
1861 } else {
1862 return nni0;
1863 }
1864 }
1865
1866 /** @brief ??? Not sure */
evalNNIForSubtree(pllInstance * tr,partitionList * pr,nodeptr p,nniMove * nniList,int * cnt,int * cnt_nni,double curLH)1867 void evalNNIForSubtree(pllInstance* tr, partitionList *pr, nodeptr p,
1868 nniMove* nniList, int* cnt, int* cnt_nni, double curLH) {
1869 if (!isTip(p->number, tr->mxtips)) {
1870 nniList[*cnt] = getBestNNIForBran(tr, pr, p, curLH);
1871 if (nniList[*cnt].deltaLH != 0.0) {
1872 *cnt_nni = *cnt_nni + 1;
1873 }
1874 *cnt = *cnt + 1;
1875 nodeptr q = p->next;
1876 while (q != p) {
1877 evalNNIForSubtree(tr, pr, q->back, nniList, cnt, cnt_nni, curLH);
1878 q = q->next;
1879 }
1880 }
1881 }
1882
1883 /** @brief Perform an NNI search
1884
1885 Modify the tree topology of instance and model parameters \a tr by performing a NNI (Neighbour Neighbor
1886 Interchange) moves \a p.
1887
1888 @param tr
1889 PLL instance
1890
1891 @param pr
1892 List of partitions
1893
1894 @param estimateModel
1895 Determine wheter the model parameters should be optimized
1896
1897 @return
1898 In case of success \b PLL_TRUE, otherwise \b PLL_FALSE
1899
1900 */
pllNniSearch(pllInstance * tr,partitionList * pr,int estimateModel)1901 int pllNniSearch(pllInstance * tr, partitionList *pr, int estimateModel) {
1902
1903 double curScore = tr->likelihood;
1904
1905 /* Initialize the NNI list */
1906 nniMove* nniList = (nniMove*) malloc((tr->mxtips - 3) * sizeof(nniMove));
1907 int i;
1908 /* fill up the NNI list */
1909 nodeptr p = tr->start->back;
1910 nodeptr q = p->next;
1911 int cnt = 0; // number of visited internal branches during NNI evaluation
1912 int cnt_nni = 0; // number of positive NNI found
1913 while (q != p) {
1914 evalNNIForSubtree(tr, pr, q->back, nniList, &cnt, &cnt_nni, curScore);
1915 q = q->next;
1916 }
1917 if (cnt_nni == 0)
1918 return 0.0;
1919
1920 nniMove* impNNIList = (nniMove*) malloc(cnt_nni * sizeof(nniMove));
1921 int j = 0;
1922 for (i = 0; i < tr->mxtips - 3; i++) {
1923 if (nniList[i].deltaLH > 0.0) {
1924 impNNIList[j] = nniList[i];
1925 j++;
1926 }
1927 }
1928 // sort impNNIList
1929 qsort(impNNIList, cnt_nni, sizeof(nniMove), cmp_nni);
1930
1931 // creating a list of non-conflicting positive NNI
1932 nniMove* nonConfNNIList = (nniMove*) calloc(cnt_nni, sizeof(nniMove));
1933
1934 // the best NNI will always be taken
1935 nonConfNNIList[0] = impNNIList[cnt_nni - 1];
1936
1937 // Filter out conflicting NNI
1938 int numNonConflictNNI = 1; // size of the non-conflicting NNI list;
1939 int k;
1940 for (k = cnt_nni - 2; k >= 0; k--) {
1941 int conflict = PLL_FALSE;
1942 int j;
1943 for (j = 0; j < numNonConflictNNI; j++) {
1944 if (impNNIList[k].p->number == nonConfNNIList[j].p->number
1945 || impNNIList[k].p->number
1946 == nonConfNNIList[j].p->back->number) {
1947 conflict = PLL_TRUE;
1948 break;
1949 }
1950 }
1951 if (conflict) {
1952 continue;
1953 } else {
1954 nonConfNNIList[numNonConflictNNI] = impNNIList[k];
1955 numNonConflictNNI++;
1956 }
1957 }
1958
1959 // Applying non-conflicting NNI moves
1960 double delta = 1.0; // portion of NNI moves to apply
1961 int notImproved;
1962 do {
1963 notImproved = PLL_FALSE;
1964 int numNNI2Apply = ceil(numNonConflictNNI * delta);
1965 for (i = 0; i < numNNI2Apply; i++) {
1966 // Just do the topological change
1967 pllTopologyPerformNNI(tr, nonConfNNIList[i].p, nonConfNNIList[i].nniType);
1968 pllUpdatePartials(tr, pr, nonConfNNIList[i].p, PLL_FALSE);
1969 pllUpdatePartials(tr, pr, nonConfNNIList[i].p->back, PLL_FALSE);
1970 // Apply the store branch length
1971 int j;
1972 for (j = 0; j < pr->numberOfPartitions; j++) {
1973 nonConfNNIList[i].p->z[j] = nonConfNNIList[i].z[j];
1974 nonConfNNIList[i].p->back->z[j] = nonConfNNIList[i].z[j];
1975 }
1976 }
1977 // Re-optimize all branches
1978 smoothTree(tr, pr, 2);
1979 pllEvaluateLikelihood (tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
1980 if (estimateModel) {
1981 modOpt(tr, pr, 0.1);
1982 }
1983 pllEvaluateLikelihood (tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
1984 if (tr->likelihood < curScore) {
1985 #ifdef _DEBUG_NNI
1986 printf("Tree likelihood gets worse after applying NNI\n");
1987 printf("curScore = %30.20f\n", curScore);
1988 printf("newScore = %30.20f\n", tr->likelihood);
1989 printf("Rolling back the tree\n");
1990 #endif
1991 for (i = 0; i < numNNI2Apply; i++) {
1992 pllTopologyPerformNNI(tr, nonConfNNIList[i].p, nonConfNNIList[i].nniType);
1993 // Restore the branch length
1994 int j;
1995 for (j = 0; j < pr->numberOfPartitions; j++) {
1996 nonConfNNIList[i].p->z[j] = nonConfNNIList[i].z0[j];
1997 nonConfNNIList[i].p->back->z[j] = nonConfNNIList[i].z0[j];
1998 }
1999 }
2000 pllEvaluateLikelihood (tr, pr, tr->start, PLL_FALSE, PLL_FALSE);
2001 #ifdef _DEBUG_NNI
2002 printf("Tree likelihood after rolling back = %f \n",
2003 tr->likelihood);
2004 #endif
2005 notImproved = PLL_TRUE & (numNNI2Apply > 1);
2006 delta = delta * 0.5;
2007 }
2008 } while (notImproved);
2009 free(nniList);
2010 free(impNNIList);
2011 free(nonConfNNIList);
2012
2013 return PLL_TRUE;
2014 }
2015
2016
2017 /** @defgroup rearrangementGroup Topological rearrangements
2018
2019 This set of functions handles the rearrangement of the tree topology
2020 */
2021
2022
2023 /** @ingroup rearrangementGroup
2024 @brief Create a list for storing topology rearrangements
2025
2026 Allocates space and initializes a structure that will hold information
2027 of \a max topological rearrangements
2028
2029 @param max
2030 Maximum number of elements that the structure should hold
2031
2032 @note This should be called for creating a storage space (list) for
2033 routines such as ::pllRearrangeSearch which compute the best NNI/PR/TBR rearrangements.
2034 */
pllCreateRearrangeList(int max)2035 pllRearrangeList * pllCreateRearrangeList (int max)
2036 {
2037 pllRearrangeList * bl;
2038
2039 bl = (pllRearrangeList *) malloc (sizeof (pllRearrangeList));
2040
2041 bl->max_entries = max;
2042 bl->entries = 0;
2043 bl->rearr = (pllRearrangeInfo *) malloc (max * sizeof (pllRearrangeInfo));
2044
2045 return bl;
2046 }
2047
2048 /** @ingroup rearrangementGroup
2049 @brief Deallocator for topology rearrangements list
2050
2051 Call this to destroy (deallocate) the memory taken by the \a bestList which holds
2052 topological rearrangements
2053
2054 @param bestList
2055 Pointer to the list to be deallocated
2056 */
pllDestroyRearrangeList(pllRearrangeList ** bestList)2057 void pllDestroyRearrangeList (pllRearrangeList ** bestList)
2058 {
2059 pllRearrangeList * bl;
2060
2061 bl = *bestList;
2062
2063 rax_free (bl->rearr);
2064 rax_free (bl);
2065
2066 *bestList = NULL;
2067 }
2068
2069
2070 /** @ingroup rearrangementGroup
2071 @brief Store a rearrangement move to the list of best rearrangement moves
2072
2073 Checks if the likelihood yielded by the rearrangement move described in \a rearr
2074 is better than any in the sorted list \a bestList. If it is, or
2075 if there is still space in \a bestList, the info about the
2076 move is inserted in the list.
2077
2078 @param bestList
2079 The list of information about the best rearrangement moves
2080
2081 @param rearr
2082 Info about the current rearrangement move
2083
2084 @return
2085 Returns \b PLL_FALSE if the rearrangement move doesn't make it in the list, otherwise \b PLL_TRUE
2086 */
pllStoreRearrangement(pllRearrangeList * bestList,pllRearrangeInfo * rearr)2087 static int pllStoreRearrangement (pllRearrangeList * bestList, pllRearrangeInfo * rearr)
2088 {
2089 /* naive implementation of saving rearrangement moves */
2090 int i;
2091
2092 for (i = 0; i < bestList->entries; ++ i)
2093 {
2094 /* Does the new rearrangement yield a better likelihood that the current in the list */
2095 if (rearr->likelihood > bestList->rearr[i].likelihood)
2096 {
2097 /* is there enough space in the array ? */
2098 if (bestList->entries < bestList->max_entries)
2099 {
2100 /* slide the entries to the right and overwrite the i-th element with the new item */
2101 memmove (&(bestList->rearr[i + 1]), &(bestList->rearr[i]), (bestList->entries - i ) * sizeof (pllRearrangeInfo));
2102 ++ bestList->entries;
2103 }
2104 else
2105 {
2106 memmove (&(bestList->rearr[i + 1]), &(bestList->rearr[i]), (bestList->entries - i - 1 ) * sizeof (pllRearrangeInfo));
2107 }
2108 memcpy (&(bestList->rearr[i]), rearr, sizeof (pllRearrangeInfo));
2109 return (PLL_TRUE);
2110 }
2111 }
2112 if (bestList->entries < bestList->max_entries)
2113 {
2114 memcpy (&(bestList->rearr[bestList->entries]), rearr, sizeof (pllRearrangeInfo));
2115 ++ bestList->entries;
2116 return (PLL_TRUE);
2117 }
2118
2119 return (PLL_FALSE);
2120 }
2121
2122 /** @ingroup rearrangementGroup
2123 @brief Internal function for testing and saving an SPR move
2124
2125 Checks the likelihood of the placement of the pruned subtree specified by \a p
2126 to node \a q. If the likelihood is better than some in the sorted list
2127 \a bestList, or if there is still free space in \a bestList, then the SPR
2128 move is recorded (in \a bestList)
2129
2130 @param tr
2131 PLL instance
2132
2133 @param pr
2134 List of partitions
2135
2136 @param p
2137 Root of the subtree that is to be pruned
2138
2139 @param q
2140 Where to place the pruned subtree (between \a q and \a q->back
2141
2142 @param bestList
2143 Where to store the SPR move
2144
2145 @note Internal function which is not part of the PLL API and therefore should not be
2146 called by the user
2147
2148 @return
2149 */
2150 static int
pllTestInsertBIG(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q,pllRearrangeList * bestList)2151 pllTestInsertBIG (pllInstance * tr, partitionList * pr, nodeptr p, nodeptr q, pllRearrangeList * bestList)
2152 {
2153 int numBranches = pr->perGeneBranchLengths?pr->numberOfPartitions:1;
2154 pllRearrangeInfo rearr;
2155
2156 double qz[PLL_NUM_BRANCHES], pz[PLL_NUM_BRANCHES];
2157 nodeptr r;
2158 //double startLH = tr->endLH;
2159 int i;
2160
2161 r = q->back;
2162 for(i = 0; i < numBranches; i++)
2163 {
2164 qz[i] = q->z[i];
2165 pz[i] = p->z[i];
2166 }
2167
2168 if (! insertBIG(tr, pr, p, q)) return PLL_FALSE;
2169
2170 pllEvaluateLikelihood (tr, pr, p->next->next, PLL_FALSE, PLL_FALSE);
2171
2172 rearr.rearrangeType = PLL_REARRANGE_SPR;
2173 rearr.likelihood = tr->likelihood;
2174 rearr.SPR.removeNode = p;
2175 rearr.SPR.insertNode = q;
2176 for (i = 0; i < numBranches; ++ i)
2177 {
2178 rearr.SPR.zqr[i] = tr->zqr[i];
2179 }
2180
2181 pllStoreRearrangement (bestList, &rearr);
2182
2183 /*
2184 if(tr->likelihood > tr->bestOfNode)
2185 {
2186 pllStoreRearrangement (bestList, rearr)
2187 tr->bestOfNode = tr->likelihood;
2188 tr->insertNode = q;
2189 tr->removeNode = p;
2190 for(i = 0; i < numBranches; i++)
2191 {
2192 tr->currentZQR[i] = tr->zqr[i];
2193 tr->currentLZR[i] = tr->lzr[i];
2194 tr->currentLZQ[i] = tr->lzq[i];
2195 tr->currentLZS[i] = tr->lzs[i];
2196 }
2197 }
2198
2199 if(tr->likelihood > tr->endLH)
2200 {
2201
2202 tr->insertNode = q;
2203 tr->removeNode = p;
2204 for(i = 0; i < numBranches; i++)
2205 tr->currentZQR[i] = tr->zqr[i];
2206 tr->endLH = tr->likelihood;
2207 }
2208 */
2209 /* reset the topology so that it is the same as it was before calling insertBIG */
2210 hookup(q, r, qz, numBranches);
2211
2212 p->next->next->back = p->next->back = (nodeptr) NULL;
2213
2214 if(tr->thoroughInsertion)
2215 {
2216 nodeptr s = p->back;
2217 hookup(p, s, pz, numBranches);
2218 }
2219
2220 /*
2221 if((tr->doCutoff) && (tr->likelihood < startLH))
2222 {
2223 tr->lhAVG += (startLH - tr->likelihood);
2224 tr->lhDEC++;
2225 if((startLH - tr->likelihood) >= tr->lhCutoff)
2226 return PLL_FALSE;
2227 else
2228 return PLL_TRUE;
2229 }
2230 else
2231 return PLL_TRUE;
2232 */
2233 return (PLL_TRUE);
2234 }
2235
2236 /** @ingroup rearrangementGroup
2237 @brief Internal function for recursively traversing a tree and testing a possible subtree insertion
2238
2239 Recursively traverses the tree rooted at \a q in the direction of \a q->next->back and \a q->next->next->back
2240 and at each node tests the placement of the pruned subtree rooted at \a p by calling the function
2241 \a pllTestInsertBIG, which in turn saves the computed SPR in \a bestList if a) there is still space in
2242 the \a bestList or b) if the likelihood of the SPR is better than any of the ones in \a bestList.
2243
2244 @note This function is not part of the API and should not be called by the user.
2245 */
pllTraverseUpdate(pllInstance * tr,partitionList * pr,nodeptr p,nodeptr q,int mintrav,int maxtrav,pllRearrangeList * bestList)2246 static void pllTraverseUpdate (pllInstance *tr, partitionList *pr, nodeptr p, nodeptr q, int mintrav, int maxtrav, pllRearrangeList * bestList)
2247 {
2248 if (--mintrav <= 0)
2249 {
2250 if (! pllTestInsertBIG(tr, pr, p, q, bestList)) return;
2251
2252 }
2253
2254 if ((!isTip(q->number, tr->mxtips)) && (--maxtrav > 0))
2255 {
2256 pllTraverseUpdate(tr, pr, p, q->next->back, mintrav, maxtrav, bestList);
2257 pllTraverseUpdate(tr, pr, p, q->next->next->back, mintrav, maxtrav, bestList);
2258 }
2259 }
2260
2261
2262 /** @ingroup rearrangementGroup
2263 @brief Internal function for computing SPR moves
2264
2265 Compute a list of at most \a max SPR moves that can be performed by pruning
2266 the subtree rooted at node \a p and testing all possible placements in a
2267 radius of at least \a mintrav nodes and at most \a maxtrav nodes from \a p.
2268 Note that \a tr->thoroughInsertion affects the behaviour of the function (see note).
2269
2270 @param tr
2271 PLL instance
2272
2273 @param pr
2274 List of partitions
2275
2276 @param p
2277 Node specifying the root of the pruned subtree, i.e. where to prune.
2278
2279 @param mintrav
2280 Minimum distance from \a p where to try inserting the pruned subtree
2281
2282 @param maxtrav
2283 Maximum distance from \a p where to try inserting the pruned subtree
2284
2285 @param bestList
2286 The list of best topological rearrangements
2287
2288 @note This function is not part of the API and should not be called by the user
2289 as it is called internally by the API function \a pllComputeSPR.
2290 Also, setting \a tr->thoroughInsertion affects this function. For each tested SPR
2291 the new branch lengths will also be optimized. This computes better likelihoods
2292 but also slows down the method considerably.
2293 */
pllTestSPR(pllInstance * tr,partitionList * pr,nodeptr p,int mintrav,int maxtrav,pllRearrangeList * bestList)2294 static int pllTestSPR (pllInstance * tr, partitionList * pr, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList)
2295 {
2296 nodeptr
2297 p1, p2, q, q1, q2;
2298 double
2299 p1z[PLL_NUM_BRANCHES], p2z[PLL_NUM_BRANCHES], q1z[PLL_NUM_BRANCHES], q2z[PLL_NUM_BRANCHES];
2300 int
2301 mintrav2, i;
2302 int numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
2303
2304 if (maxtrav < 1 || mintrav > maxtrav) return (PLL_FALSE);
2305 q = p->back;
2306
2307 if (!isTip (p->number, tr->mxtips))
2308 {
2309 p1 = p->next->back;
2310 p2 = p->next->next->back;
2311
2312 if (!isTip (p1->number, tr->mxtips) || !isTip (p2->number, tr->mxtips))
2313 {
2314 /* save branch lengths before splitting the tree in two components */
2315 for (i = 0; i < numBranches; ++ i)
2316 {
2317 p1z[i] = p1->z[i];
2318 p2z[i] = p2->z[i];
2319 }
2320
2321 /* split the tree in two components */
2322 if (! removeNodeBIG (tr, pr, p, numBranches)) return PLL_BADREAR;
2323
2324 /* recursively traverse and perform SPR on the subtree rooted at p1 */
2325 if (!isTip (p1->number, tr->mxtips))
2326 {
2327 pllTraverseUpdate (tr, pr, p, p1->next->back, mintrav, maxtrav, bestList);
2328 pllTraverseUpdate (tr, pr, p, p1->next->next->back, mintrav, maxtrav, bestList);
2329 }
2330
2331 /* recursively traverse and perform SPR on the subtree rooted at p2 */
2332 if (!isTip (p2->number, tr->mxtips))
2333 {
2334 pllTraverseUpdate (tr, pr, p, p2->next->back, mintrav, maxtrav, bestList);
2335 pllTraverseUpdate (tr, pr, p, p2->next->next->back, mintrav, maxtrav, bestList);
2336 }
2337
2338 /* restore the topology as it was before the split */
2339 hookup (p->next, p1, p1z, numBranches);
2340 hookup (p->next->next, p2, p2z, numBranches);
2341 pllUpdatePartials (tr, pr, p, PLL_FALSE);
2342 }
2343 }
2344
2345 if (!isTip (q->number, tr->mxtips) && maxtrav > 0)
2346 {
2347 q1 = q->next->back;
2348 q2 = q->next->next->back;
2349
2350 /* why so many conditions? Why is it not analogous to the previous if for node p? */
2351 if (
2352 (
2353 ! isTip(q1->number, tr->mxtips) &&
2354 (! isTip(q1->next->back->number, tr->mxtips) || ! isTip(q1->next->next->back->number, tr->mxtips))
2355 )
2356 ||
2357 (
2358 ! isTip(q2->number, tr->mxtips) &&
2359 (! isTip(q2->next->back->number, tr->mxtips) || ! isTip(q2->next->next->back->number, tr->mxtips))
2360 )
2361 )
2362 {
2363 for (i = 0; i < numBranches; ++ i)
2364 {
2365 q1z[i] = q1->z[i];
2366 q2z[i] = q2->z[i];
2367 }
2368
2369 if (! removeNodeBIG (tr, pr, q, numBranches)) return PLL_BADREAR;
2370
2371 mintrav2 = mintrav > 2 ? mintrav : 2;
2372
2373 if (!isTip (q1->number, tr->mxtips))
2374 {
2375 pllTraverseUpdate (tr, pr, q, q1->next->back, mintrav2, maxtrav, bestList);
2376 pllTraverseUpdate (tr, pr, q, q1->next->next->back, mintrav2, maxtrav, bestList);
2377 }
2378
2379 if (!isTip (q2->number, tr->mxtips))
2380 {
2381 pllTraverseUpdate (tr, pr, q, q2->next->back, mintrav2, maxtrav, bestList);
2382 pllTraverseUpdate (tr, pr, q, q2->next->next->back, mintrav2, maxtrav, bestList);
2383 }
2384
2385 hookup (q->next, q1, q1z, numBranches);
2386 hookup (q->next->next, q2, q2z, numBranches);
2387 pllUpdatePartials (tr, pr, q, PLL_FALSE);
2388 }
2389 }
2390 return (PLL_TRUE);
2391 }
2392
2393 /** @ingroup rearrangementGroup
2394 @brief Compute a list of possible SPR moves
2395
2396 Iteratively tries all possible SPR moves that can be performed by
2397 pruning the subtree rooted at \a p and testing all possible placements
2398 in a radius of at least \a mintrav nodea and at most \a maxtrav nodes from
2399 \a p. Note that \a tr->thoroughInsertion affects the behaviour of the function (see note).
2400
2401 @param tr
2402 PLL instance
2403
2404 @param pr
2405 List of partitions
2406
2407 @param p
2408 Node specifying the root of the pruned subtree, i.e. where to prune.
2409
2410 @param mintrav
2411 Minimum distance from \a p where to try inserting the pruned subtree
2412
2413 @param maxtrav
2414 Maximum distance from \a p where to try inserting the pruned subtree
2415
2416 @note
2417 Setting \a tr->thoroughInsertion affects this function. For each tested SPR
2418 the new branch lengths will also be optimized. This computes better likelihoods
2419 but also slows down the method considerably.
2420 */
2421 static void
pllComputeSPR(pllInstance * tr,partitionList * pr,nodeptr p,int mintrav,int maxtrav,pllRearrangeList * bestList)2422 pllComputeSPR (pllInstance * tr, partitionList * pr, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList)
2423 {
2424
2425 tr->startLH = tr->endLH = tr->likelihood;
2426
2427 /* TODO: Add cutoff code */
2428
2429 tr->bestOfNode = PLL_UNLIKELY;
2430
2431 pllTestSPR (tr, pr, p, mintrav, maxtrav, bestList);
2432 }
2433
2434 /** @ingroup rearrangementGroup
2435 @brief Return the yielded likelihood of an NNI move, without altering the topology
2436
2437 This function performs the NNI move of type \a swapType at node \a p, optimizes
2438 the branch with endpoints \a p and \a p->back and evalutes the resulting likelihood.
2439 It then restores the topology to the origin and returns the likelihood that the NNI
2440 move yielded.
2441
2442 @param tr
2443 PLL instance
2444
2445 @param pr
2446 List of partitions
2447
2448 @param p
2449 Where to perform the NNI move
2450
2451 @param swapType
2452 What type of NNI move to perform
2453
2454 @return
2455 The likelihood yielded from the NNI
2456 */
2457 static double
pllTestNNILikelihood(pllInstance * tr,partitionList * pr,nodeptr p,int swapType)2458 pllTestNNILikelihood (pllInstance * tr, partitionList * pr, nodeptr p, int swapType)
2459 {
2460 double lh;
2461 double z0[PLL_NUM_BRANCHES];
2462 int i;
2463
2464 /* store the origin branch lengths and likelihood. The original branch lengths could
2465 be passed as a parameter in order to avoid duplicate computations because of the two
2466 NNI moves */
2467 for (i = 0; i < pr->numberOfPartitions; ++ i)
2468 {
2469 z0[i] = p->z[i];
2470 }
2471
2472 /* perform NNI */
2473 pllTopologyPerformNNI(tr, p, swapType);
2474 /* recompute the likelihood vectors of the two subtrees rooted at p and p->back,
2475 optimize the branch lengths and evaluate the likelihood */
2476 pllUpdatePartials (tr, pr, p, PLL_FALSE);
2477 pllUpdatePartials (tr, pr, p->back, PLL_FALSE);
2478 update (tr, pr, p);
2479 pllEvaluateLikelihood (tr, pr, p, PLL_FALSE, PLL_FALSE);
2480 lh = tr->likelihood;
2481
2482 /* restore topology */
2483 pllTopologyPerformNNI(tr, p, swapType);
2484 pllUpdatePartials (tr, pr, p, PLL_FALSE);
2485 pllUpdatePartials (tr, pr, p->back, PLL_FALSE);
2486 //update (tr, pr, p);
2487 pllEvaluateLikelihood (tr, pr, p, PLL_FALSE, PLL_FALSE);
2488 for (i = 0; i < pr->numberOfPartitions; ++ i)
2489 {
2490 p->z[i] = p->back->z[i] = z0[i];
2491 }
2492
2493 return lh;
2494 }
2495 /** @ingroup rearrangementGroup
2496 @brief Compares NNI likelihoods at a node and store in the rearrangement list
2497
2498 Compares the two possible NNI moves that can be performed at node \a p, and
2499 if the likelihood improves from the one of the original topology, then
2500 it picks the one that yields the highest likelihood and tries to insert it in
2501 the list of best rearrangement moves
2502
2503 @param tr
2504 PLL instance
2505
2506 @param pr
2507 List of partitions
2508
2509 @param bestList
2510 Rearrangement moves list
2511 */
pllTestNNI(pllInstance * tr,partitionList * pr,nodeptr p,pllRearrangeList * bestList)2512 static void pllTestNNI (pllInstance * tr, partitionList * pr, nodeptr p, pllRearrangeList * bestList)
2513 {
2514 double lh0, lh1, lh2;
2515 pllRearrangeInfo rearr;
2516
2517 /* store the original likelihood */
2518 lh0 = tr->likelihood;
2519
2520 lh1 = pllTestNNILikelihood (tr, pr, p, PLL_NNI_P_NEXT);
2521 lh2 = pllTestNNILikelihood (tr, pr, p, PLL_NNI_P_NEXTNEXT);
2522
2523 if (lh0 > lh1 && lh0 > lh2) return;
2524
2525 /* set the arrangement structure */
2526 rearr.rearrangeType = PLL_REARRANGE_NNI;
2527 rearr.likelihood = PLL_MAX (lh1, lh2);
2528 rearr.NNI.originNode = p;
2529 rearr.NNI.swapType = (lh1 > lh2) ? PLL_NNI_P_NEXT : PLL_NNI_P_NEXTNEXT;
2530
2531 /* try to store it in the best list */
2532 pllStoreRearrangement (bestList, &rearr);
2533 }
2534
2535 /** @ingroup rearrangementGroup
2536 @brief Recursive traversal of the tree structure for testing NNI moves
2537
2538 Recursively traverses the tree structure and tests all allowed NNI
2539 moves in the area specified by \a mintrav and \a maxtrav. For more
2540 information and details on the function arguments check ::pllSearchNNI
2541 */
2542 static void
pllTraverseNNI(pllInstance * tr,partitionList * pr,nodeptr p,int mintrav,int maxtrav,pllRearrangeList * bestList)2543 pllTraverseNNI (pllInstance *tr, partitionList *pr, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList)
2544 {
2545 if (isTip (p->number, tr->mxtips)) return;
2546
2547 /* if we are at the right radius then compute the NNIs for nodes p->next and p->next->next */
2548 if (!mintrav)
2549 {
2550 pllTestNNI (tr, pr, p->next, bestList);
2551 pllTestNNI (tr, pr, p->next->next, bestList);
2552 }
2553
2554 /* and then avoid computing the NNIs for nodes p->next->back and p->next->next->back as they are
2555 the same to the ones computed in the above two lines. This way we do not need to resolve conflicts
2556 later on as in the old code */
2557 if (maxtrav)
2558 {
2559 if (!isTip (p->next->back->number, tr->mxtips))
2560 pllTraverseNNI (tr, pr, p->next->back, mintrav ? mintrav - 1 : 0, maxtrav - 1, bestList);
2561 if (!isTip (p->next->next->back->number, tr->mxtips))
2562 pllTraverseNNI (tr, pr, p->next->next->back, mintrav ? mintrav - 1 : 0, maxtrav - 1, bestList);
2563 }
2564 }
2565
2566 /** @ingroup rearrangementGroup
2567 @brief Compute a list of possible NNI moves
2568
2569 Iteratively tries all possible NNI moves at each node that is at
2570 least \a mintrav and at most \a maxtrav nodes far from node \a p.
2571 At each NNI move, the likelihood is tested and if it is higher than
2572 the likelihood of an element in the sorted (by likelihood) list
2573 \a bestList, or if there is still empty space in \a bestList, it is
2574 inserted at the corresponding position.
2575
2576 @param tr
2577 PLL instance
2578
2579 @param pr
2580 List of partitions
2581
2582 @param p
2583 Node specifying the point where the NNI will be performed.
2584
2585 @param mintrav
2586 Minimum distance from \a p where the NNI can be tested
2587
2588 @param maxtrav
2589 Maximum distance from \a p where to try NNIs
2590 */
2591 static void
pllSearchNNI(pllInstance * tr,partitionList * pr,nodeptr p,int mintrav,int maxtrav,pllRearrangeList * bestList)2592 pllSearchNNI (pllInstance * tr, partitionList * pr, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList)
2593 {
2594 /* avoid conflicts by precomputing the NNI of the first node */
2595
2596 if (mintrav == 0)
2597 pllTestNNI (tr, pr, p, bestList);
2598
2599 pllTraverseNNI (tr, pr, p, mintrav, maxtrav, bestList);
2600 if (maxtrav)
2601 pllTraverseNNI (tr, pr, p->back, mintrav, maxtrav - 1, bestList);
2602
2603 }
2604
2605 /** @ingroup rearrangementGroup
2606 @brief Create rollback information for an SPR move
2607
2608 Creates a structure of type ::pllRollbackInfo and fills it with rollback
2609 information about the SPR move described in \a rearr. The rollback info
2610 is stored in the PLL instance in a LIFO manner.
2611
2612 @param tr
2613 PLL instance
2614
2615 @param rearr
2616 Description of the SPR move
2617
2618 @param numBranches
2619 Number of partitions
2620 */
2621 static void
pllCreateSprInfoRollback(pllInstance * tr,pllRearrangeInfo * rearr,int numBranches)2622 pllCreateSprInfoRollback (pllInstance * tr, pllRearrangeInfo * rearr, int numBranches)
2623 {
2624 pllRollbackInfo * sprRb;
2625 nodeptr p, q;
2626 int i;
2627
2628 p = rearr->SPR.removeNode;
2629 q = rearr->SPR.insertNode;
2630
2631 sprRb = (pllRollbackInfo *) rax_malloc (sizeof (pllRollbackInfo) + 4 * numBranches * sizeof (double));
2632 sprRb->SPR.zp = (double *) ((char *)sprRb + sizeof (pllRollbackInfo));
2633 sprRb->SPR.zpn = (double *) ((char *)sprRb + sizeof (pllRollbackInfo) + numBranches * sizeof (double));
2634 sprRb->SPR.zpnn = (double *) ((char *)sprRb + sizeof (pllRollbackInfo) + 2 * numBranches * sizeof (double));
2635 sprRb->SPR.zqr = (double *) ((char *)sprRb + sizeof (pllRollbackInfo) + 3 * numBranches * sizeof (double));
2636
2637 for (i = 0; i < numBranches; ++ i)
2638 {
2639 sprRb->SPR.zp[i] = p->z[i];
2640 sprRb->SPR.zpn[i] = p->next->z[i];
2641 sprRb->SPR.zpnn[i] = p->next->next->z[i];
2642 sprRb->SPR.zqr[i] = q->z[i];
2643 }
2644
2645 sprRb->SPR.pn = p->next->back;
2646 sprRb->SPR.pnn = p->next->next->back;
2647 sprRb->SPR.r = q->back;
2648 sprRb->SPR.q = q;
2649 sprRb->SPR.p = p;
2650
2651 sprRb->rearrangeType = PLL_REARRANGE_SPR;
2652
2653 pllStackPush (&(tr->rearrangeHistory), (void *) sprRb);
2654 }
2655
2656 /** @ingroup rearrangementGroup
2657 @brief Create rollback information for an NNI move
2658
2659 Creates a structure of type ::pllRollbackInfo and fills it with rollback
2660 information about the SPR move described in \a rearr. The rollback info
2661 is stored in the PLL instance in a LIFO manner
2662
2663 @param tr
2664 PLL instance
2665
2666 @param rearr
2667 Description of the NNI move
2668 */
2669 static void
pllCreateNniInfoRollback(pllInstance * tr,pllRearrangeInfo * rearr)2670 pllCreateNniInfoRollback (pllInstance * tr, pllRearrangeInfo * rearr)
2671 {
2672 /*TODO: add the branches ? */
2673 pllRollbackInfo * ri;
2674
2675 ri = (pllRollbackInfo *) rax_malloc (sizeof (pllRollbackInfo));
2676
2677 ri->rearrangeType = PLL_REARRANGE_NNI;
2678
2679 ri->NNI.origin = rearr->NNI.originNode;
2680 ri->NNI.swapType = rearr->NNI.swapType;
2681
2682 pllStackPush (&(tr->rearrangeHistory), (void *) ri);
2683
2684 }
2685
2686
2687 /** @ingroup rearrangementGroup
2688 @brief Generic function for creating rollback information
2689
2690 Creates a structure of type ::pllRollbackInfo and fills it with rollback
2691 information about the move described in \a rearr. The rollback info
2692 is stored in the PLL instance in a LIFO manner
2693
2694 @param tr
2695 PLL instance
2696
2697 @param rearr
2698 Description of the NNI move
2699
2700 @param numBranches
2701 Number of partitions
2702 */
2703 static void
pllCreateRollbackInfo(pllInstance * tr,pllRearrangeInfo * rearr,int numBranches)2704 pllCreateRollbackInfo (pllInstance * tr, pllRearrangeInfo * rearr, int numBranches)
2705 {
2706 switch (rearr->rearrangeType)
2707 {
2708 case PLL_REARRANGE_NNI:
2709 pllCreateNniInfoRollback (tr, rearr);
2710 break;
2711 case PLL_REARRANGE_SPR:
2712 pllCreateSprInfoRollback (tr, rearr, numBranches);
2713 break;
2714 default:
2715 break;
2716 }
2717
2718 }
2719
2720
2721 /** @ingroup rearrangementGroup
2722 @brief Rollback an SPR move
2723
2724 Perform a rollback (undo) on the last SPR move.
2725
2726 @param tr
2727 PLL instance
2728
2729 @param pr
2730 List of partitions
2731
2732 @param ri
2733 Rollback information
2734 */
2735 static void
pllRollbackSPR(partitionList * pr,pllRollbackInfo * ri)2736 pllRollbackSPR (partitionList * pr, pllRollbackInfo * ri)
2737 {
2738 int numBranches;
2739
2740 numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
2741
2742 hookup (ri->SPR.p->next, ri->SPR.pn, ri->SPR.zpn, numBranches);
2743 hookup (ri->SPR.p->next->next, ri->SPR.pnn, ri->SPR.zpnn, numBranches);
2744 hookup (ri->SPR.p, ri->SPR.p->back, ri->SPR.zp, numBranches);
2745 hookup (ri->SPR.q, ri->SPR.r, ri->SPR.zqr, numBranches);
2746
2747 rax_free (ri);
2748 }
2749
2750 /** @ingroup rearrangementGroup
2751 @brief Rollback an NNI move
2752
2753 Perform a rollback (undo) on the last NNI move.
2754
2755 @param tr
2756 PLL instance
2757
2758 @param pr
2759 List of partitions
2760
2761 @param ri
2762 Rollback information
2763 */
2764 static void
pllRollbackNNI(pllInstance * tr,partitionList * pr,pllRollbackInfo * ri)2765 pllRollbackNNI (pllInstance * tr, partitionList * pr, pllRollbackInfo * ri)
2766 {
2767 nodeptr p = ri->NNI.origin;
2768
2769 pllTopologyPerformNNI(tr, p, ri->NNI.swapType);
2770 pllUpdatePartials (tr, pr, p, PLL_FALSE);
2771 pllUpdatePartials (tr, pr, p->back, PLL_FALSE);
2772 update (tr, pr, p);
2773 pllEvaluateLikelihood (tr, pr, p, PLL_FALSE, PLL_FALSE);
2774
2775 rax_free (ri);
2776 }
2777
2778 /** @ingroup rearrangementGroup
2779 @brief Rollback the last committed rearrangement move
2780
2781 Perform a rollback (undo) on the last committed rearrangement move.
2782
2783 @param tr
2784 PLL instance
2785
2786 @param pr
2787 List of partitions
2788
2789 @return
2790 Returns \b PLL_TRUE is the rollback was successful, otherwise \b PLL_FALSE
2791 (if no rollback was done)
2792 */
2793 int
pllRearrangeRollback(pllInstance * tr,partitionList * pr)2794 pllRearrangeRollback (pllInstance * tr, partitionList * pr)
2795 {
2796 pllRollbackInfo * ri;
2797
2798 ri = (pllRollbackInfo *) pllStackPop (&(tr->rearrangeHistory));
2799 if (!ri) return (PLL_FALSE);
2800
2801 switch (ri->rearrangeType)
2802 {
2803 case PLL_REARRANGE_NNI:
2804 pllRollbackNNI (tr, pr, ri);
2805 break;
2806 case PLL_REARRANGE_SPR:
2807 pllRollbackSPR (pr, ri);
2808 break;
2809 default:
2810 rax_free (ri);
2811 return (PLL_FALSE);
2812 }
2813
2814 return (PLL_TRUE);
2815
2816 }
2817
2818
2819 /** @ingroup rearrangementGroup
2820 @brief Commit a rearrangement move
2821
2822 Applies the rearrangement move specified in \a rearr to the tree topology in \a tr.
2823 In case of SPR moves, if
2824 \a tr->thoroughInsertion is set to \b PLL_TRUE, the new branch lengths are also optimized.
2825 The function stores rollback information in pllInstance::rearrangeHistory if \a saveRollbackInfo
2826 is set to \b PLL_TRUE. This way, the rearrangement move can be rolled back (undone) by calling
2827 ::pllRearrangeRollback
2828
2829 @param tr
2830 PLL instance
2831
2832 @param pr
2833 List of partitions
2834
2835 @param rearr
2836 An element of a \a pllRearrangeInfo structure that contains information about the rearrangement move
2837
2838 @param saveRollbackInfo
2839 If set to \b PLL_TRUE, rollback info will be kept for undoing the rearrangement move
2840 */
2841 void
pllRearrangeCommit(pllInstance * tr,partitionList * pr,pllRearrangeInfo * rearr,int saveRollbackInfo)2842 pllRearrangeCommit (pllInstance * tr, partitionList * pr, pllRearrangeInfo * rearr, int saveRollbackInfo)
2843 {
2844 int numBranches;
2845
2846 numBranches = pr->perGeneBranchLengths ? pr->numberOfPartitions : 1;
2847
2848 if (saveRollbackInfo)
2849 pllCreateRollbackInfo (tr, rearr, numBranches);
2850
2851 switch (rearr->rearrangeType)
2852 {
2853 case PLL_REARRANGE_NNI:
2854 pllTopologyPerformNNI(tr, rearr->NNI.originNode, rearr->NNI.swapType);
2855 pllUpdatePartials (tr, pr, rearr->NNI.originNode, PLL_FALSE);
2856 pllUpdatePartials (tr, pr, rearr->NNI.originNode->back, PLL_FALSE);
2857 update (tr, pr, rearr->NNI.originNode);
2858 pllEvaluateLikelihood (tr, pr, rearr->NNI.originNode, PLL_FALSE, PLL_FALSE);
2859 break;
2860 case PLL_REARRANGE_SPR:
2861 removeNodeBIG (tr, pr, rearr->SPR.removeNode, numBranches);
2862 insertBIG (tr, pr, rearr->SPR.removeNode, rearr->SPR.insertNode);
2863 break;
2864 default:
2865 break;
2866 }
2867 }
2868
2869
2870 /******** new rearrangement functions ****************/
2871
2872 /* change this to return the number of new elements in the list */
2873 /** @ingroup rearrangementGroup
2874 @brief Search for rearrangement topologies
2875
2876 Search for possible rearrangement moves of type \a rearrangeType in the
2877 annular area defined by the minimal resp. maximal radii \a mintrav resp.
2878 \a maxtrav. If the resulting likelihood is better than the current, try
2879 to insert the move specification in \a bestList, which is a sorted list
2880 that holds the rearrange info of the best moves sorted by likelihood
2881 (desccending order).
2882
2883 @param tr
2884 PLL instance
2885
2886 @param pr
2887 List of partitions
2888
2889 @param rearrangeType
2890 Type of rearrangement. Can be \b PLL_REARRANGE_SPR or \b PLL_REARRANGE_NNI
2891
2892 @param p
2893 Point of origin, i.e. where to start searching from
2894
2895 @param mintrav
2896 The minimal radius of the annulus
2897
2898 @param maxtrav
2899 The maximal radius of the annulus
2900
2901 @param bestList
2902 List that holds the details of the best rearrangement moves found
2903
2904 @note
2905 If \a bestList is not empty, the existing entries will not be altered unless
2906 better rearrangement moves (that means yielding better likelihood) are found
2907 and the list is full, in which case the entries with the worst likelihood will be
2908 thrown away.
2909 */
2910 void
pllRearrangeSearch(pllInstance * tr,partitionList * pr,int rearrangeType,nodeptr p,int mintrav,int maxtrav,pllRearrangeList * bestList)2911 pllRearrangeSearch (pllInstance * tr, partitionList * pr, int rearrangeType, nodeptr p, int mintrav, int maxtrav, pllRearrangeList * bestList)
2912 {
2913 switch (rearrangeType)
2914 {
2915 case PLL_REARRANGE_SPR:
2916 pllComputeSPR (tr, pr, p, mintrav, maxtrav, bestList);
2917 break;
2918
2919 case PLL_REARRANGE_NNI:
2920 pllSearchNNI (tr, pr, p, mintrav, maxtrav, bestList);
2921 break;
2922
2923 case PLL_REARRANGE_TBR:
2924 break;
2925 default:
2926 break;
2927 }
2928 }
2929
2930
2931 static int
determineRearrangementSetting(pllInstance * tr,partitionList * pr,bestlist * bestT,bestlist * bt)2932 determineRearrangementSetting(pllInstance *tr, partitionList *pr,
2933 bestlist *bestT, bestlist *bt)
2934 {
2935 int i, mintrav, maxtrav, bestTrav, impr, index, MaxFast, *perm = (int*) NULL;
2936 double startLH;
2937 pllBoolean cutoff;
2938
2939 MaxFast = 26;
2940
2941 startLH = tr->likelihood;
2942
2943 cutoff = tr->doCutoff;
2944 tr->doCutoff = PLL_FALSE;
2945
2946 mintrav = 1;
2947 maxtrav = 5;
2948
2949 bestTrav = maxtrav = 5;
2950
2951 impr = 1;
2952
2953 resetBestTree(bt);
2954
2955 if (tr->permuteTreeoptimize)
2956 {
2957 int n = tr->mxtips + tr->mxtips - 2;
2958 perm = (int *) rax_malloc(sizeof(int) * (n + 1));
2959 makePermutation(perm, n, tr);
2960 }
2961
2962 while (impr && maxtrav < MaxFast)
2963 {
2964 recallBestTree(bestT, 1, tr, pr);
2965 nodeRectifier(tr);
2966
2967 if (maxtrav > tr->ntips - 3)
2968 maxtrav = tr->ntips - 3;
2969
2970 tr->startLH = tr->endLH = tr->likelihood;
2971
2972 for (i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
2973 {
2974
2975 if (tr->permuteTreeoptimize)
2976 index = perm[i];
2977 else
2978 index = i;
2979
2980 tr->bestOfNode = PLL_UNLIKELY;
2981 if (rearrangeBIG(tr, pr, tr->nodep[index], mintrav, maxtrav))
2982 {
2983 if (tr->endLH > tr->startLH)
2984 {
2985 restoreTreeFast(tr, pr);
2986 tr->startLH = tr->endLH = tr->likelihood;
2987 }
2988 }
2989 }
2990
2991 pllOptimizeBranchLengths(tr, pr, 8);
2992 saveBestTree(bt, tr,
2993 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
2994
2995 if (tr->likelihood > startLH)
2996 {
2997 startLH = tr->likelihood;
2998 bestTrav = maxtrav;
2999 impr = 1;
3000 }
3001 else
3002 {
3003 impr = 0;
3004 }
3005 maxtrav += 5;
3006
3007 if (tr->doCutoff)
3008 {
3009 tr->lhCutoff = (tr->lhAVG) / ((double) (tr->lhDEC));
3010
3011 tr->itCount = tr->itCount + 1;
3012 tr->lhAVG = 0;
3013 tr->lhDEC = 0;
3014 }
3015 }
3016
3017 recallBestTree(bt, 1, tr, pr);
3018 tr->doCutoff = cutoff;
3019
3020 if (tr->permuteTreeoptimize)
3021 rax_free(perm);
3022
3023 return bestTrav;
3024 }
3025
3026
hash_dealloc_bipentry(void * entry)3027 static void hash_dealloc_bipentry (void * entry)
3028 {
3029 pllBipartitionEntry * e = (pllBipartitionEntry *)entry;
3030
3031 if(e->bitVector) rax_free(e->bitVector);
3032 if(e->treeVector) rax_free(e->treeVector);
3033 if(e->supportVector) rax_free(e->supportVector);
3034
3035 }
3036
3037 /** @ingroup rearrangementGroup
3038 @brief RAxML algorithm for ML search
3039
3040 RAxML algorithm for searching the Maximum Likelihood tree and model.
3041
3042 @param tr
3043 PLL instance
3044
3045 @param pr
3046 List of partitions
3047
3048 @param estimateModel
3049 If true, model parameters are optimized in a ML framework.
3050
3051 @note
3052 For datasets with a large number of taxa, setting tr->searchConvergenceCriterion to
3053 PLL_TRUE can improve the execution time in up to 50% looking for topology convergence.
3054 */
3055 int
pllRaxmlSearchAlgorithm(pllInstance * tr,partitionList * pr,pllBoolean estimateModel)3056 pllRaxmlSearchAlgorithm(pllInstance * tr, partitionList * pr,
3057 pllBoolean estimateModel)
3058 {
3059 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
3060 pllOptimizeBranchLengths(tr, pr, 32);
3061
3062 unsigned int vLength = 0;
3063 int i, impr, bestTrav, rearrangementsMax = 0, rearrangementsMin = 0,
3064 thoroughIterations = 0, fastIterations = 0;
3065
3066 double lh, previousLh, difference, epsilon;
3067 bestlist *bestT, *bt;
3068 infoList iList;
3069 pllOptimizeBranchLengths(tr, pr, 32);
3070
3071 pllHashTable *h = NULL;
3072 //hashtable *h = NULL;
3073 unsigned int **bitVectors = (unsigned int**) NULL;
3074
3075 /* Security check... These variables might have not been initialized! */
3076 if (tr->stepwidth == 0) tr->stepwidth = 5;
3077 if (tr->max_rearrange == 0) tr->max_rearrange = 21;
3078
3079 if (tr->searchConvergenceCriterion)
3080 {
3081 bitVectors = initBitVector(tr->mxtips, &vLength);
3082 //h = initHashTable(tr->mxtips * 4);
3083 h = pllHashInit (tr->mxtips * 4);
3084 }
3085
3086 bestT = (bestlist *) rax_malloc(sizeof(bestlist));
3087 bestT->ninit = 0;
3088 initBestTree(bestT, 1, tr->mxtips);
3089
3090 bt = (bestlist *) rax_malloc(sizeof(bestlist));
3091 bt->ninit = 0;
3092 initBestTree(bt, 20, tr->mxtips);
3093
3094 initInfoList(&iList, 50);
3095
3096 difference = 10.0;
3097 epsilon = tr->likelihoodEpsilon;
3098
3099 tr->thoroughInsertion = 0;
3100
3101 if (estimateModel)
3102 {
3103 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
3104 pllOptimizeModelParameters(tr, pr, 10.0);
3105 }
3106 else
3107 pllOptimizeBranchLengths(tr, pr, 64);
3108
3109 saveBestTree(bestT, tr,
3110 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3111
3112 if (!tr->initialSet)
3113 bestTrav = tr->bestTrav = determineRearrangementSetting(tr, pr, bestT, bt);
3114 else
3115 bestTrav = tr->bestTrav = tr->initial;
3116
3117 if (estimateModel)
3118 {
3119 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
3120 pllOptimizeModelParameters(tr, pr, 5.0);
3121 }
3122 else
3123 pllOptimizeBranchLengths(tr, pr, 32);
3124
3125 saveBestTree(bestT, tr,
3126 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3127 impr = 1;
3128 if (tr->doCutoff)
3129 tr->itCount = 0;
3130
3131 while (impr)
3132 {
3133 recallBestTree(bestT, 1, tr, pr);
3134
3135 if (tr->searchConvergenceCriterion)
3136 {
3137 int bCounter = 0;
3138
3139 if (fastIterations > 1)
3140 cleanupHashTable(h, (fastIterations % 2));
3141
3142 bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back, tr->mxtips,
3143 vLength, h, fastIterations % 2, PLL_BIPARTITIONS_RF,
3144 (branchInfo *) NULL, &bCounter, 1, PLL_FALSE, PLL_FALSE, 0);
3145
3146 assert(bCounter == tr->mxtips - 3);
3147
3148 if (fastIterations > 0)
3149 {
3150 double rrf = convergenceCriterion(h, tr->mxtips);
3151
3152 if (rrf <= 0.01) /* 1% cutoff */
3153 {
3154 cleanupHashTable(h, 0);
3155 cleanupHashTable(h, 1);
3156 goto cleanup_fast;
3157 }
3158 }
3159 }
3160
3161 fastIterations++;
3162
3163 pllOptimizeBranchLengths(tr, pr, 32);
3164
3165 saveBestTree(bestT, tr,
3166 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3167
3168 lh = previousLh = tr->likelihood;
3169
3170 treeOptimizeRapid(tr, pr, 1, bestTrav, bt, &iList);
3171
3172 impr = 0;
3173
3174 for (i = 1; i <= bt->nvalid; i++)
3175 {
3176 recallBestTree(bt, i, tr, pr);
3177
3178 pllOptimizeBranchLengths(tr, pr, 8);
3179
3180 difference = (
3181 (tr->likelihood > previousLh) ?
3182 tr->likelihood - previousLh : previousLh - tr->likelihood);
3183 if (tr->likelihood > lh && difference > epsilon)
3184 {
3185 impr = 1;
3186 lh = tr->likelihood;
3187 saveBestTree(bestT, tr,
3188 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3189 }
3190 }
3191 }
3192
3193 if (tr->searchConvergenceCriterion)
3194 {
3195 cleanupHashTable(h, 0);
3196 cleanupHashTable(h, 1);
3197 }
3198
3199 cleanup_fast:
3200
3201 tr->thoroughInsertion = 1;
3202 impr = 1;
3203
3204 recallBestTree(bestT, 1, tr, pr);
3205 if (estimateModel)
3206 {
3207 pllEvaluateLikelihood(tr, pr, tr->start, PLL_TRUE, PLL_FALSE);
3208 pllOptimizeModelParameters(tr, pr, 1.0);
3209 }
3210 else
3211 pllOptimizeBranchLengths(tr, pr, 32);
3212
3213 while (1)
3214 {
3215 recallBestTree(bestT, 1, tr, pr);
3216 if (impr)
3217 {
3218 rearrangementsMin = 1;
3219 rearrangementsMax = tr->stepwidth;
3220
3221 if (tr->searchConvergenceCriterion)
3222 {
3223 int bCounter = 0;
3224
3225 if (thoroughIterations > 1)
3226 cleanupHashTable(h, (thoroughIterations % 2));
3227
3228 bitVectorInitravSpecial(bitVectors, tr->nodep[1]->back,
3229 tr->mxtips, vLength, h, thoroughIterations % 2,
3230 PLL_BIPARTITIONS_RF, (branchInfo *) NULL, &bCounter, 1,
3231 PLL_FALSE, PLL_FALSE, 0);
3232
3233 assert(bCounter == tr->mxtips - 3);
3234
3235 if (thoroughIterations > 0)
3236 {
3237 double rrf = convergenceCriterion(h, tr->mxtips);
3238
3239 if (rrf <= 0.01) /* 1% cutoff */
3240 {
3241 goto cleanup;
3242 }
3243 }
3244 }
3245
3246 thoroughIterations++;
3247 }
3248 else
3249 {
3250 rearrangementsMax += tr->stepwidth;
3251 rearrangementsMin += tr->stepwidth;
3252 if (rearrangementsMax > tr->max_rearrange)
3253 goto cleanup;
3254 }
3255 pllOptimizeBranchLengths(tr, pr, 32);
3256
3257 previousLh = lh = tr->likelihood;
3258 saveBestTree(bestT, tr,
3259 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3260
3261 treeOptimizeRapid(tr, pr, rearrangementsMin, rearrangementsMax, bt,
3262 &iList);
3263
3264 impr = 0;
3265
3266 for (i = 1; i <= bt->nvalid; i++)
3267 {
3268 recallBestTree(bt, i, tr, pr);
3269
3270 pllOptimizeBranchLengths(tr, pr, 8);
3271
3272 difference = (
3273 (tr->likelihood > previousLh) ?
3274 tr->likelihood - previousLh : previousLh - tr->likelihood);
3275 if (tr->likelihood > lh && difference > epsilon)
3276 {
3277 impr = 1;
3278 lh = tr->likelihood;
3279 saveBestTree(bestT, tr,
3280 pr->perGeneBranchLengths ? pr->numberOfPartitions : 1);
3281 }
3282 }
3283
3284 }
3285
3286 cleanup:
3287 if (tr->searchConvergenceCriterion)
3288 {
3289 freeBitVectors(bitVectors, 2 * tr->mxtips);
3290 rax_free(bitVectors);
3291 //freeHashTable(h);
3292 //rax_free(h);
3293 pllHashDestroy(&h, hash_dealloc_bipentry);
3294 }
3295
3296 freeBestTree(bestT);
3297 rax_free(bestT);
3298 freeBestTree(bt);
3299 rax_free(bt);
3300
3301 freeInfoList(&iList);
3302
3303 if (estimateModel) {
3304 pllOptimizeModelParameters(tr, pr, epsilon);
3305 }
3306 pllOptimizeBranchLengths(tr, pr, 64);
3307
3308 return 0;
3309 }
3310
3311