source: branches/profile/GDE/PHYML/simu.c

Last change on this file was 4073, checked in by westram, 13 years ago
  • phyml 2.4.5
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.4 KB
Line 
1/*
2
3PHYML :  a program that  computes maximum likelihood  phylogenies from
4DNA or AA homologous sequences
5
6Copyright (C) Stephane Guindon. Oct 2003 onward
7
8All parts of  the source except where indicated  are distributed under
9the GNU public licence.  See http://www.opensource.org for details.
10
11*/
12
13#include "utilities.h"
14#include "lk.h"
15#include "optimiz.h"
16#include "models.h"
17#include "free.h"
18#include "simu.h"
19
20
21/* int    BRENT_ITMAX; */
22/* double BRENT_CGOLD; */
23/* double BRENT_ZEPS; */
24/* double MNBRAK_GOLD; */
25/* double MNBRAK_GLIMIT; */
26/* double MNBRAK_TINY; */
27/* double ALPHA_MIN; */
28/* double ALPHA_MAX; */
29/* double BL_MIN; */
30/* double BL_START; */
31/* double BL_MAX; */
32/* double MIN_DIFF_LK; */
33/* double GOLDEN_R; */
34/* double GOLDEN_C; */
35/* int    T_MAX_FILE; */
36/* int    T_MAX_LINE; */
37/* int    T_MAX_NAME; */
38/* int    T_MAX_SEQ; */
39/* int    N_MAX_INSERT; */
40/* int    N_MAX_OTU; */
41/* double UNLIKELY; */
42/* double NJ_SEUIL; */
43/* int    MAX_TOPO_DIST; */
44/* double DIST_MAX; */
45/* int    LIM_SCALE; */
46/* double AROUND_LK; */
47/* double PROP_STEP; */
48
49
50/*********************************************************/
51
52void Simu(arbre *tree, int n_step_max)
53{
54    double old_loglk,n_iter,lambda,diff_lk;
55    int i,n_neg,n_tested,n_without_swap,step,it_lim_without_swap;
56    edge **sorted_b,**tested_b;
57    int each,each_invar;
58    int opt_free_param;
59   
60    sorted_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
61    tested_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
62   
63   
64   
65    old_loglk = tree->tot_loglk = UNLIKELY;
66    n_iter              = 1.0;
67    it_lim_without_swap = (tree->mod->invar)?(8):(5);
68    n_tested            = 0;
69    n_without_swap      = 0;
70    step                = 0;
71    each                = 4;
72    lambda              = 0.75;
73    each_invar          = 2;
74    old_loglk           = tree->tot_loglk;
75    opt_free_param      = 0;
76
77    do
78        {
79            ++step;
80           
81            if(step > n_step_max) break;
82           
83            each--;
84            each_invar--;
85           
86            tree->mod->s_opt->opt_bl = 0;
87            tree->both_sides    = 1;
88            Lk(tree,tree->data);
89           
90           
91            if(tree->mod->s_opt->print)
92                {
93                    if(old_loglk < UNLIKELY+1)
94                        printf("\n. Log(lk) :               * -> %15.6f ",tree->tot_loglk);
95                    else
96                        printf("\n. Log(lk) : %15.6f -> %15.6f ",old_loglk,tree->tot_loglk);
97                   
98                    if(old_loglk > UNLIKELY+1)
99                        {
100                            if(n_tested > 1) printf("%3d swaps done ",n_tested);
101                            else             printf("%3d swap  done", n_tested);
102                        }
103                }
104           
105            fflush(NULL);
106           
107            if((fabs(old_loglk-tree->tot_loglk) < 1.E-03) || (n_without_swap > it_lim_without_swap)) break;
108           
109            diff_lk = old_loglk-tree->tot_loglk;
110           
111            if(tree->tot_loglk < old_loglk)
112                {
113                    if(tree->mod->s_opt->print)
114                        printf("\n\n. Moving backward (topology + branch lengths) \n");
115                    fflush(NULL);
116                    if(!Mov_Backward_Topo_Bl(tree,old_loglk,tested_b,n_tested))
117                        Exit("\n. Err: mov_back failed\n");
118                    if(!tree->n_swap) n_neg = 0;
119                   
120                    For(i,2*tree->n_otu-3) tree->t_edges[i]->l_old = tree->t_edges[i]->l;
121                   
122                    Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
123                }
124            else 
125                {
126                    if(!each)
127                        {
128                            opt_free_param = 1;
129                            each = 4;
130                            if(tree->mod->s_opt->print) printf("\n");
131                            Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
132                            tree->mod->s_opt->opt_bl = 0;
133                            tree->both_sides    = 1;
134                            Lk(tree,tree->data);
135                        }
136                   
137                    old_loglk = tree->tot_loglk;
138                   
139                    Fix_All(tree);
140                   
141                    n_neg = 0;
142                    For(i,2*tree->n_otu-3)
143                        if((!tree->t_edges[i]->left->tax) && 
144                           (!tree->t_edges[i]->rght->tax)) 
145                            NNI(tree,tree->t_edges[i],0);
146                   
147                    Select_Edges_To_Swap(tree,sorted_b,&n_neg);
148                   
149                    Sort_Edges_Diff_Lk(tree,sorted_b,n_neg);
150                   
151                    Optimiz_Ext_Br(tree);         
152                   
153                    Update_Bl(tree,lambda);
154                   
155                    n_tested = 0;
156                    For(i,(int)ceil((double)n_neg*(lambda)))
157                        tested_b[n_tested++] = sorted_b[i];
158                   
159                    Make_N_Swap(tree,tested_b,0,n_tested);
160                   
161                    if(n_tested > 0) n_without_swap = 0;
162                    else             n_without_swap++;
163                   
164                    fflush(NULL);
165                }
166            n_iter+=1.0;
167        }
168    while(1);
169   
170   
171    Free(sorted_b);
172    Free(tested_b);
173   
174   
175    if((n_without_swap > it_lim_without_swap) || (tree->mod->s_opt->last_opt))
176      {
177        if(tree->mod->s_opt->print)
178          printf("\n\n. Last optimization step...\n"); fflush(NULL);
179        Round_Optimize(tree,tree->data);
180      }
181}
182
183/*********************************************************/
184
185void Select_Edges_To_Swap(arbre *tree, edge **sorted_b, int *n_neg)
186{
187  int i;
188  edge *b;
189  int min;
190
191
192  *n_neg = 0;
193  tree->min_diff_lk = .0;
194  min = 0;
195  For(i,2*tree->n_otu-3)
196    {
197      b = tree->t_edges[i];
198
199      if((!b->left->tax) 
200      && (!b->rght->tax) 
201      && (b->diff_lk < 0.0-MDBL_MIN)) 
202        {
203
204          if((b->left->b[b->l_v1]->diff_lk < b->diff_lk) ||
205             (b->left->b[b->l_v2]->diff_lk < b->diff_lk) ||
206             (b->rght->b[b->r_v1]->diff_lk < b->diff_lk) ||
207             (b->rght->b[b->r_v2]->diff_lk < b->diff_lk)) continue;
208
209          if(b->diff_lk < tree->min_diff_lk) 
210            {
211              tree->min_diff_lk = b->diff_lk;
212              min = i;
213            }
214
215          sorted_b[*n_neg] = b;
216          (*n_neg)++;
217        }
218    }
219}
220
221/*********************************************************/
222
223void Fix_All(arbre *tree)
224{
225  int i;
226  edge *b;
227
228  tree->mod->pinvar_old = tree->mod->pinvar;
229  tree->mod->alpha_old = tree->mod->alpha;
230  tree->mod->kappa_old = tree->mod->kappa;
231  tree->mod->lambda_old = tree->mod->lambda;
232 
233  For(i,2*tree->n_otu-3)
234    {
235      b = tree->t_edges[i];
236      b->l_old = b->l;
237    }
238}
239
240/*********************************************************/
241
242void Update_Bl(arbre *tree, double fact)
243{
244  int i;
245  edge *b;
246
247  For(i,2*tree->n_otu-3)
248    {
249      b = tree->t_edges[i];
250      b->l = b->l_old + (b->ql[0]-b->l_old)*fact;
251    }
252}
253
254/*********************************************************/
255
256void Make_N_Swap(arbre *tree,edge **b, int beg, int end)
257{
258  int i;
259
260  tree->n_swap = 0;
261  for(i=beg;i<end;i++)
262    {
263/*       printf("make swap on %3d d->%10f\n",b[i]->num,b[i]->diff_lk); */
264/*       if(drand48()>0.75) */
265/*      { */
266        (b[i]->best_conf == 2)?
267          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)):
268          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree));
269       
270        b[i]->l = b[i]->ql[b[i]->best_conf-1];
271        tree->n_swap++;
272/*      } */
273    }
274}
275
276/*********************************************************/
277
278int Make_Best_Swap(arbre *tree)
279{
280  int i,j,return_value;
281  edge *b,**sorted_b;
282 
283
284  sorted_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
285 
286  j=0;
287  For(i,2*tree->n_otu-3) if((!tree->t_edges[i]->left->tax) &&
288                            (!tree->t_edges[i]->rght->tax))
289                              sorted_b[j++] = tree->t_edges[i];
290
291  Sort_Edges_Diff_Lk(tree,sorted_b,tree->n_otu-3);
292
293  if(sorted_b[0]->diff_lk < -0.0)
294    {
295      b = sorted_b[0];
296      return_value = 1;
297      (b->best_conf == 2)?
298        (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v1],tree)):
299        (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v2],tree));
300     
301      b->l = b->ql[b->best_conf-1];
302    }
303  else return_value = 0;
304
305  Free(sorted_b);
306
307  return return_value;
308}
309
310/*********************************************************/
311
312int Mov_Backward_Topo_Bl(arbre *tree, double lk_old, edge **tested_b, int n_tested)
313{
314  double *l_init;
315  int i,step,n_swp,beg,end;
316  edge *b,**swp;
317
318
319  l_init = (double *)mCalloc(2*tree->n_otu-3,sizeof(double));
320  swp = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
321
322  For(i,2*tree->n_otu-3) l_init[i] = tree->t_edges[i]->l;
323 
324  step = 2;
325  tree->both_sides = 0;
326  do
327    {
328      n_swp = 0;
329      For(i,2*tree->n_otu-3) 
330        {
331          b = tree->t_edges[i];
332          b->l = b->l_old + (1./step) * (l_init[i] - b->l_old);
333        }
334
335      beg = (int)floor((double)n_tested/(step-1));
336      end = 0;
337      Unswap_N_Branch(tree,tested_b,beg,end);
338      beg = 0;
339      end = (int)floor((double)n_tested/step);
340      Swap_N_Branch(tree,tested_b,beg,end);
341     
342      if(end == n_swp) tree->n_swap = 0;
343     
344      tree->mod->s_opt->opt_bl = 0;
345      tree->both_sides    = 0;
346      Lk(tree,tree->data);
347     
348      step++;
349
350    }while((tree->tot_loglk < lk_old) && (step < 100));
351
352
353  if(step == 100)
354    {
355      For(i,2*tree->n_otu-3) 
356        {
357          b = tree->t_edges[i];
358          b->l = b->l_old;
359        }
360
361      tree->mod->s_opt->opt_bl = 0;
362      tree->both_sides    = 0;
363      Lk(tree,tree->data);
364    }
365
366  Free(l_init);
367  Free(swp);
368
369  tree->n_swap = 0;
370  For(i,2*tree->n_otu-3) 
371    {
372      if(tree->t_edges[i]->diff_lk < 0.0) tree->n_swap++;
373      tree->t_edges[i]->diff_lk = +1.0;
374    }
375
376  if(tree->tot_loglk > lk_old)                 return  1;
377  else if((tree->tot_loglk > lk_old-MIN_DIFF_LK) && 
378          (tree->tot_loglk < lk_old+MIN_DIFF_LK)) return -1;
379  else                                         return  0;
380}
381
382/*********************************************************/
383
384void Unswap_N_Branch(arbre *tree, edge **b, int beg, int end)
385{
386  int i;
387 
388  if(end>beg)
389    {
390      for(i=beg;i<end;i++)
391        {
392          (b[i]->best_conf == 2)?
393            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)):
394            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree));
395          b[i]->l = b[i]->l_old;
396        }
397    }
398  else
399    {
400      for(i=beg-1;i>=end;i--)
401        {
402          (b[i]->best_conf == 2)?
403            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)):
404            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree));
405          b[i]->l = b[i]->l_old;
406        }
407    }
408}
409
410/*********************************************************/
411
412void Swap_N_Branch(arbre *tree,edge **b, int beg, int end)
413{
414  int i;
415 
416  if(end>beg)
417    {
418      for(i=beg;i<end;i++)
419        {
420          (b[i]->best_conf == 2)?
421            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)):
422            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree));
423          b[i]->l = b[i]->ql[b[i]->best_conf-1];
424        }
425    }
426  else
427    {
428      for(i=beg-1;i>=end;i--)
429        {
430          (b[i]->best_conf == 2)?
431            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)):
432            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree));
433          b[i]->l = b[i]->ql[b[i]->best_conf-1];
434        }
435
436    }
437}
438
439/*********************************************************/
Note: See TracBrowser for help on using the repository browser.