source: branches/ali/GDE/PHYML20130708/phyml/src/simu.c

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

added most recent version of phyml
File size: 25.3 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 "simu.h"
14
15
16//////////////////////////////////////////////////////////////
17//////////////////////////////////////////////////////////////
18
19void Simu_Loop(t_tree *mixt_tree)
20{
21  phydbl lk_old;
22  int *orig_catg;
23  int n;
24  t_tree *tree,**tree_list;
25 
26
27  if((mixt_tree->mod->s_opt->print) &&
28     (!mixt_tree->io->quiet)) PhyML_Printf("\n\n. Refining the tree...\n");
29
30  /*! Get the number of classes in each mixture */
31  orig_catg = MIXT_Get_Number_Of_Classes_In_All_Mixtures(mixt_tree);
32
33
34  /*! Set the number of rate classes to (at most) 2.
35    ! Propagate this to every mixture tree in the analysis
36  */
37  tree = mixt_tree;
38  n = 0;
39  do
40    {
41      tree->mod->ras->n_catg = MIN(2,orig_catg[n]);
42      if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
43      tree = tree->next_mixt;
44      n++;
45    }
46  while(tree);
47 
48  /*! Make sure the number of trees in each mixture is at most 2
49   */
50  tree_list = MIXT_Record_All_Mixtures(mixt_tree);
51  MIXT_Break_All_Mixtures(orig_catg,mixt_tree);
52
53 
54  Set_Both_Sides(YES,mixt_tree);
55  Lk(NULL,mixt_tree);
56
57
58  mixt_tree->best_pars                     = 1E+8;
59  mixt_tree->mod->s_opt->spr_pars          = NO;
60  mixt_tree->mod->s_opt->quickdirty        = NO;
61  mixt_tree->best_lnL                      = mixt_tree->c_lnL;
62  mixt_tree->mod->s_opt->max_delta_lnL_spr = 0.;
63  mixt_tree->mod->s_opt->br_len_in_spr     = 1;
64  mixt_tree->mod->s_opt->max_depth_path    = 2*mixt_tree->n_otu-3;
65  mixt_tree->mod->s_opt->spr_lnL           = NO;
66
67
68  /*! Optimize branch lengths and substitution model parameters
69   */
70  Optimiz_All_Free_Param(mixt_tree,(mixt_tree->io->quiet)?(0):(mixt_tree->mod->s_opt->print));
71
72
73  /*! Rough SPR search
74   */
75  do
76    {
77      lk_old = mixt_tree->c_lnL;
78      Speed_Spr(mixt_tree,1);
79      if((mixt_tree->n_improvements < 20) ||
80         (mixt_tree->max_spr_depth  < 5) ||
81         (FABS(lk_old-mixt_tree->c_lnL) < 1.)) break;
82    }
83  while(1);
84
85
86  Set_Both_Sides(YES,mixt_tree);
87  Lk(NULL,mixt_tree);
88 
89  /*! Go back to the original data structure, with potentially more
90    ! than 2 trees per mixture
91   */
92  MIXT_Reconnect_All_Mixtures(tree_list,mixt_tree);
93  Free(tree_list);
94
95  /*! Set the number of rate classes to their original values
96   */
97  tree = mixt_tree;
98  n = 0;
99  do
100    {
101      tree->mod->ras->n_catg = orig_catg[n];
102      if(tree->mod->ras->invar == YES) tree->mod->ras->n_catg--;
103      tree = tree->next_mixt;
104      n++;
105    }
106  while(tree);
107 
108  Free(orig_catg);
109
110
111
112  /*! Only the first two trees for each mixture have been modified so
113    ! far -> need to update the other trees by copying the modified trees
114    ! onto them.
115   */
116  tree = mixt_tree;
117  do
118    {
119      if(tree != mixt_tree) Copy_Tree(mixt_tree,tree);
120      tree = tree->next;
121    }
122  while(tree);
123
124  if((mixt_tree->mod->s_opt->print) && (!mixt_tree->io->quiet))
125    {
126      PhyML_Printf("\n\n. End of refining stage...\n. The log-likelihood might now decrease and then increase again...\n");
127      PhyML_Printf("\n\n. Maximizing likelihood (using NNI moves)...\n");
128    }
129
130
131  Set_Both_Sides(YES,mixt_tree);
132  Lk(NULL,mixt_tree);
133  mixt_tree->best_lnL = mixt_tree->c_lnL;
134
135  int n_tot_moves = 0;
136  do
137    {
138      lk_old = mixt_tree->c_lnL;
139      Optimiz_All_Free_Param(mixt_tree,(mixt_tree->io->quiet)?(0):(mixt_tree->mod->s_opt->print));     
140      n_tot_moves = Simu(mixt_tree,10);
141      if(!n_tot_moves) break;
142    }
143  while(mixt_tree->c_lnL > lk_old + mixt_tree->mod->s_opt->min_diff_lk_local);
144
145  do
146    {
147      Round_Optimize(mixt_tree,mixt_tree->data,ROUND_MAX);
148      if(!Check_NNI_Five_Branches(mixt_tree)) break;     
149    }while(1);
150  /*****************************/
151 
152  if((mixt_tree->mod->s_opt->print) && 
153     (!mixt_tree->io->quiet)) PhyML_Printf("\n");
154}
155
156//////////////////////////////////////////////////////////////
157//////////////////////////////////////////////////////////////
158
159int Simu(t_tree *tree, int n_step_max)
160{
161  phydbl old_loglk,n_iter,lambda;
162  int i,n_neg,n_tested,n_without_swap,n_tot_swap,step,it_lim_without_swap;
163  t_edge **sorted_b,**tested_b;
164 
165  sorted_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
166  tested_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
167 
168  old_loglk           = UNLIKELY;
169  tree->c_lnL         = UNLIKELY;
170  n_iter              = 1.0;
171  /* it_lim_without_swap = (tree->mod->ras->invar)?(5):(2); */
172  it_lim_without_swap = (tree->mod->ras->invar)?(1):(1);
173  n_tested            = 0;
174  n_without_swap      = 0;
175  step                = 0;
176  lambda              = .75;
177  n_tot_swap          = 0;
178 
179  Update_Dirs(tree);
180     
181  if(tree->lock_topo)
182    {
183      PhyML_Printf("\n. The tree topology is locked.");
184      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
185      Warn_And_Exit("");
186    }
187 
188  do
189    {
190      ++step;
191                 
192      if(n_tested || step == 1) MIXT_Set_Alias_Subpatt(YES,tree);
193      old_loglk = tree->c_lnL;
194      Set_Both_Sides(YES,tree);
195      Lk(NULL,tree);
196      MIXT_Set_Alias_Subpatt(NO,tree);
197           
198      if(tree->c_lnL < old_loglk)
199        {
200          if((tree->mod->s_opt->print) && (!tree->io->quiet)) PhyML_Printf("\n\n. Moving backward\n");
201          if(!Mov_Backward_Topo_Bl(tree,old_loglk,tested_b,n_tested))
202            Exit("\n== Err. mov_back failed\n");
203          if(!tree->n_swap) n_neg = 0;
204         
205          /* For(i,2*tree->n_otu-3) tree->a_edges[i]->l_old->v = tree->a_edges[i]->l->v;             */
206          Record_Br_Len(tree);
207          Set_Both_Sides(YES,tree);
208          Lk(NULL,tree);
209        }
210
211      if(step > n_step_max) break;
212
213      if(tree->io->print_trace)
214        {
215          char *s = Write_Tree(tree,NO);
216          PhyML_Fprintf(tree->io->fp_out_trace,"[%f]%s\n",tree->c_lnL,s); fflush(tree->io->fp_out_trace);
217          if(tree->io->print_site_lnl) Print_Site_Lk(tree,tree->io->fp_out_lk); fflush(tree->io->fp_out_lk);
218          Free(s);
219        }
220
221      if((tree->mod->s_opt->print) && (!tree->io->quiet)) Print_Lk(tree,"[Topology           ]");
222     
223
224/*       if(((tree->c_lnL > old_loglk) && (FABS(old_loglk-tree->c_lnL) < tree->mod->s_opt->min_diff_lk_local)) || (n_without_swap > it_lim_without_swap)) break; */
225      if((FABS(old_loglk-tree->c_lnL) < tree->mod->s_opt->min_diff_lk_global) || (n_without_swap > it_lim_without_swap)) break;
226     
227      Fill_Dir_Table(tree);
228      Fix_All(tree);
229      n_neg = 0;
230      For(i,2*tree->n_otu-3)
231        if((!tree->a_edges[i]->left->tax) && 
232           (!tree->a_edges[i]->rght->tax))
233          NNI(tree,tree->a_edges[i],0);
234     
235     
236
237      Select_Edges_To_Swap(tree,sorted_b,&n_neg);                 
238      Sort_Edges_NNI_Score(tree,sorted_b,n_neg);
239      Optimiz_Ext_Br(tree);                 
240      Update_Bl(tree,lambda);
241                       
242      n_tested = 0;
243      For(i,(int)CEIL((phydbl)n_neg*(lambda)))
244        tested_b[n_tested++] = sorted_b[i];
245     
246
247      Make_N_Swap(tree,tested_b,0,n_tested);
248
249      if((tree->mod->s_opt->print) && (!tree->io->quiet)) PhyML_Printf("[# nnis=%3d]",n_tested);
250
251      n_tot_swap += n_tested;
252     
253      if(n_tested > 0) n_without_swap = 0;
254      else             n_without_swap++;
255     
256      n_iter+=1.0;
257    }
258  while(1);
259   
260/*   Round_Optimize(tree,tree->data); */
261
262  Free(sorted_b);
263  Free(tested_b);
264
265  return n_tot_swap;
266}
267
268//////////////////////////////////////////////////////////////
269//////////////////////////////////////////////////////////////
270
271
272void Simu_Pars(t_tree *tree, int n_step_max)
273{
274  phydbl old_pars,n_iter,lambda;
275  int i,n_neg,n_tested,n_without_swap,n_tot_swap,step;
276  t_edge **sorted_b,**tested_b;
277  int each;
278 
279  sorted_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
280  tested_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
281 
282  old_pars            = 0;
283  tree->c_pars        = 0;
284  n_iter              = 1.0;
285  n_tested            = 0;
286  n_without_swap      = 0;
287  step                = 0;
288  each                = 4;
289  lambda              = .75;
290  n_tot_swap          = 0;
291 
292  Update_Dirs(tree);
293 
294  if((tree->mod->s_opt->print) && (!tree->io->quiet)) PhyML_Printf("\n\n. Starting simultaneous NNI moves (parsimony criterion)...\n");
295 
296  do
297    {
298      ++step;
299     
300      if(step > n_step_max) break;
301     
302      each--;
303     
304      Set_Both_Sides(YES,tree);
305      Pars(NULL,tree);
306     
307      if((tree->mod->s_opt->print) && (!tree->io->quiet)) 
308        {
309          Print_Pars(tree);
310          if(step > 1) (n_tested > 1)?(printf("[%4d NNIs]",n_tested)):(printf("[%4d NNI ]",n_tested));
311        }
312     
313
314      if(FABS(old_pars - tree->c_pars) < SMALL) break;
315             
316      if((tree->c_pars > old_pars) && (step > 1))
317        {
318          if((tree->mod->s_opt->print) && (!tree->io->quiet))
319            PhyML_Printf("\n\n. Moving backward (topoLlogy) \n");
320          if(!Mov_Backward_Topo_Pars(tree,old_pars,tested_b,n_tested))
321            Exit("\n. Err: mov_back failed\n");
322          if(!tree->n_swap) n_neg = 0;
323         
324         
325          Set_Both_Sides(YES,tree);
326          Pars(NULL,tree);
327        }
328      else 
329        {
330         
331          old_pars = tree->c_pars;         
332          Fill_Dir_Table(tree);
333
334          n_neg = 0;
335          For(i,2*tree->n_otu-3)
336            if((!tree->a_edges[i]->left->tax) && 
337               (!tree->a_edges[i]->rght->tax)) 
338              NNI_Pars(tree,tree->a_edges[i],NO);
339         
340          Select_Edges_To_Swap(tree,sorted_b,&n_neg);             
341          Sort_Edges_NNI_Score(tree,sorted_b,n_neg);       
342         
343          n_tested = 0;
344          For(i,(int)CEIL((phydbl)n_neg*(lambda)))
345            tested_b[n_tested++] = sorted_b[i];
346         
347          Make_N_Swap(tree,tested_b,0,n_tested);
348         
349          n_tot_swap += n_tested;
350         
351          if(n_tested > 0) n_without_swap = 0;
352          else             n_without_swap++;
353        }
354      n_iter+=1.0;
355    }
356  while(1);
357 
358  Free(sorted_b);
359  Free(tested_b);
360}
361
362//////////////////////////////////////////////////////////////
363//////////////////////////////////////////////////////////////
364
365
366void Select_Edges_To_Swap(t_tree *tree, t_edge **sorted_b, int *n_neg)
367{
368  int i;
369  t_edge *b;
370  phydbl best_score;
371
372  *n_neg = 0;
373
374  For(i,2*tree->n_otu-3)
375    {
376      b = tree->a_edges[i];
377      best_score = b->nni->score;
378
379      if((!b->left->tax) && (!b->rght->tax) && (b->nni->score < -tree->mod->s_opt->min_diff_lk_move)) 
380        {
381          Check_NNI_Scores_Around(b->left,b->rght,b,&best_score,tree);
382          Check_NNI_Scores_Around(b->rght,b->left,b,&best_score,tree);
383          if(best_score < b->nni->score) continue;
384          sorted_b[*n_neg] = b;
385          (*n_neg)++;
386        }
387    }
388}
389
390//////////////////////////////////////////////////////////////
391//////////////////////////////////////////////////////////////
392
393
394void Update_Bl(t_tree *tree, phydbl fact)
395{
396  int i;
397  t_edge *b,*orig;
398
399  For(i,2*tree->n_otu-3)
400    {
401      b = tree->a_edges[i];
402      /* b->l->v = b->l_old->v + (b->nni->l0 - b->l_old->v)*fact;       */
403     
404      orig = b;
405      do
406        {
407          b->l->v = b->l_old->v + (b->nni->l0 - b->l_old->v)*fact;     
408          if(b->next) b = b->next;
409          else         b = b->next;
410        }
411      while(b);
412      b = orig;
413
414    }
415}
416
417//////////////////////////////////////////////////////////////
418//////////////////////////////////////////////////////////////
419
420
421void Make_N_Swap(t_tree *tree,t_edge **b, int beg, int end)
422{
423  int i;
424  int dim;
425  t_edge *orig;
426
427  dim = 2*tree->n_otu-2;
428
429  /* PhyML_Printf("\n. Beg Actually performing swaps\n"); */
430  tree->n_swap = 0;
431  for(i=beg;i<end;i++)
432    {
433      /* we use t_dir here to take into account previous modifications of the topology */
434      /* printf("\n. Swap on edge %d [%d %d %d %d]",b[i]->num, */
435      /*             b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]]->num, */
436      /*             b[i]->nni->swap_node_v2->num, */
437      /*             b[i]->nni->swap_node_v3->num, */
438      /*             b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]]->num); */
439
440      Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
441           b[i]->nni->swap_node_v2,
442           b[i]->nni->swap_node_v3,
443           b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
444           tree);
445
446      if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
447        {
448          /* Undo this swap as it violates one of the topological constraints
449             defined in the input constraint tree
450          */
451          Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
452               b[i]->nni->swap_node_v2,
453               b[i]->nni->swap_node_v3,
454               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
455               tree);   
456        }
457
458      if(tree->n_root)
459        {
460          tree->n_root->v[2] = tree->e_root->left;
461          tree->n_root->v[1] = tree->e_root->rght;
462        }
463
464      orig = b[i];
465      do
466        {
467          b[i]->l->v = b[i]->nni->best_l;
468          if(b[i]->next) b[i] = b[i]->next;
469          else            b[i] = b[i]->next;
470        }
471      while(b[i]);
472      b[i] = orig;
473
474      tree->n_swap++;
475    }
476
477
478  /* PhyML_Printf("\n. End Actually performing swaps\n"); */
479
480}
481
482//////////////////////////////////////////////////////////////
483//////////////////////////////////////////////////////////////
484
485
486int Make_Best_Swap(t_tree *tree)
487{
488  int i,j,return_value;
489  t_edge *b,**sorted_b;
490  int dim;
491  t_edge *orig;
492
493  dim = 2*tree->n_otu-2;
494
495  sorted_b = (t_edge **)mCalloc(tree->n_otu-3,sizeof(t_edge *));
496 
497  j=0;
498  For(i,2*tree->n_otu-3) if((!tree->a_edges[i]->left->tax) &&
499                            (!tree->a_edges[i]->rght->tax))
500                              sorted_b[j++] = tree->a_edges[i];
501
502  Sort_Edges_NNI_Score(tree,sorted_b,tree->n_otu-3);
503
504  if(sorted_b[0]->nni->score < -0.0)
505    {
506      b = sorted_b[0];
507      return_value = 1;
508
509      Swap(b->nni->swap_node_v2->v[tree->t_dir[b->nni->swap_node_v2->num*dim+b->nni->swap_node_v1->num]],
510           b->nni->swap_node_v2,
511           b->nni->swap_node_v3,
512           b->nni->swap_node_v3->v[tree->t_dir[b->nni->swap_node_v3->num*dim+b->nni->swap_node_v4->num]],
513           tree);
514
515      if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
516        {
517          /* Undo this swap as it violates one of the topological constraints
518             defined in the input constraint tree
519          */
520          Swap(b->nni->swap_node_v2->v[tree->t_dir[b->nni->swap_node_v2->num*dim+b->nni->swap_node_v1->num]],
521               b->nni->swap_node_v2,
522               b->nni->swap_node_v3,
523               b->nni->swap_node_v3->v[tree->t_dir[b->nni->swap_node_v3->num*dim+b->nni->swap_node_v4->num]],
524               tree);   
525        }
526
527      /* b->l->v = b->nni->best_l; */
528
529      orig = b;
530      do
531        {
532          b->l->v = b->nni->best_l;
533          if(b->next) b = b->next;
534          else         b = b->next;
535        }
536      while(b);
537      b = orig;
538
539
540/*       (b->nni->best_conf == 1)? */
541/*      (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v1],tree)): */
542/*      (Swap(b->left->v[b->l_v2],b->left,b->rght,b->rght->v[b->r_v2],tree)); */
543     
544/*       b->l->v =  */
545/*      (b->nni->best_conf == 1)? */
546/*      (b->nni->l1): */
547/*      (b->nni->l2); */
548
549
550    }
551  else return_value = 0;
552
553  Free(sorted_b);
554
555  return return_value;
556}
557
558//////////////////////////////////////////////////////////////
559//////////////////////////////////////////////////////////////
560
561
562int Mov_Backward_Topo_Bl(t_tree *tree, phydbl lk_old, t_edge **tested_b, int n_tested)
563{
564  phydbl **l_init;
565  int i,j,step,beg,end;
566  t_edge *b,*orig;
567
568  l_init = (phydbl **)mCalloc(2*tree->n_otu-3,sizeof(phydbl *));
569
570  For(i,2*tree->n_otu-3) l_init[i] = MIXT_Get_Lengths_Of_This_Edge(tree->a_edges[i],tree);
571 
572  step = 2;
573  do
574    {
575      For(i,2*tree->n_otu-3) 
576        {
577          b = tree->a_edges[i];
578
579          /* b->l->v = b->l_old->v + (1./step) * (l_init[i] - b->l_old->v); */
580
581          j = 0;
582          orig = b;
583          do
584            {
585              b->l->v = b->l_old->v + (1./step) * (l_init[i][j] - b->l_old->v);
586              if(b->next) b = b->next;
587              else         b = b->next;
588              j++;
589            }
590          while(b);
591          b = orig;
592        }
593
594      beg = (int)FLOOR((phydbl)n_tested/(step-1));
595      end = 0;
596      Unswap_N_Branch(tree,tested_b,beg,end);
597      beg = 0;
598      end = (int)FLOOR((phydbl)n_tested/step);
599      Swap_N_Branch(tree,tested_b,beg,end);
600     
601      if(!end) tree->n_swap = 0;
602           
603      Set_Both_Sides(NO,tree);
604      Lk(NULL,tree);
605     
606      step++;
607
608    }while((tree->c_lnL < lk_old) && (step < 1000));
609
610 
611  if(step == 1000)
612    {
613      if(tree->n_swap)  Exit("\n== Err. in Mov_Backward_Topo_Bl (n_swap > 0)\n");
614
615      For(i,2*tree->n_otu-3) 
616        {
617          b = tree->a_edges[i];
618
619          orig = b;
620          do
621            {
622              b->l->v = b->l_old->v;
623              if(b->next) b = b->next;
624              else         b = b->next;
625            }
626          while(b);
627          b = orig;
628        }
629     
630     
631      Set_Both_Sides(NO,tree);
632      Lk(NULL,tree);
633    }
634 
635  For(i,2*tree->n_otu-3) Free(l_init[i]);
636  Free(l_init);
637
638  tree->n_swap = 0;
639  For(i,2*tree->n_otu-3) 
640    {
641      if(tree->a_edges[i]->nni->score < 0.0) tree->n_swap++;
642      tree->a_edges[i]->nni->score = +1.0;
643    }
644
645 
646  if(tree->c_lnL > lk_old)                                return  1;
647  else if((tree->c_lnL > lk_old-tree->mod->s_opt->min_diff_lk_local) && 
648          (tree->c_lnL < lk_old+tree->mod->s_opt->min_diff_lk_local)) return -1;
649  else                                                    return  0;
650}
651
652//////////////////////////////////////////////////////////////
653//////////////////////////////////////////////////////////////
654
655
656int Mov_Backward_Topo_Pars(t_tree *tree, int pars_old, t_edge **tested_b, int n_tested)
657{
658  int i,step,beg,end;
659 
660  step = 2;
661  do
662    {
663      beg = (int)FLOOR((phydbl)n_tested/(step-1));
664      end = 0;
665      Unswap_N_Branch(tree,tested_b,beg,end);
666      beg = 0;
667      end = (int)FLOOR((phydbl)n_tested/step);
668      Swap_N_Branch(tree,tested_b,beg,end);
669     
670      if(!end) tree->n_swap = 0;
671     
672      Set_Both_Sides(NO,tree);     
673      Pars(NULL,tree);
674     
675      step++;
676
677    }while((tree->c_pars > pars_old) && (step < 1000));
678
679 
680  if(step == 1000)
681    {
682      if(tree->n_swap)  Exit("\n. Err. in Mov_Backward_Topo_Bl (n_swap > 0)\n");
683
684      Set_Both_Sides(NO,tree);
685      Pars(NULL,tree);
686    }
687
688  tree->n_swap = 0;
689  For(i,2*tree->n_otu-3) 
690    {
691      if(tree->a_edges[i]->nni->score < 0.0) tree->n_swap++;
692      tree->a_edges[i]->nni->score = +1.0;
693    }
694
695 
696  if(tree->c_pars < pars_old)       return  1;
697  else if(tree->c_pars == pars_old) return -1;
698  else                              return  0;
699}
700
701//////////////////////////////////////////////////////////////
702//////////////////////////////////////////////////////////////
703
704
705void Unswap_N_Branch(t_tree *tree, t_edge **b, int beg, int end)
706{
707  int i;
708  int dim;
709  t_edge *orig;
710
711  dim = 2*tree->n_otu-2;
712
713  if(end>beg)
714    {
715      for(i=beg;i<end;i++)
716        {
717         
718/*        PhyML_Printf("MOV BACK UNSWAP Edge %d Swap nodes %d(%d) %d %d %d(%d)\n", */
719/*               b[i]->num, */
720/*               b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num][b[i]->nni->swap_node_v1->num]]->num, */
721/*               b[i]->nni->swap_node_v4->num, */
722/*               b[i]->nni->swap_node_v2->num, */
723/*               b[i]->nni->swap_node_v3->num, */
724/*               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num][b[i]->nni->swap_node_v4->num]]->num, */
725/*               b[i]->nni->swap_node_v1->num */
726/*               ); */
727
728          Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
729               b[i]->nni->swap_node_v2,
730               b[i]->nni->swap_node_v3,
731               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
732               tree);
733
734          if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
735            {
736              /* Undo this swap as it violates one of the topological constraints
737                 defined in the input constraint tree
738              */
739              Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
740                   b[i]->nni->swap_node_v2,
741                   b[i]->nni->swap_node_v3,
742                   b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
743                   tree);       
744            }
745
746
747/*        (b[i]->nni->best_conf == 1)? */
748/*          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v1],tree)): */
749/*          (Swap(b[i]->left->v[b[i]->l_v2],b[i]->left,b[i]->rght,b[i]->rght->v[b[i]->r_v2],tree)); */
750
751          /* b[i]->l->v = b[i]->l_old->v; */
752
753
754          orig = b[i];
755          do
756            {
757              b[i]->l->v = b[i]->l_old->v;
758              if(b[i]->next) b[i] = b[i]->next;
759              else            b[i] = b[i]->next;
760            }
761          while(b[i]);
762          b[i] = orig;
763
764        }
765    }
766  else
767    {
768      for(i=beg-1;i>=end;i--)
769        {       
770          Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
771               b[i]->nni->swap_node_v2,
772               b[i]->nni->swap_node_v3,
773               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
774               tree);
775
776          if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
777            {
778              /* Undo this swap as it violates one of the topological constraints
779                 defined in the input constraint tree
780              */
781              Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
782                   b[i]->nni->swap_node_v2,
783                   b[i]->nni->swap_node_v3,
784                   b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
785                   tree);       
786            }
787         
788
789          /* b[i]->l->v = b[i]->l_old->v; */
790
791          orig = b[i];
792          do
793            {
794              b[i]->l->v = b[i]->l_old->v;
795              if(b[i]->next) b[i] = b[i]->next;
796              else            b[i] = b[i]->next;
797            }
798          while(b[i]);
799          b[i] = orig;
800
801        }
802    }
803}
804
805//////////////////////////////////////////////////////////////
806//////////////////////////////////////////////////////////////
807
808
809void Swap_N_Branch(t_tree *tree,t_edge **b, int beg, int end)
810{
811  int i;
812  int dim;
813  t_edge *orig;
814
815  dim = 2*tree->n_otu-2;
816
817  if(end>beg)
818    {
819      for(i=beg;i<end;i++)
820        {
821          Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
822               b[i]->nni->swap_node_v2,
823               b[i]->nni->swap_node_v3,
824               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
825               tree);
826
827          if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
828            {
829              /* Undo this swap as it violates one of the topological constraints
830                 defined in the input constraint tree
831              */
832              Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
833                   b[i]->nni->swap_node_v2,
834                   b[i]->nni->swap_node_v3,
835                   b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
836                   tree);       
837            }
838
839          /* b[i]->l->v = b[i]->nni->best_l; */
840
841          orig = b[i];
842          do
843            {
844              b[i]->l->v = b[i]->nni->best_l;
845              if(b[i]->next) b[i] = b[i]->next;
846              else            b[i] = b[i]->next;
847            }
848          while(b[i]);
849          b[i] = orig;
850
851        }
852    }
853  else
854    {
855      for(i=beg-1;i>=end;i--)
856        {
857          Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
858               b[i]->nni->swap_node_v2,
859               b[i]->nni->swap_node_v3,
860               b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
861               tree);
862
863          if(!Check_Topo_Constraints(tree,tree->io->cstr_tree))
864            {
865              /* Undo this swap as it violates one of the topological constraints
866                 defined in the input constraint tree
867              */
868              Swap(b[i]->nni->swap_node_v2->v[tree->t_dir[b[i]->nni->swap_node_v2->num*dim+b[i]->nni->swap_node_v1->num]],
869                   b[i]->nni->swap_node_v2,
870                   b[i]->nni->swap_node_v3,
871                   b[i]->nni->swap_node_v3->v[tree->t_dir[b[i]->nni->swap_node_v3->num*dim+b[i]->nni->swap_node_v4->num]],
872                   tree);       
873            }
874
875          /* b[i]->l->v = b[i]->nni->best_l; */
876
877          orig = b[i];
878          do
879            {
880              b[i]->l->v = b[i]->nni->best_l;
881              if(b[i]->next) b[i] = b[i]->next;
882              else            b[i] = b[i]->next;
883            }
884          while(b[i]);
885          b[i] = orig;
886
887        }
888    }
889}
890
891//////////////////////////////////////////////////////////////
892//////////////////////////////////////////////////////////////
893
894
895void Check_NNI_Scores_Around(t_node *a, t_node *d, t_edge *b, phydbl *best_score, t_tree *tree)
896{
897
898  int i;
899  For(i,3)
900    {
901      if((d->v[i] != a) && (!d->v[i]->tax))
902        {
903          if((d->b[i]->nni->score > *best_score-1.E-10) &&
904             (d->b[i]->nni->score < *best_score+1.E-10)) /* ties */
905             {
906               d->b[i]->nni->score = *best_score+1.;
907             }
908
909          if(d->b[i]->nni->score < *best_score)
910            {
911              *best_score = d->b[i]->nni->score;
912            }
913        }
914    }
915}
916
917//////////////////////////////////////////////////////////////
918//////////////////////////////////////////////////////////////
919
920//////////////////////////////////////////////////////////////
921//////////////////////////////////////////////////////////////
922
923//////////////////////////////////////////////////////////////
924//////////////////////////////////////////////////////////////
925
926//////////////////////////////////////////////////////////////
927//////////////////////////////////////////////////////////////
928
929//////////////////////////////////////////////////////////////
930//////////////////////////////////////////////////////////////
931
932//////////////////////////////////////////////////////////////
933//////////////////////////////////////////////////////////////
934
935//////////////////////////////////////////////////////////////
936//////////////////////////////////////////////////////////////
937
938//////////////////////////////////////////////////////////////
939//////////////////////////////////////////////////////////////
940
941//////////////////////////////////////////////////////////////
942//////////////////////////////////////////////////////////////
943
944//////////////////////////////////////////////////////////////
945//////////////////////////////////////////////////////////////
946
947//////////////////////////////////////////////////////////////
948//////////////////////////////////////////////////////////////
949
Note: See TracBrowser for help on using the repository browser.