source: trunk/GDE/PHYML/simu.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.1 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;
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   
59    sorted_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
60    tested_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
61   
62   
63   
64    old_loglk = tree->tot_loglk = UNLIKELY;
65    n_iter              = 1.0;
66    it_lim_without_swap = (tree->mod->invar)?(8):(5);
67    n_tested            = 0;
68    n_without_swap      = 0;
69    step                = 0;
70    each                = 4;
71    lambda              = 0.75;
72    each_invar          = 2;
73    old_loglk           = tree->tot_loglk;
74
75    do
76        {
77            ++step;
78           
79            if(step > n_step_max) break;
80           
81            each--;
82            each_invar--;
83           
84            tree->mod->s_opt->opt_bl = 0;
85            tree->both_sides    = 1;
86            Lk(tree,tree->data);
87           
88           
89            if(tree->mod->s_opt->print)
90                {
91                    if(old_loglk < UNLIKELY+1)
92                        printf("\n. Log(lk) :               * -> %15.6f ",tree->tot_loglk);
93                    else
94                        printf("\n. Log(lk) : %15.6f -> %15.6f ",old_loglk,tree->tot_loglk);
95                   
96                    if(old_loglk > UNLIKELY+1)
97                        {
98                            if(n_tested > 1) printf("%3d swaps done ",n_tested);
99                            else             printf("%3d swap  done", n_tested);
100                        }
101                }
102           
103            fflush(NULL);
104           
105            if((fabs(old_loglk-tree->tot_loglk) < 1.E-03) || (n_without_swap > it_lim_without_swap)) break;
106           
107            if(tree->tot_loglk < old_loglk)
108                {
109                    if(tree->mod->s_opt->print)
110                        printf("\n\n. Moving backward (topology + branch lengths) \n");
111                    fflush(NULL);
112                    if(!Mov_Backward_Topo_Bl(tree,old_loglk,tested_b,n_tested))
113                        Exit("\n. Err: mov_back failed\n");
114                    if(!tree->n_swap) n_neg = 0;
115                   
116                    For(i,2*tree->n_otu-3) tree->t_edges[i]->l_old = tree->t_edges[i]->l;
117                   
118                    Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
119                }
120            else 
121                {
122                    if(!each)
123                        {
124                            each = 4;
125                            if(tree->mod->s_opt->print) printf("\n");
126                            Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
127                            tree->mod->s_opt->opt_bl = 0;
128                            tree->both_sides    = 1;
129                            Lk(tree,tree->data);
130                        }
131                   
132                    old_loglk = tree->tot_loglk;
133                   
134                    Fix_All(tree);
135                   
136                    n_neg = 0;
137                    For(i,2*tree->n_otu-3)
138                        if((!tree->t_edges[i]->left->tax) && 
139                           (!tree->t_edges[i]->rght->tax)) 
140                            NNI(tree,tree->t_edges[i],0);
141                   
142                    Select_Edges_To_Swap(tree,sorted_b,&n_neg);
143                   
144                    Sort_Edges_Diff_Lk(sorted_b,n_neg);
145                   
146                    Optimiz_Ext_Br(tree);         
147                   
148                    Update_Bl(tree,lambda);
149                   
150                    n_tested = 0;
151                    For(i,(int)ceil((double)n_neg*(lambda)))
152                        tested_b[n_tested++] = sorted_b[i];
153                   
154                    Make_N_Swap(tree,tested_b,0,n_tested);
155                   
156                    if(n_tested > 0) n_without_swap = 0;
157                    else             n_without_swap++;
158                   
159                    fflush(NULL);
160                }
161            n_iter+=1.0;
162        }
163    while(1);
164   
165   
166    Free(sorted_b);
167    Free(tested_b);
168   
169   
170    if((n_without_swap > it_lim_without_swap) || (tree->mod->s_opt->last_opt))
171      {
172        if(tree->mod->s_opt->print)
173          printf("\n\n. Last optimization step...\n"); fflush(NULL);
174        Round_Optimize(tree,tree->data);
175      }
176}
177
178/*********************************************************/
179
180void Select_Edges_To_Swap(arbre *tree, edge **sorted_b, int *n_neg)
181{
182  int i;
183  edge *b;
184
185
186  *n_neg = 0;
187  tree->min_diff_lk = .0;
188  For(i,2*tree->n_otu-3)
189    {
190      b = tree->t_edges[i];
191
192      if((!b->left->tax) 
193      && (!b->rght->tax) 
194      && (b->diff_lk < 0.0-MDBL_MIN)) 
195        {
196
197          if((b->left->b[b->l_v1]->diff_lk < b->diff_lk) ||
198             (b->left->b[b->l_v2]->diff_lk < b->diff_lk) ||
199             (b->rght->b[b->r_v1]->diff_lk < b->diff_lk) ||
200             (b->rght->b[b->r_v2]->diff_lk < b->diff_lk)) continue;
201
202          if(b->diff_lk < tree->min_diff_lk) 
203            {
204              tree->min_diff_lk = b->diff_lk;
205            }
206
207          sorted_b[*n_neg] = b;
208          (*n_neg)++;
209        }
210    }
211}
212
213/*********************************************************/
214
215void Fix_All(arbre *tree)
216{
217  int i;
218  edge *b;
219
220  tree->mod->pinvar_old = tree->mod->pinvar;
221  tree->mod->alpha_old = tree->mod->alpha;
222  tree->mod->kappa_old = tree->mod->kappa;
223  tree->mod->lambda_old = tree->mod->lambda;
224 
225  For(i,2*tree->n_otu-3)
226    {
227      b = tree->t_edges[i];
228      b->l_old = b->l;
229    }
230}
231
232/*********************************************************/
233
234void Update_Bl(arbre *tree, double fact)
235{
236  int i;
237  edge *b;
238
239  For(i,2*tree->n_otu-3)
240    {
241      b = tree->t_edges[i];
242      b->l = b->l_old + (b->ql[0]-b->l_old)*fact;
243    }
244}
245
246/*********************************************************/
247
248void Make_N_Swap(arbre *tree,edge **b, int beg, int end)
249{
250  int i;
251
252  tree->n_swap = 0;
253  for(i=beg;i<end;i++)
254    {
255/*       printf("make swap on %3d d->%10f\n",b[i]->num,b[i]->diff_lk); */
256/*       if(drand48()>0.75) */
257/*      { */
258        (b[i]->best_conf == 2)?
259          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1])):
260          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2]));
261       
262        b[i]->l = b[i]->ql[b[i]->best_conf-1];
263        tree->n_swap++;
264/*      } */
265    }
266}
267
268/*********************************************************/
269
270int Make_Best_Swap(arbre *tree)
271{
272  int i,j,return_value;
273  edge *b,**sorted_b;
274 
275
276  sorted_b = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
277 
278  j=0;
279  For(i,2*tree->n_otu-3) if((!tree->t_edges[i]->left->tax) &&
280                            (!tree->t_edges[i]->rght->tax))
281                              sorted_b[j++] = tree->t_edges[i];
282
283  Sort_Edges_Diff_Lk(sorted_b,tree->n_otu-3);
284
285  if(sorted_b[0]->diff_lk < -0.0)
286    {
287      b = sorted_b[0];
288      return_value = 1;
289      (b->best_conf == 2)?
290        (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v1])):
291        (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v2]));
292     
293      b->l = b->ql[b->best_conf-1];
294    }
295  else return_value = 0;
296
297  Free(sorted_b);
298
299  return return_value;
300}
301
302/*********************************************************/
303
304int Mov_Backward_Topo_Bl(arbre *tree, double lk_old, edge **tested_b, int n_tested)
305{
306  double *l_init;
307  int i,step,n_swp,beg,end;
308  edge *b,**swp;
309
310
311  l_init = (double *)mCalloc(2*tree->n_otu-3,sizeof(double));
312  swp = (edge **)mCalloc(tree->n_otu-3,sizeof(edge *));
313
314  For(i,2*tree->n_otu-3) l_init[i] = tree->t_edges[i]->l;
315 
316  step = 2;
317  tree->both_sides = 0;
318  do
319    {
320      n_swp = 0;
321      For(i,2*tree->n_otu-3) 
322        {
323          b = tree->t_edges[i];
324          b->l = b->l_old + (1./step) * (l_init[i] - b->l_old);
325        }
326
327      beg = (int)floor((double)n_tested/(step-1));
328      end = 0;
329      Unswap_N_Branch(tested_b,beg,end);
330      beg = 0;
331      end = (int)floor((double)n_tested/step);
332      Swap_N_Branch(tested_b,beg,end);
333     
334      if(end == n_swp) tree->n_swap = 0;
335     
336      tree->mod->s_opt->opt_bl = 0;
337      tree->both_sides    = 0;
338      Lk(tree,tree->data);
339     
340      step++;
341
342    }while((tree->tot_loglk < lk_old) && (step < 100));
343
344
345  if(step == 100)
346    {
347      For(i,2*tree->n_otu-3) 
348        {
349          b = tree->t_edges[i];
350          b->l = b->l_old;
351        }
352
353      tree->mod->s_opt->opt_bl = 0;
354      tree->both_sides    = 0;
355      Lk(tree,tree->data);
356    }
357
358  Free(l_init);
359  Free(swp);
360
361  tree->n_swap = 0;
362  For(i,2*tree->n_otu-3) 
363    {
364      if(tree->t_edges[i]->diff_lk < 0.0) tree->n_swap++;
365      tree->t_edges[i]->diff_lk = +1.0;
366    }
367
368  if(tree->tot_loglk > lk_old)                 return  1;
369  else if((tree->tot_loglk > lk_old-MIN_DIFF_LK) && 
370          (tree->tot_loglk < lk_old+MIN_DIFF_LK)) return -1;
371  else                                         return  0;
372}
373
374/*********************************************************/
375
376void Unswap_N_Branch(edge **b, int beg, int end)
377{
378  int i;
379 
380  if(end>beg)
381    {
382      for(i=beg;i<end;i++)
383        {
384          (b[i]->best_conf == 2)?
385            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1])):
386            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2]));
387          b[i]->l = b[i]->l_old;
388        }
389    }
390  else
391    {
392      for(i=beg-1;i>=end;i--)
393        {
394          (b[i]->best_conf == 2)?
395            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1])):
396            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2]));
397          b[i]->l = b[i]->l_old;
398        }
399    }
400}
401
402/*********************************************************/
403
404void Swap_N_Branch(edge **b, int beg, int end)
405{
406  int i;
407 
408  if(end>beg)
409    {
410      for(i=beg;i<end;i++)
411        {
412          (b[i]->best_conf == 2)?
413            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1])):
414            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2]));
415          b[i]->l = b[i]->ql[b[i]->best_conf-1];
416        }
417    }
418  else
419    {
420      for(i=beg-1;i>=end;i--)
421        {
422          (b[i]->best_conf == 2)?
423            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1])):
424            (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2]));
425          b[i]->l = b[i]->ql[b[i]->best_conf-1];
426        }
427
428    }
429}
430
431/*********************************************************/
Note: See TracBrowser for help on using the repository browser.