source: trunk/GDE/RAxML/rapidBootstrap.c

Last change on this file was 10409, checked in by aboeckma, 12 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: 9.2 KB
Line 
1/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright August 2006 by Alexandros Stamatakis
3 *
4 *  Partially derived from
5 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6 * 
7 *  and
8 *
9 *  Programs of the PHYLIP package by Joe Felsenstein.
10 *
11 *  This program is free software; you may redistribute it and/or modify its
12 *  under the terms of the GNU General Public License as published by the Free
13 *  Software Foundation; either version 2 of the License, or (at your option)
14 *  any later version.
15 *
16 *  This program is distributed in the hope that it will be useful, but
17 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19 *  for more details.
20 *
21 *
22 *  For any other enquiries send an Email to Alexandros Stamatakis
23 *  Alexandros.Stamatakis@epfl.ch
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 *
27 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
28 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
29 */
30
31#ifndef WIN32
32#include <sys/times.h>
33#include <sys/types.h>
34#include <sys/time.h>
35#include <unistd.h>
36#endif
37
38#include <math.h>
39#include <time.h>
40#include <stdlib.h>
41#include <stdio.h>
42#include <ctype.h>
43#include <string.h>
44
45#include "axml.h"
46
47extern int Thorough;
48extern infoList iList;
49extern char seq_file[1024];
50
51
52
53
54
55
56
57
58
59
60
61void computeBOOTRAPID (tree *tr, analdef *adef, long *radiusSeed) 
62{ 
63  int i, bestTrav, impr;               
64  double lh, previousLh, difference, epsilon;             
65  bestlist *bestT, *bt; 
66  int countIT;
67 
68  bestT = (bestlist *) rax_malloc(sizeof(bestlist));
69  bestT->ninit = 0;
70  initBestTree(bestT, 1, tr->mxtips);
71  saveBestTree(bestT, tr);
72
73  bt = (bestlist *) rax_malloc(sizeof(bestlist));     
74  bt->ninit = 0; 
75  initBestTree(bt, 5, tr->mxtips);
76
77  initInfoList(10);
78 
79  difference = 10.0;
80  epsilon = 0.01;         
81
82  bestTrav = adef->bestTrav = 5 + 11 * randum(radiusSeed);
83   
84  Thorough = 1; 
85 
86  impr = 1;
87
88  if(tr->doCutoff)
89    tr->itCount = 0;
90
91  tr->bigCutoff = TRUE;
92 
93  for(countIT = 0; countIT < 2 && impr; countIT++) 
94    {             
95      recallBestTree(bestT, 1, tr);       
96      treeEvaluate(tr, 1);                                 
97      saveBestTree(bestT, tr); 
98           
99      lh = previousLh = tr->likelihood;
100         
101      treeOptimizeRapid(tr, 1, bestTrav, adef, bt);       
102
103      impr = 0;
104         
105      for(i = 1; i <= bt->nvalid; i++)
106        {                                 
107          recallBestTree(bt, i, tr);                     
108         
109
110          treeEvaluate(tr, 0.25);                 
111             
112          difference = ((tr->likelihood > previousLh)? 
113                        tr->likelihood - previousLh: 
114                        previousLh - tr->likelihood);       
115          if(tr->likelihood > lh && difference > epsilon)
116            {
117              impr = 1;       
118              lh = tr->likelihood;                   
119              saveBestTree(bestT, tr);
120            }             
121        }
122
123    } 
124 
125  tr->bigCutoff = FALSE;
126
127  recallBestTree(bestT, 1, tr);   
128  freeBestTree(bestT);
129  rax_free(bestT);
130  freeBestTree(bt);
131  rax_free(bt);
132  freeInfoList(); 
133}
134
135
136
137
138
139
140
141void optimizeRAPID(tree *tr, analdef *adef) 
142{ 
143  int i,  impr, bestTrav;
144   
145  double lh, previousLh, difference, epsilon;             
146  bestlist *bestT, *bt;   
147 
148  bestT = (bestlist *) rax_malloc(sizeof(bestlist));
149  bestT->ninit = 0;
150  initBestTree(bestT, 1, tr->mxtips);
151     
152  bt = (bestlist *) rax_malloc(sizeof(bestlist));     
153  bt->ninit = 0;
154  initBestTree(bt, 20, tr->mxtips); 
155
156  initInfoList(50);
157 
158  difference = 10.0;
159  epsilon = 0.01;   
160   
161  Thorough = 0; 
162
163 
164         
165  saveBestTree(bestT, tr); 
166  bestTrav = adef->bestTrav = determineRearrangementSetting(tr, adef, bestT, bt);                   
167  saveBestTree(bestT, tr); 
168
169  impr = 1;
170  if(tr->doCutoff)
171    tr->itCount = 0;
172
173  while(impr)
174    {             
175      recallBestTree(bestT, 1, tr);     
176      treeEvaluate(tr, 1);
177
178
179      saveBestTree(bestT, tr);     
180         
181      lh = previousLh = tr->likelihood;
182         
183      treeOptimizeRapid(tr, 1, bestTrav, adef, bt);   
184     
185      impr = 0;
186         
187      for(i = 1; i <= bt->nvalid; i++)
188        {                                 
189          recallBestTree(bt, i, tr);       
190          treeEvaluate(tr, 0.25);                       
191
192          difference = ((tr->likelihood > previousLh)? 
193                        tr->likelihood - previousLh: 
194                        previousLh - tr->likelihood);       
195          if(tr->likelihood > lh && difference > epsilon)
196            {
197              impr = 1;       
198              lh = tr->likelihood;                   
199              saveBestTree(bestT, tr);
200            }             
201        }         
202    }
203
204  recallBestTree(bestT, 1, tr);
205  freeBestTree(bestT);
206  rax_free(bestT);
207  freeBestTree(bt);
208  rax_free(bt);
209  freeInfoList(); 
210}
211
212
213
214void thoroughOptimization(tree *tr, analdef *adef, topolRELL_LIST *rl, int index) 
215{ 
216  int i, impr;               
217  int rearrangementsMin = 1, rearrangementsMax = adef->stepwidth;
218  double lh, previousLh, difference, epsilon;             
219  bestlist *bestT, *bt; 
220   
221  bestT = (bestlist *) rax_malloc(sizeof(bestlist));
222  bestT->ninit = 0;
223  initBestTree(bestT, 1, tr->mxtips);
224     
225  bt = (bestlist *) rax_malloc(sizeof(bestlist));     
226  bt->ninit = 0;   
227  initBestTree(bt, 20, tr->mxtips);
228
229  initInfoList(50);
230 
231  difference = 10.0;
232  epsilon = 0.01;       
233 
234 
235  saveBestTree(bestT, tr); 
236 
237  impr = 1;
238  if(tr->doCutoff)
239    tr->itCount = 0;
240 
241  Thorough = 1;
242  impr = 1;
243
244  while(1)
245    {           
246      recallBestTree(bestT, 1, tr);   
247      if(impr)
248        {               
249          rearrangementsMin = 1;
250          rearrangementsMax = adef->stepwidth;           
251        }                                               
252      else
253        {                         
254          rearrangementsMax += adef->stepwidth;
255          rearrangementsMin += adef->stepwidth;                               
256          if(rearrangementsMax > adef->max_rearrange)                   
257            goto cleanup;         
258        }
259         
260      treeEvaluate(tr, 1.0);         
261      previousLh = lh = tr->likelihood;       
262      saveBestTree(bestT, tr);                           
263     
264      treeOptimizeRapid(tr, rearrangementsMin, rearrangementsMax, adef, bt);
265     
266      impr = 0;                         
267               
268      for(i = 1; i <= bt->nvalid; i++)
269        {       
270          recallBestTree(bt, i, tr);       
271                 
272          treeEvaluate(tr, 0.25);               
273         
274          difference = ((tr->likelihood > previousLh)? 
275                        tr->likelihood - previousLh: 
276                        previousLh - tr->likelihood);       
277          if(tr->likelihood > lh && difference > epsilon)
278            {
279              impr = 1;       
280              lh = tr->likelihood;                   
281              saveBestTree(bestT, tr);
282            }             
283        }
284
285       
286    }
287
288 cleanup: 
289  saveTL(rl, tr, index); 
290  freeBestTree(bestT);
291  rax_free(bestT);
292  freeBestTree(bt);
293  rax_free(bt);
294  freeInfoList();
295}
296
297
298
299
300/*********************************************************************************************************************/
301
302
303
304
305
306
307
308
309
310
311static boolean qupdate (tree *tr, nodeptr p)
312{   
313  nodeptr  q;
314  double   z0[NUM_BRANCHES], z[NUM_BRANCHES];   
315  int i;
316     
317  q = p->back;
318  for(i = 0; i < tr->numBranches; i++)
319    z0[i] = q->z[i];
320     
321  makenewzGeneric(tr, p, q, z0, 1, z, FALSE);   
322     
323  for(i = 0; i < tr->numBranches; i++)
324    p->z[i] = q->z[i] = z[i]; 
325       
326  return TRUE;   
327} 
328
329
330
331
332
333static boolean qsmoothLocal(tree *tr, nodeptr p, int n)
334{
335  nodeptr  q;
336 
337  if(n == 0)
338    return TRUE;
339  else
340    {
341      if (! qupdate(tr, p))               return FALSE; /*  Adjust branch */
342      if (!isTip(p->number, tr->rdta->numsp)) 
343        {                                  /*  Adjust descendants */
344          q = p->next;
345          while (q != p) 
346            {
347              if (! qsmoothLocal(tr, q->back, n - 1))   return FALSE;
348              q = q->next;
349            }   
350         
351          newviewGeneric(tr, p);
352        }
353     
354      return TRUE;
355    }
356} 
357
358static void quickSmoothLocal(tree *tr, int n)
359{
360  nodeptr p = tr->insertNode;
361  nodeptr q;
362
363  if(n == 0)
364    {
365      evaluateGeneric(tr, p);
366    }
367  else
368    {
369      qsmoothLocal(tr, p->back, n - 1);
370      if(!isTip(p->number, tr->rdta->numsp))
371        {
372          q = p->next;
373          while(q != p)
374            {
375              qsmoothLocal(tr, q->back, n - 1);
376              q = q->next;
377            }
378        }
379      evaluateGeneric(tr, p);
380    }
381}
382
383
384
385
386int treeOptimizeThorough(tree *tr, int mintrav, int maxtrav)
387{
388  int i;   
389  bestlist *bestT; 
390
391  nodeRectifier(tr);
392
393  bestT = (bestlist *) rax_malloc(sizeof(bestlist));
394  bestT->ninit = 0;
395  initBestTree(bestT, 1, tr->mxtips);
396 
397  if (maxtrav > tr->ntips - 3) 
398    maxtrav = tr->ntips - 3;     
399 
400  tr->startLH = tr->endLH = tr->likelihood; 
401
402  for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++)
403    {     
404
405     
406      tr->bestOfNode = unlikely;     
407      if(rearrangeBIG(tr, tr->nodep[i], mintrav, maxtrav))
408        {         
409         
410          if((tr->endLH > tr->startLH) && (tr->bestOfNode != unlikely))
411            {                       
412              restoreTreeFast(tr);           
413              quickSmoothLocal(tr, 3);
414              tr->startLH = tr->endLH = tr->likelihood;                     
415            }   
416          else
417            {           
418              if(tr->bestOfNode != unlikely)
419                {                   
420                  resetBestTree(bestT);                                           
421                  saveBestTree(bestT, tr);
422                  restoreTreeFast(tr);           
423                  quickSmoothLocal(tr, 3);                                 
424                  if(tr->likelihood < tr->startLH)                                 
425                    {
426                      int res;
427                      res = recallBestTree(bestT, 1, tr);                                   
428                      assert(res > 0);
429                    }
430                  else                             
431                    tr->startLH = tr->endLH = tr->likelihood;             
432                }
433            }
434        }
435
436       
437    }   
438
439  freeBestTree(bestT);
440  rax_free(bestT);
441
442  return 1;     
443}
Note: See TracBrowser for help on using the repository browser.