source: trunk/GDE/RAxML/topologies.c

Last change on this file was 10409, checked in by aboeckma, 11 years ago

Updated raxml to current version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1
2/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
3 *  Copyright August 2006 by Alexandros Stamatakis
4 *
5 *  Partially derived from
6 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
7 * 
8 *  and
9 *
10 *  Programs of the PHYLIP package by Joe Felsenstein.
11 *
12 *  This program is free software; you may redistribute it and/or modify its
13 *  under the terms of the GNU General Public License as published by the Free
14 *  Software Foundation; either version 2 of the License, or (at your option)
15 *  any later version.
16 *
17 *  This program is distributed in the hope that it will be useful, but
18 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
19 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
20 *  for more details.
21 *
22 *
23 *  For any other enquiries send an Email to Alexandros Stamatakis
24 *  Alexandros.Stamatakis@epfl.ch
25 *
26 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
27 *
28 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
29 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30 */
31
32#ifndef WIN32
33#include <sys/times.h>
34#include <sys/types.h>
35#include <sys/time.h>
36#include <unistd.h>
37#endif
38
39#include <math.h>
40#include <time.h>
41#include <stdlib.h>
42#include <stdio.h>
43#include <ctype.h>
44#include <string.h>
45
46#include "axml.h"
47
48
49
50
51
52static void saveTopolRELLRec(tree *tr, nodeptr p, topolRELL *tpl, int *i, int numsp, int numBranches)
53{
54  int k;
55  if(isTip(p->number, numsp))
56    return;
57  else
58    {
59      nodeptr q = p->next;     
60      while(q != p)
61        {         
62          tpl->connect[*i].p = q;
63          tpl->connect[*i].q = q->back; 
64         
65          if(tr->grouped ||  tr->constrained)
66            {
67              tpl->connect[*i].cp = tr->constraintVector[q->number];
68              tpl->connect[*i].cq = tr->constraintVector[q->back->number]; 
69            }
70         
71          for(k = 0; k < numBranches; k++)
72            tpl->connect[*i].z[k] = q->z[k];
73          *i = *i + 1;
74         
75          saveTopolRELLRec(tr, q->back, tpl, i, numsp, numBranches);
76          q = q->next;
77        }
78    }
79}
80
81static void saveTopolRELL(tree *tr, topolRELL *tpl)
82{
83  nodeptr p = tr->start;
84  int k, i = 0;
85     
86  tpl->likelihood = tr->likelihood;
87  tpl->start      = 1;
88     
89  tpl->connect[i].p = p;
90  tpl->connect[i].q = p->back;
91 
92  if(tr->grouped ||  tr->constrained)
93    {
94      tpl->connect[i].cp = tr->constraintVector[p->number];
95      tpl->connect[i].cq = tr->constraintVector[p->back->number]; 
96    }
97
98  for(k = 0; k < tr->numBranches; k++)
99    tpl->connect[i].z[k] = p->z[k];
100  i++;
101     
102  saveTopolRELLRec(tr, p->back, tpl, &i, tr->rdta->numsp, tr->numBranches);   
103
104  assert(i == 2 * tr->ntips - 3);
105}
106
107
108static void restoreTopolRELL(tree *tr, topolRELL *tpl)
109{
110  int i;
111 
112  for (i = 0; i < 2 * tr->mxtips - 3; i++) 
113    {
114      hookup(tpl->connect[i].p, tpl->connect[i].q, tpl->connect[i].z,  tr->numBranches);   
115      tr->constraintVector[tpl->connect[i].p->number] = tpl->connect[i].cp;
116      tr->constraintVector[tpl->connect[i].q->number] = tpl->connect[i].cq;
117    }
118 
119
120  tr->likelihood = tpl->likelihood;
121  tr->start      = tr->nodep[tpl->start];
122  /* TODO */
123}
124
125
126
127
128void initTL(topolRELL_LIST *rl, tree *tr, int n)
129{
130  int i;
131
132  rl->max = n; 
133  rl->t = (topolRELL **)rax_malloc(sizeof(topolRELL *) * n);
134
135  for(i = 0; i < n; i++)
136    {
137      rl->t[i] = (topolRELL *)rax_malloc(sizeof(topolRELL));
138      rl->t[i]->connect = (connectRELL *)rax_malloc((2 * tr->mxtips - 3) * sizeof(connectRELL));
139      rl->t[i]->likelihood = unlikely;     
140    }
141}
142
143
144void freeTL(topolRELL_LIST *rl)
145{
146  int i;
147  for(i = 0; i < rl->max; i++)   
148    {
149      rax_free(rl->t[i]->connect);         
150      rax_free(rl->t[i]);
151    }
152  rax_free(rl->t);
153}
154
155
156void restoreTL(topolRELL_LIST *rl, tree *tr, int n)
157{
158  assert(n >= 0 && n < rl->max);   
159
160  restoreTopolRELL(tr, rl->t[n]); 
161}
162
163
164
165
166void resetTL(topolRELL_LIST *rl)
167{
168  int i;
169
170  for(i = 0; i < rl->max; i++)   
171    rl->t[i]->likelihood = unlikely;         
172}
173
174
175
176
177void saveTL(topolRELL_LIST *rl, tree *tr, int index)
178{ 
179  assert(index >= 0 && index < rl->max);   
180   
181  if(tr->likelihood > rl->t[index]->likelihood)       
182    saveTopolRELL(tr, rl->t[index]); 
183}
184
185
186static void  *tipValPtr (nodeptr p)
187{ 
188  return  (void *) & p->number;
189}
190
191
192static int  cmpTipVal (void *v1, void *v2)
193{
194  int  i1, i2;
195 
196  i1 = *((int *) v1);
197  i2 = *((int *) v2);
198  return  (i1 < i2) ? -1 : ((i1 == i2) ? 0 : 1);
199}
200
201
202/*  These are the only routines that need to UNDERSTAND topologies */
203
204static topol  *setupTopol (int maxtips)
205{
206  topol   *tpl;
207
208  if (! (tpl = (topol *) rax_malloc(sizeof(topol))) || 
209      ! (tpl->links = (connptr) rax_malloc((2*maxtips-3) * sizeof(connect))))
210    {
211      printf("ERROR: Unable to get topology memory");
212      tpl = (topol *) NULL;
213    }
214  else 
215    {
216      tpl->likelihood  = unlikely;
217      tpl->start       = (node *) NULL;
218      tpl->nextlink    = 0;
219      tpl->ntips       = 0;
220      tpl->nextnode    = 0;   
221      tpl->scrNum      = 0;     /* position in sorted list of scores */
222      tpl->tplNum      = 0;     /* position in sorted list of trees */       
223    }
224 
225  return  tpl;
226} 
227
228
229static void  freeTopol (topol *tpl)
230{
231  rax_free(tpl->links);
232  rax_free(tpl);
233} 
234
235
236static int saveSubtree (nodeptr p, topol *tpl, int numsp, int numBranches) 
237{
238  connptr  r, r0;
239  nodeptr  q, s;
240  int      t, t0, t1, k;
241
242  r0 = tpl->links;
243  r = r0 + (tpl->nextlink)++;
244  r->p = p;
245  r->q = q = p->back;
246
247  for(k = 0; k < numBranches; k++)
248    r->z[k] = p->z[k];
249
250  r->descend = 0;                     /* No children (yet) */
251
252  if (isTip(q->number, numsp)) 
253    {
254      r->valptr = tipValPtr(q);         /* Assign value */
255    }
256  else 
257    {                              /* Internal node, look at children */
258      s = q->next;                      /* First child */
259      do 
260        {
261          t = saveSubtree(s, tpl, numsp, numBranches);        /* Generate child's subtree */
262
263          t0 = 0;                         /* Merge child into list */
264          t1 = r->descend;
265          while (t1 && (cmpTipVal(r0[t1].valptr, r0[t].valptr) < 0)) {
266            t0 = t1;
267            t1 = r0[t1].sibling;
268          }
269          if (t0) r0[t0].sibling = t;  else  r->descend = t;
270          r0[t].sibling = t1;
271
272          s = s->next;                    /* Next child */
273        } while (s != q);
274
275      r->valptr = r0[r->descend].valptr;   /* Inherit first child's value */
276      }                                 /* End of internal node processing */
277
278  return  (r - r0);
279}
280
281
282static nodeptr minSubtreeTip (nodeptr  p0, int numsp)
283{ 
284  nodeptr  minTip, p, testTip;
285
286  if (isTip(p0->number, numsp)) 
287    return p0;
288
289  p = p0->next;
290
291  minTip = minSubtreeTip(p->back, numsp);
292
293  while ((p = p->next) != p0) 
294    {
295      testTip = minSubtreeTip(p->back, numsp);
296      if (cmpTipVal(tipValPtr(testTip), tipValPtr(minTip)) < 0)
297        minTip = testTip;
298    }
299  return minTip;
300} 
301
302
303static nodeptr  minTreeTip (nodeptr  p, int numsp)
304{
305  nodeptr  minp, minpb;
306
307  minp  = minSubtreeTip(p, numsp);
308  minpb = minSubtreeTip(p->back, numsp);
309  return (cmpTipVal(tipValPtr(minp), tipValPtr(minpb)) < 0 ? minp : minpb);
310}
311
312
313static void saveTree (tree *tr, topol *tpl)
314/*  Save a tree topology in a standard order so that first branches
315 *  from a node contain lower value tips than do second branches from
316 *  the node.  The root tip should have the lowest value of all.
317 */
318{
319  connptr  r; 
320 
321  tpl->nextlink = 0;                             /* Reset link pointer */
322  r = tpl->links + saveSubtree(minTreeTip(tr->start, tr->rdta->numsp), tpl, tr->rdta->numsp, tr->numBranches);  /* Save tree */
323  r->sibling = 0;
324 
325  tpl->likelihood = tr->likelihood;
326  tpl->start      = tr->start;
327  tpl->ntips      = tr->ntips;
328  tpl->nextnode   = tr->nextnode;   
329 
330} /* saveTree */
331
332
333static boolean restoreTree (topol *tpl, tree *tr)
334{ 
335  connptr  r;
336  nodeptr  p, p0;   
337  int  i;
338
339  for (i = 1; i <= 2*(tr->mxtips) - 2; i++) 
340    { 
341      /* Uses p = p->next at tip */
342      p0 = p = tr->nodep[i];
343      do 
344        {
345          p->back = (nodeptr) NULL;
346          p = p->next;
347        } 
348      while (p != p0);
349    }
350
351  /*  Copy connections from topology */
352
353  for (r = tpl->links, i = 0; i < tpl->nextlink; r++, i++)     
354    hookup(r->p, r->q, r->z, tr->numBranches);     
355
356  tr->likelihood = tpl->likelihood;
357  tr->start      = tpl->start;
358  tr->ntips      = tpl->ntips;
359 
360  tr->nextnode   = tpl->nextnode;   
361
362  onlyInitrav(tr, tr->start);
363  return TRUE;
364}
365
366
367
368
369int initBestTree (bestlist *bt, int newkeep, int numsp)
370{ /* initBestTree */
371  int  i;
372
373  bt->nkeep = 0;
374
375  if (bt->ninit <= 0) 
376    {
377      if (! (bt->start = setupTopol(numsp)))  return  0;
378      bt->ninit = -1;
379      bt->nvalid = 0;
380      bt->numtrees = 0;
381      bt->best = unlikely;
382      bt->improved = FALSE;
383      bt->byScore = (topol **) rax_malloc((newkeep+1) * sizeof(topol *));
384      bt->byTopol = (topol **) rax_malloc((newkeep+1) * sizeof(topol *));
385      if (! bt->byScore || ! bt->byTopol) {
386        printf( "initBestTree: rax_malloc failure\n");
387        return 0;
388      }
389    }
390  else if (ABS(newkeep) > bt->ninit) {
391    if (newkeep <  0) newkeep = -(bt->ninit);
392    else newkeep = bt->ninit;
393  }
394
395  if (newkeep < 1) {    /*  Use negative newkeep to clear list  */
396    newkeep = -newkeep;
397    if (newkeep < 1) newkeep = 1;
398    bt->nvalid = 0;
399    bt->best = unlikely;
400  }
401 
402  if (bt->nvalid >= newkeep) {
403    bt->nvalid = newkeep;
404    bt->worst = bt->byScore[newkeep]->likelihood;
405  }
406  else 
407    {
408      bt->worst = unlikely;
409    }
410 
411  for (i = bt->ninit + 1; i <= newkeep; i++) 
412    {   
413      if (! (bt->byScore[i] = setupTopol(numsp)))  break;
414      bt->byTopol[i] = bt->byScore[i];
415      bt->ninit = i;
416    }
417 
418  return  (bt->nkeep = MIN(newkeep, bt->ninit));
419} /* initBestTree */
420
421
422
423void resetBestTree (bestlist *bt)
424{ /* resetBestTree */
425  bt->best     = unlikely;
426  bt->worst    = unlikely;
427  bt->nvalid   = 0;
428  bt->improved = FALSE;
429} /* resetBestTree */
430
431
432boolean  freeBestTree(bestlist *bt)
433{ /* freeBestTree */
434  while (bt->ninit >= 0)  freeTopol(bt->byScore[(bt->ninit)--]);
435   
436  /* VALGRIND */
437
438  rax_free(bt->byScore);
439  rax_free(bt->byTopol);
440
441  /* VALGRIND END */
442
443  freeTopol(bt->start);
444  return TRUE;
445} /* freeBestTree */
446
447
448/*  Compare two trees, assuming that each is in standard order.  Return
449 *  -1 if first preceeds second, 0 if they are identical, or +1 if first
450 *  follows second in standard order.  Lower number tips preceed higher
451 *  number tips.  A tip preceeds a corresponding internal node.  Internal
452 *  nodes are ranked by their lowest number tip.
453 */
454
455static int  cmpSubtopol (connptr p10, connptr p1, connptr p20, connptr p2)
456{
457  connptr  p1d, p2d;
458  int  cmp;
459 
460  if (! p1->descend && ! p2->descend)          /* Two tips */
461    return cmpTipVal(p1->valptr, p2->valptr);
462 
463  if (! p1->descend) return -1;                /* p1 = tip, p2 = node */
464  if (! p2->descend) return  1;                /* p2 = tip, p1 = node */
465 
466  p1d = p10 + p1->descend;
467  p2d = p20 + p2->descend;
468  while (1) {                                  /* Two nodes */
469    if ((cmp = cmpSubtopol(p10, p1d, p20, p2d)))  return cmp; /* Subtrees */
470    if (! p1d->sibling && ! p2d->sibling)  return  0; /* Lists done */
471    if (! p1d->sibling) return -1;             /* One done, other not */
472    if (! p2d->sibling) return  1;             /* One done, other not */
473    p1d = p10 + p1d->sibling;                  /* Neither done */
474    p2d = p20 + p2d->sibling;
475  }
476}
477
478
479
480static int  cmpTopol (void *tpl1, void *tpl2)
481{ 
482  connptr  r1, r2;
483  int      cmp;   
484 
485  r1 = ((topol *) tpl1)->links;
486  r2 = ((topol *) tpl2)->links;
487  cmp = cmpTipVal(tipValPtr(r1->p), tipValPtr(r2->p));
488  if (cmp)             
489    return cmp;     
490  return  cmpSubtopol(r1, r1, r2, r2);
491} 
492
493
494
495static int  cmpTplScore (void *tpl1, void *tpl2)
496{ 
497  double  l1, l2;
498 
499  l1 = ((topol *) tpl1)->likelihood;
500  l2 = ((topol *) tpl2)->likelihood;
501  return  (l1 > l2) ? -1 : ((l1 == l2) ? 0 : 1);
502}
503
504
505
506/*  Find an item in a sorted list of n items.  If the item is in the list,
507 *  return its index.  If it is not in the list, return the negative of the
508 *  position into which it should be inserted.
509 */
510
511static int  findInList (void *item, void *list[], int n, int (* cmpFunc)(void *, void *))
512{
513  int  mid, hi, lo, cmp = 0;
514 
515  if (n < 1) return  -1;                    /*  No match; first index  */
516 
517  lo = 1;
518  mid = 0;
519  hi = n;
520  while (lo < hi) {
521    mid = (lo + hi) >> 1;
522    cmp = (* cmpFunc)(item, list[mid-1]);
523    if (cmp) {
524      if (cmp < 0) hi = mid;
525      else lo = mid + 1;
526    }
527    else  return  mid;                        /*  Exact match  */
528  }
529 
530  if (lo != mid) {
531    cmp = (* cmpFunc)(item, list[lo-1]);
532    if (cmp == 0) return lo;
533  }
534  if (cmp > 0) lo++;                         /*  Result of step = 0 test  */
535  return  -lo;
536} 
537
538
539
540static int  findTreeInList (bestlist *bt, tree *tr)
541{
542  topol  *tpl;
543 
544  tpl = bt->byScore[0];
545  saveTree(tr, tpl);
546  return  findInList((void *) tpl, (void **) (& (bt->byTopol[1])),
547                     bt->nvalid, cmpTopol);
548} 
549
550
551int  saveBestTree (bestlist *bt, tree *tr)
552{   
553  topol  *tpl, *reuse;
554  int  tplNum, scrNum, reuseScrNum, reuseTplNum, i, oldValid, newValid;
555 
556  tplNum = findTreeInList(bt, tr);
557  tpl = bt->byScore[0];
558  oldValid = newValid = bt->nvalid;
559 
560  if (tplNum > 0) {                      /* Topology is in list  */
561    reuse = bt->byTopol[tplNum];         /* Matching topol  */
562    reuseScrNum = reuse->scrNum;
563    reuseTplNum = reuse->tplNum;
564  }
565  /* Good enough to keep? */
566  else if (tr->likelihood < bt->worst)  return 0;
567 
568  else {                                 /* Topology is not in list */
569    tplNum = -tplNum;                    /* Add to list (not replace) */
570    if (newValid < bt->nkeep) bt->nvalid = ++newValid;
571    reuseScrNum = newValid;              /* Take worst tree */
572    reuse = bt->byScore[reuseScrNum];
573    reuseTplNum = (newValid > oldValid) ? newValid : reuse->tplNum;
574    if (tr->likelihood > bt->start->likelihood) bt->improved = TRUE;
575  }
576 
577  scrNum = findInList((void *) tpl, (void **) (& (bt->byScore[1])),
578                      oldValid, cmpTplScore);
579  scrNum = ABS(scrNum);
580 
581  if (scrNum < reuseScrNum)
582    for (i = reuseScrNum; i > scrNum; i--)
583      (bt->byScore[i] = bt->byScore[i-1])->scrNum = i;
584 
585  else if (scrNum > reuseScrNum) {
586    scrNum--;
587    for (i = reuseScrNum; i < scrNum; i++)
588      (bt->byScore[i] = bt->byScore[i+1])->scrNum = i;
589  }
590 
591  if (tplNum < reuseTplNum)
592    for (i = reuseTplNum; i > tplNum; i--)
593      (bt->byTopol[i] = bt->byTopol[i-1])->tplNum = i;
594 
595  else if (tplNum > reuseTplNum) {
596    tplNum--;
597    for (i = reuseTplNum; i < tplNum; i++)
598      (bt->byTopol[i] = bt->byTopol[i+1])->tplNum = i;
599  }
600 
601 
602 
603  tpl->scrNum = scrNum;
604  tpl->tplNum = tplNum;
605  bt->byTopol[tplNum] = bt->byScore[scrNum] = tpl;
606  bt->byScore[0] = reuse;
607 
608  if (scrNum == 1)  bt->best = tr->likelihood;
609  if (newValid == bt->nkeep) bt->worst = bt->byScore[newValid]->likelihood;
610 
611  return  scrNum;
612} 
613
614
615int  recallBestTree (bestlist *bt, int rank, tree *tr)
616{ 
617  if (rank < 1)  rank = 1;
618  if (rank > bt->nvalid)  rank = bt->nvalid;
619  if (rank > 0)  if (! restoreTree(bt->byScore[rank], tr)) return FALSE;
620  return  rank;
621}
622
623
624
625
Note: See TracBrowser for help on using the repository browser.