source: trunk/GDE/PHYML20130708/phyml/src/mg.c

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

added most recent version of phyml
File size: 78.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 <config.h>
14
15#include "mg.h"
16#include "free.h"
17#include "help.h"
18#include "utilities.h"
19#include "optimiz.h"
20#include "models.h"
21#include "simu.h"
22#include "lk.h"
23#include "pars.h"
24#include "interface.h"
25
26//////////////////////////////////////////////////////////////
27//////////////////////////////////////////////////////////////
28
29
30int PART_main(int argc, char **argv)
31{
32  option *io;
33  char *s_tree;
34  FILE *fp_phyml_tree,*fp_phyml_stats,*fp_phyml_lk;
35  int part;
36  time_t t_beg,t_end;
37  div_t hour,min;
38  int r_seed;
39  int i;
40
41  fflush(NULL);
42  io = (option *)Get_Input(argc,argv);
43  r_seed = (io->r_seed < 0)?(time(NULL)):(io->r_seed);
44  srand(r_seed);
45  fp_phyml_stats = Openfile(io->out_stats_file,io->out_stats_file_open_mode);
46  PhyML_Fprintf(fp_phyml_stats,"\n- PHYML %s -\n\n", VERSION);
47  fp_phyml_tree = Openfile(io->out_tree_file,io->out_tree_file_open_mode);
48  fp_phyml_lk = fopen(io->out_lk_file,"w");
49
50  time(&t_beg);
51
52  if(io->multigene)
53    {
54      align ***data;
55      calign **cdata;
56      t_mod **mod;
57      matrix **mat;
58      t_treelist *treelist;
59      supert_tree *st;
60
61      data  = (align ***)  mCalloc(io->n_part,sizeof(align **));
62      cdata = (calign **)mCalloc(io->n_part,sizeof(calign *));
63      mod   = (t_mod **) mCalloc(io->n_part,sizeof(t_mod *));
64      mat   = (matrix **)mCalloc(io->n_part,sizeof(matrix *));
65
66      /* Read the sequences (for each partition) */
67      For(part,io->n_part)
68        {
69          Make_Model_Complete(io->st->optionlist[part]->mod); /* Complete model for each data part */
70          data[part] = Get_Seq(io->st->optionlist[part]);
71          Make_Model_Complete(io->st->optionlist[part]->mod);
72          Set_Model_Name(io->st->optionlist[part]->mod);
73          mod[part]  = io->st->optionlist[part]->mod;
74
75          PhyML_Printf("\n. Data part [#%d]\n",part+1);
76          PhyML_Printf("\n. Compressing sequences...\n");
77          cdata[part] = Compact_Data(data[part],io->st->optionlist[part]);
78          fclose(io->st->optionlist[part]->fp_in_align);
79          Free_Seq(data[part],cdata[part]->n_otu);
80          Init_Model(cdata[part],mod[part],io->st->optionlist[part]);
81          Check_Ambiguities(cdata[part],
82                            io->st->optionlist[part]->mod->io->datatype,
83                            io->st->optionlist[part]->mod->io->state_len);
84        }
85
86      PART_Make_Supert_tree_Full(io->st,io,cdata);
87      st = io->st;
88      treelist = st->treelist;
89      Fill_Dir_Table(st->tree);
90      Update_Dirs(st->tree);
91
92      For(part,io->n_part)
93        {
94          st->curr_cdata = cdata[part];
95          if(!PART_Get_Species_Found_In_St(st,cdata[part])) break;
96          treelist->tree[part] = Make_Tree_From_Scratch(st->tree->n_otu,NULL);
97          Copy_Tree_Topology_With_Labels(st->tree,treelist->tree[part]);
98          treelist->tree[part]->num_curr_branch_available = 0;
99          Connect_Edges_To_Nodes_Recur(treelist->tree[part]->a_nodes[0],
100                                       treelist->tree[part]->a_nodes[0]->v[0],
101                                       treelist->tree[part]);
102          PART_Prune_St_Topo(treelist->tree[part],cdata[part],st);
103
104          if(treelist->tree[part]->n_otu != cdata[part]->n_otu)
105            {
106              PhyML_Printf("\n. Problem with sequence file '%s'\n",io->st->optionlist[part]->in_align_file);
107              PhyML_Printf("\n. # taxa found in supertree restricted to '%s' taxa = %d\n",
108                     io->st->optionlist[part]->in_align_file,
109                     treelist->tree[part]->n_otu);
110              PhyML_Printf("\n. # sequences in '%s' = %d\n",
111                     io->st->optionlist[part]->in_align_file,
112                     cdata[part]->n_otu);
113              Exit("\n");
114            }
115
116          treelist->tree[part]->dp         = part;
117          treelist->tree[part]->n_otu      = cdata[part]->n_otu;
118          treelist->tree[part]->mod        = mod[part];
119          treelist->tree[part]->io         = io->st->optionlist[part];
120          treelist->tree[part]->data       = cdata[part];
121          treelist->tree[part]->n_pattern  = treelist->tree[part]->data->crunch_len/
122                                             treelist->tree[part]->io->state_len;
123
124          Set_Both_Sides(YES,treelist->tree[part]);
125          Connect_CSeqs_To_Nodes(cdata[part],treelist->tree[part]);
126          Fill_Dir_Table(treelist->tree[part]);
127          Update_Dirs(treelist->tree[part]);
128          Make_Tree_4_Lk(treelist->tree[part],cdata[part],cdata[part]->init_len);
129          Make_Tree_4_Pars(treelist->tree[part],cdata[part],cdata[part]->init_len);
130          treelist->tree[part]->triplet_struct = Make_Triplet_Struct(treelist->tree[part]->mod);
131          Init_Triplet_Struct(treelist->tree[part]->triplet_struct);
132        }
133
134      if(part != io->n_part)
135        {
136          PhyML_Printf("\n. Sequence data part found in '%s' has one or more taxa not found in the '%s' tree file\n",
137                 io->st->optionlist[part]->in_align_file,
138                 io->in_tree_file);
139          Exit("\n");
140        }
141
142      PART_Check_Extra_Taxa(st);
143       
144      st->tree->c_lnL = .0;
145      For(part,io->n_part)
146        {
147          Lk(NULL,treelist->tree[part]);
148/*        Get_List_Of_Reachable_Tips(treelist->tree[part]); */
149          PART_Match_St_Nodes_In_Gt(treelist->tree[part],st);
150          PART_Match_St_Edges_In_Gt(treelist->tree[part],st);
151          PART_Map_St_Nodes_In_Gt(treelist->tree[part],st);
152          PART_Map_St_Edges_In_Gt(treelist->tree[part],st);
153          PART_Map_Gt_Edges_In_St(treelist->tree[part],st);
154          st->tree->c_lnL += treelist->tree[part]->c_lnL;
155        }
156
157      PART_Initialise_Bl_Partition(st);
158      PART_Set_Bl(st->bl,st);
159
160      time(&(st->tree->t_beg));
161      time(&(st->tree->t_current));
162      PhyML_Printf("\n. (%5d sec) [00] [%10.2f] [%5d]\n",
163             (int)(st->tree->t_current-st->tree->t_beg),
164             PART_Lk(st),
165             PART_Pars(st));
166
167
168      int n_iter=0;
169      do
170        {
171          PART_Optimize_Br_Len_Serie(st->tree->a_nodes[0],
172                                   st->tree->a_nodes[0]->v[0],
173                                   st->tree->a_nodes[0]->b[0],
174                                   st);
175
176          Set_Both_Sides(YES,st->tree);
177          PART_Lk(st);
178          PhyML_Printf("\n. %f",st->tree->c_lnL);
179/*        For(part,st->n_part) PhyML_Printf("\n. %s",Write_Tree(st->treelist->tree[part],NO)); */
180          n_iter++;
181        }while(n_iter < 5);
182/*       Exit("\n"); */
183
184
185/*       PART_Lk(st); */
186/*       PhyML_Printf("\n> %f",st->tree->c_lnL); */
187/*       For(i,2*st->tree->n_otu-3) */
188/*      { */
189/*        if((!st->tree->a_edges[i]->left->tax) && (!st->tree->a_edges[i]->rght->tax)) */
190/*          { */
191/*            PART_NNI(st->tree->a_edges[i],st); */
192/*          } */
193/*      } */
194
195     
196/*       if(io->mod->s_opt->topo_search == NNI_MOVE) */
197      PART_Simu(st);
198/*       else */
199/*       PART_Speed_Spr(st); */
200
201
202      time(&t_end);
203
204      hour = div(t_end-t_beg,3600);
205      min  = div(t_end-t_beg,60  );
206
207      min.quot -= hour.quot*60;
208
209
210      PhyML_Fprintf(fp_phyml_stats,"\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n");
211      PhyML_Fprintf(fp_phyml_stats,"\n. Number of partitions = %d\n\n",st->n_part);
212      PhyML_Fprintf(fp_phyml_stats,"\n. Full data set -- lnL = %f\n\n",st->tree->c_lnL);
213      PhyML_Fprintf(fp_phyml_stats,"\n. Tree search algorithm : %s\n\n",(io->mod->s_opt->topo_search == NNI_MOVE)?("NNIs"):("SPRs"));
214      PhyML_Fprintf(fp_phyml_stats,"\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n\n");
215      For(part,io->n_part)
216        {
217          Print_Fp_Out(fp_phyml_stats,t_beg,t_end,st->treelist->tree[part],
218                       io->st->optionlist[part],1,1,YES);
219        }
220
221
222      PhyML_Printf("\n\n. Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
223      PhyML_Printf("\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n");
224
225      For(i,2*st->tree->n_otu-3) st->tree->a_edges[i]->l->v = 0.1;
226      s_tree = Write_Tree(st->tree,NO);
227      PhyML_Fprintf(fp_phyml_tree,"Supertree\n");
228      PhyML_Fprintf(fp_phyml_tree,"%s\n",s_tree);
229      Free(s_tree);
230      For(part,st->n_part)
231        {
232          PhyML_Fprintf(fp_phyml_tree,"Gene tree number %d\n",part+1);
233          s_tree = Write_Tree(st->treelist->tree[part],NO);
234          PhyML_Fprintf(fp_phyml_tree,"%s\n",s_tree);
235          Free(s_tree);
236        }
237
238
239      For(part,st->n_part)
240        {
241          if(io->mod->s_opt->topo_search == SPR_MOVE) Free_Spr_List(treelist->tree[part]);
242          Free_Tree_Lk(treelist->tree[part]);
243          Free_Tree_Pars(treelist->tree[part]);
244          Free_Triplet(treelist->tree[part]->triplet_struct);
245          Free_Tree(treelist->tree[part]);
246          Free_Cseq(cdata[part]);
247          Free_Model(mod[part]);
248          Free_Input(io->st->optionlist[part]);
249
250        }
251
252      if(io->mod->s_opt->topo_search == SPR_MOVE) Free_Spr_List(st->tree);
253      Free(mat);
254      Free(mod);
255      Free(data);
256      Free(cdata);
257      Free(treelist->tree);
258      Free(treelist);
259      Free_St(st);
260    }
261
262
263  if(io->fp_in_align ) fclose(io->fp_in_align);
264  if(io->fp_in_tree) fclose(io->fp_in_tree);
265
266
267  Free_Model(io->mod);
268  Free_Input(io);
269
270  fclose(fp_phyml_lk);
271  fclose(fp_phyml_tree);
272  fclose(fp_phyml_stats);
273
274
275  return 0;
276}
277
278//////////////////////////////////////////////////////////////
279//////////////////////////////////////////////////////////////
280
281
282void PART_Print_Nodes(t_node *a, t_node *d, supert_tree *st)
283{
284  int i;
285  PhyML_Printf(">>>>>>>>>>>>>>>>>>>>\n");
286  PhyML_Printf("num t_node = %d\n",d->num);
287  if(d->tax) PhyML_Printf("name='%s'\n",d->name);
288  else
289    {
290      PhyML_Printf("n_of_reachable_tips : \n");
291      For(i,3)
292        {
293/*        PhyML_Printf("dir%d=%d; ",i,st->n_of_reachable_tips[st->num_data_of_interest][d->num][i]); */
294/*        For(j,st->n_of_reachable_tips[st->num_data_of_interest][d->num][i]) */
295/*          { */
296/*            PhyML_Printf("%s ", */
297/*                   st->list_of_reachable_tips[st->num_data_of_interest][d->num][i][j]->name); */
298/*          } */
299         
300          PhyML_Printf("\n");
301        }
302    }
303  PhyML_Printf("<<<<<<<<<<<<<<<<<<<\n\n");
304  if(d->tax) return;
305  else
306    {
307      For(i,3) if(d->v[i] != a) PART_Print_Nodes(d,d->v[i],st);
308    }
309}
310
311//////////////////////////////////////////////////////////////
312//////////////////////////////////////////////////////////////
313
314
315supert_tree *PART_Make_Supert_tree_Light(option *input)
316{
317  supert_tree *st;
318
319  st = (supert_tree *)mCalloc(1,sizeof(supert_tree));
320  st->optionlist = (option **)mCalloc(input->n_part,sizeof(option *));
321  st->bl_partition = (int *)mCalloc(input->n_part,sizeof(int ));
322  return st;
323}
324
325//////////////////////////////////////////////////////////////
326//////////////////////////////////////////////////////////////
327
328
329void PART_Make_Supert_tree_Full(supert_tree *st, option *io, calign **data)
330{
331  int i,j,k;
332
333  if(io->in_tree == 2)
334    {
335      PhyML_Printf("\n. Reading user tree...\n");
336      rewind(io->fp_in_tree);
337     
338      st->tree = Read_Tree_File_Phylip(io->fp_in_tree);
339     
340      if(!st->tree->has_branch_lengths)
341        {
342          PhyML_Printf("\n. Branch lengths are all set to 0.1...\n");
343          For(i,2*st->tree->n_otu-3) st->tree->a_edges[i]->l->v = 0.1;
344        }
345    }
346  else
347    {
348      Warn_And_Exit("\n. A user-defined input tree is needed\n");
349    }
350
351  st->tree->io      = io;
352  st->treelist      = (t_treelist *)Make_Treelist(io->n_part);
353  st->n_part          = io->n_part;
354  st->tree->mod     = io->mod;
355  st->lock_br_len   = 0;
356
357  st->map_st_node_in_gt = (t_node *****)mCalloc(st->n_part,sizeof(t_node ****));
358  For(i,st->n_part) 
359    {
360      st->map_st_node_in_gt[i] = (t_node ****)mCalloc(2*st->tree->n_otu-2,sizeof(t_node ***));
361      For(j,2*st->tree->n_otu-2) 
362        {
363          st->map_st_node_in_gt[i][j] = (t_node ***)mCalloc(3,sizeof(t_node **));
364          For(k,3) st->map_st_node_in_gt[i][j][k] = (t_node **)mCalloc(2,sizeof(t_node *));
365        }
366    }
367
368  st->map_st_edge_in_gt = (t_edge ***)mCalloc(st->n_part,sizeof(t_edge **));
369  For(i,st->n_part) st->map_st_edge_in_gt[i] = (t_edge **)mCalloc(2*st->tree->n_otu-3,sizeof(t_edge *));
370
371  st->map_gt_edge_in_st = (t_edge ****)mCalloc(st->n_part,sizeof(t_edge ***));
372  For(i,st->n_part)
373    {
374      st->map_gt_edge_in_st[i] = (t_edge ***)mCalloc(2*st->tree->n_otu-3,sizeof(t_edge **));
375      For(j,2*st->tree->n_otu-3) st->map_gt_edge_in_st[i][j] = (t_edge **)mCalloc(2*st->tree->n_otu-3,sizeof(t_edge *));
376    }
377
378  st->size_map_gt_edge_in_st = (int **)mCalloc(st->n_part,sizeof(int *));
379  For(i,st->n_part) st->size_map_gt_edge_in_st[i] = (int *)mCalloc(2*st->tree->n_otu-3,sizeof(int));
380
381
382  st->match_st_edge_in_gt = (t_edge ***)mCalloc(st->n_part,sizeof(t_edge **));
383  For(i,st->n_part) st->match_st_edge_in_gt[i] = (t_edge **)mCalloc(2*st->tree->n_otu-3,sizeof(t_edge *));
384
385  st->match_gt_edge_in_st = (t_edge ***)mCalloc(st->n_part,sizeof(t_edge **));
386  For(i,st->n_part) st->match_gt_edge_in_st[i] = (t_edge **)mCalloc(2*st->tree->n_otu-3,sizeof(t_edge *));
387
388  st->bl = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
389  For(i,st->n_part) st->bl[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl));
390
391  st->bl_cpy = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
392  For(i,st->n_part) st->bl_cpy[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl));
393
394  st->bl0 = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
395  For(i,st->n_part) st->bl0[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl));
396
397  st->bl1 = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
398  For(i,st->n_part) st->bl1[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl));
399
400  st->bl2 = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
401  For(i,st->n_part) st->bl2[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl));
402
403  st->s_mod = (t_mod **)mCalloc(st->n_part,sizeof(t_mod *));
404
405  For(i,2*st->tree->n_otu-3) Make_Edge_NNI(st->tree->a_edges[i]);
406
407  st->match_st_node_in_gt = (t_node ***)mCalloc(io->n_part,sizeof(t_node **));
408  For(i,io->n_part) st->match_st_node_in_gt[i] = (t_node **)mCalloc(2*st->tree->n_otu-2,sizeof(t_node *));
409}
410
411//////////////////////////////////////////////////////////////
412//////////////////////////////////////////////////////////////
413
414
415void PART_Prune_St_Topo(t_tree *tree, calign *data, supert_tree *st)
416{
417  int i,j,not_found;
418  int curr_ext_node, curr_int_node, curr_br, n_pruned_nodes;;
419  t_node **pruned_nodes;
420  t_edge **residual_edges;
421
422  pruned_nodes   = (t_node **)mCalloc(st->tree->n_otu,sizeof(t_node *));
423  residual_edges = (t_edge **)mCalloc(st->tree->n_otu,sizeof(t_edge *));
424
425  n_pruned_nodes = 0;
426  For(i,st->tree->n_otu)
427    {
428      For(j,data->n_otu)
429        {
430          if(!strcmp(data->c_seq[j]->name,st->tree->a_nodes[i]->name))
431            break;
432        }
433
434      not_found = 1;
435      if(j == data->n_otu)
436        {
437          For(j,tree->n_otu)
438            {
439              if(!strcmp(tree->a_nodes[j]->name,st->tree->a_nodes[i]->name))
440                {
441                  Prune_Subtree(tree->a_nodes[j]->v[0],
442                                tree->a_nodes[j],
443                                NULL,&(residual_edges[n_pruned_nodes]),
444                                tree);
445
446                  pruned_nodes[n_pruned_nodes] = tree->a_nodes[j];
447                  n_pruned_nodes++;
448                  not_found = 0;
449                  break;
450                }             
451            }
452
453
454          if(not_found)     
455            {
456              PhyML_Printf("\n. Taxon '%s'",st->tree->a_nodes[i]->name);
457              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
458              Warn_And_Exit("");
459            }
460        }
461    }
462
463  Free(tree->t_dir);
464
465  tree->n_otu -= n_pruned_nodes;
466
467  curr_ext_node = 0;
468  curr_int_node = tree->n_otu; 
469  curr_br = 0;
470  For(i,st->tree->n_otu)
471    {
472      For(j,n_pruned_nodes)
473        {
474          if(!strcmp(pruned_nodes[j]->name,st->tree->a_nodes[i]->name))
475            break;
476        }
477      if(j == n_pruned_nodes) /* That t_node still belongs to the tree */
478        {
479          Reassign_Node_Nums(tree->a_nodes[i],tree->a_nodes[i]->v[0], 
480                             &curr_ext_node, &curr_int_node,tree);
481          break;
482        }
483    }
484 
485  Reassign_Edge_Nums(tree->a_nodes[0],tree->a_nodes[0]->v[0],&curr_br,tree);
486
487  tree->t_dir = (short int *)mCalloc((2*tree->n_otu-2)*(2*tree->n_otu-2),sizeof(short int));
488
489  For(i,n_pruned_nodes) 
490    {
491      Free_Edge(residual_edges[i]);
492      Free_Edge(pruned_nodes[i]->b[0]);
493      Free_Node(pruned_nodes[i]->v[0]);
494      Free_Node(pruned_nodes[i]);
495    }
496
497  Free(pruned_nodes);
498  Free(residual_edges);
499
500}
501
502//////////////////////////////////////////////////////////////
503//////////////////////////////////////////////////////////////
504
505
506void PART_Match_St_Nodes_In_Gt(t_tree *gt, supert_tree *st)
507{
508  int i,j;
509
510  For(i,2*st->tree->n_otu-2) st->match_st_node_in_gt[gt->dp][i] = NULL; /* don't forget that step ! */
511
512  /* Map tips */
513  For(i,st->tree->n_otu)
514    {
515      For(j,gt->n_otu)
516        {
517          if(!strcmp(st->tree->a_nodes[i]->name,gt->a_nodes[j]->name))
518            {
519              st->match_st_node_in_gt[gt->dp][st->tree->a_nodes[i]->num] = gt->a_nodes[j];
520              break;
521            }
522        }
523    }
524
525#ifdef DEBUG
526  /* Checking that the results are correct so far */
527  int n_matches;
528  n_matches = 0;
529  For(i,2*st->tree->n_otu-2)
530    if(st->match_st_node_in_gt[gt->dp][i])
531      n_matches++;
532
533  if(n_matches != gt->n_otu)
534    {
535      PhyML_Printf("\n");
536      PhyML_Printf("\n. n_matches = %d 2*gt->n_otu-2 = %d\n",n_matches,2*gt->n_otu-2);
537      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
538      Warn_And_Exit("");
539    }
540#endif
541
542
543  /* Map internal nodes */
544  For(i,st->tree->n_otu)
545    {
546      if(st->match_st_node_in_gt[gt->dp][st->tree->a_nodes[i]->num])
547        {
548          PART_Match_St_Nodes_In_Gt_Recurr(st->match_st_node_in_gt[gt->dp][st->tree->a_nodes[i]->num],
549                                         st->match_st_node_in_gt[gt->dp][st->tree->a_nodes[i]->num]->v[0],
550                                         st->tree->a_nodes[i],
551                                         st->tree->a_nodes[i]->v[0],
552                                         gt,
553                                         st);
554          break;
555        }
556    }
557 
558
559
560#ifdef DEBUG
561  /* Checking that the results are correct */
562  n_matches = 0;
563  For(i,2*st->tree->n_otu-2) 
564    if(st->match_st_node_in_gt[gt->dp][st->tree->a_nodes[i]->num])
565        n_matches++;
566
567  if(n_matches != 2*gt->n_otu-2)
568    {
569      int j;
570      PhyML_Printf("\n");
571      PhyML_Printf("\n. n_matches = %d 2*gt->n_otu-2 = %d\n",n_matches,2*gt->n_otu-2);
572      For(j,2*gt->n_otu-2)
573        {
574          For(i,2*st->tree->n_otu-2) 
575            if(st->match_st_node_in_gt[gt->dp][i] == gt->a_nodes[j])
576              break;
577
578          if(i == 2*st->tree->n_otu-2)
579            {
580              PhyML_Printf("\n. Gt %3d t_node %3d (%3d %3d %3d) (%s %s %s) (%f %f %f) does not match\n",
581                     gt->dp,
582                     gt->a_nodes[j]->num,
583                     gt->a_nodes[j]->v[0] ? gt->a_nodes[j]->v[0]->num : -1,
584                     gt->a_nodes[j]->v[1] ? gt->a_nodes[j]->v[1]->num : -1,
585                     gt->a_nodes[j]->v[2] ? gt->a_nodes[j]->v[2]->num : -1,
586                     gt->a_nodes[j]->v[0]->tax ? gt->a_nodes[j]->v[0]->name : NULL,
587                     gt->a_nodes[j]->v[1]->tax ? gt->a_nodes[j]->v[1]->name : NULL,
588                     gt->a_nodes[j]->v[2]->tax ? gt->a_nodes[j]->v[2]->name : NULL,
589                     gt->a_nodes[j]->v[0] ? gt->a_nodes[j]->b[0]->l->v : -1.,
590                     gt->a_nodes[j]->v[1] ? gt->a_nodes[j]->b[1]->l->v : -1.,
591                     gt->a_nodes[j]->v[2] ? gt->a_nodes[j]->b[2]->l->v : -1.);
592            }
593        }
594
595      PhyML_Printf("oooooooo\n");
596      Print_Node(st->tree->a_nodes[0],
597                 st->tree->a_nodes[0]->v[0],
598                 st->tree);
599      PhyML_Printf(">>>>>>>\n");
600      For(i,st->n_part)
601        {
602          Print_Node(st->treelist->tree[i]->a_nodes[0],
603                     st->treelist->tree[i]->a_nodes[0]->v[0],
604                     st->treelist->tree[i]);
605          PhyML_Printf("<<<<<<<\n");
606        }
607      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
608      Warn_And_Exit("");
609    }
610#endif
611}
612
613//////////////////////////////////////////////////////////////
614//////////////////////////////////////////////////////////////
615
616
617void PART_Match_St_Nodes_In_Gt_Recurr(t_node *a_gt, t_node *d_gt, t_node *a_st, t_node *d_st, t_tree *gt, supert_tree *st)
618{
619  int i,j,k;
620  int *score_d_st;
621
622
623  if((d_gt->tax) || (d_st->tax)) return;
624  else
625    {
626      score_d_st = (int *)mCalloc(3,sizeof(int));
627     
628      /* Might be wrong. Check function Match_Nodes_In_Small_Tree */
629      For(i,3)
630        {
631          For(j,3)
632            {
633              For(k,d_st->bip_size[j])
634                {
635                  if(!strcmp(d_gt->bip_node[i][0]->name,d_st->bip_node[j][k]->name))
636                    {
637                      score_d_st[j] += 1;
638                      break;
639                    }
640                }
641            }
642        }
643
644
645      if((score_d_st[0] == 2) && (score_d_st[1] == 2) && (score_d_st[2] == 2))
646        {
647          st->match_st_node_in_gt[gt->dp][d_st->num] = d_gt;
648
649          For(i,3)
650            {
651              if(d_gt->v[i] != a_gt)
652                {
653                  For(j,3)
654                    {                 
655
656                      if(score_d_st[j] != 3)
657                        {
658                          PART_Match_St_Nodes_In_Gt_Recurr(d_gt,d_gt->v[i],d_st,d_st->v[j],gt,st);
659                          break;
660                        }
661
662/*                    For(k,d_st->n_of_reachable_tips[j]) */
663/*                      if(!strcmp(d_gt->list_of_reachable_tips[i][0]->name, */
664/*                                 d_st->list_of_reachable_tips[j][k]->name)) */
665/*                        { */
666/*                          PART_Match_St_Nodes_In_Gt_Recurr(d_gt,d_gt->v[i],d_st,d_st->v[j],gt,st); */
667/*                          break; */
668/*                        } */
669/*                    if(k != d_st->n_of_reachable_tips[j]) break; */
670                    }
671                }
672            }
673        }
674      else
675        {
676          For(i,3)
677            if(d_st->v[i] != a_st)
678              PART_Match_St_Nodes_In_Gt_Recurr(a_gt,d_gt,d_st,d_st->v[i],gt,st);
679        }
680      Free(score_d_st); 
681    }
682}
683
684//////////////////////////////////////////////////////////////
685//////////////////////////////////////////////////////////////
686
687
688void PART_Match_St_Edges_In_Gt(t_tree *gt, supert_tree *st)
689{
690  int i;
691
692  For(i,2*st->tree->n_otu-3) 
693    {
694      st->match_st_edge_in_gt[gt->dp][i] = NULL; 
695      st->match_gt_edge_in_st[gt->dp][i] = NULL;
696    }
697
698  For(i,st->tree->n_otu) 
699    if(st->match_st_node_in_gt[gt->dp][i])
700      {
701        PART_Match_St_Edges_In_Gt_Recurr(st->match_st_node_in_gt[gt->dp][i],
702                                       st->match_st_node_in_gt[gt->dp][i]->v[0],
703                                       st->tree->a_nodes[i],
704                                       st->tree->a_nodes[i]->v[0],
705                                       gt,st);
706        break;
707      }
708
709}
710
711//////////////////////////////////////////////////////////////
712//////////////////////////////////////////////////////////////
713
714
715void PART_Match_St_Edges_In_Gt_Recurr(t_node *a_gt, t_node *d_gt, t_node *a_st, t_node *d_st, t_tree *gt, supert_tree *st)
716{
717  t_edge *b_gt, *b_st;
718  int i,j;
719
720  b_gt = b_st = NULL;
721
722  if((st->match_st_node_in_gt[gt->dp][a_st->num] == a_gt) &&
723     (st->match_st_node_in_gt[gt->dp][d_st->num] == d_gt))
724    {
725      For(i,3) if((a_st->v[i]) && (a_st->v[i] == d_st)) {b_st = a_st->b[i]; break;}
726      For(i,3) if((a_gt->v[i]) && (a_gt->v[i] == d_gt)) {b_gt = a_gt->b[i]; break;}
727
728      st->match_st_edge_in_gt[gt->dp][b_st->num] = b_gt;
729      st->match_gt_edge_in_st[gt->dp][b_gt->num] = b_st;
730    }
731
732
733  if(!d_gt)
734    {
735      PhyML_Printf("\n");
736      PhyML_Printf("\n. a_gt->num = %d\n",a_gt->num);
737      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
738      Warn_And_Exit("");
739    }
740
741  if(d_gt->tax || d_st->tax) return;
742  else
743    {
744      if(st->match_st_node_in_gt[gt->dp][d_st->num] == d_gt)
745        {
746          For(i,3)
747            {
748              if(d_gt->v[i] != a_gt)
749                {
750                  For(j,3)
751                    {
752/*                    For(k,d_st->n_of_reachable_tips[j]) */
753/*                      if(!strcmp(d_gt->list_of_reachable_tips[i][0]->name,d_st->list_of_reachable_tips[j][k]->name)) */
754                          {
755                            PART_Match_St_Edges_In_Gt_Recurr(d_gt,d_gt->v[i],d_st,d_st->v[j],gt,st);
756                            break;
757                          }
758/*                    if(k != d_st->n_of_reachable_tips[j]) break; */
759                    }
760                }
761            }
762        }
763      else
764        {
765          For(i,3)
766            if(d_st->v[i] != a_st)
767              PART_Match_St_Edges_In_Gt_Recurr(a_gt,d_gt,d_st,d_st->v[i],gt,st);
768        }
769    }
770}
771
772//////////////////////////////////////////////////////////////
773//////////////////////////////////////////////////////////////
774
775
776void PART_Simu(supert_tree *st)
777{
778  int i,j,step,n_without_swap,it_lim_without_swap;
779  t_edge **sorted_b,*st_b,**tested_b;
780  int n_neg,n_tested,each;
781  phydbl lambda,old_loglk;
782
783  sorted_b = (t_edge **)mCalloc(st->tree->n_otu-3,sizeof(t_edge *));
784  tested_b = (t_edge **)mCalloc(st->tree->n_otu-3,sizeof(t_edge *));
785
786  For(i,st->n_part) Update_Dirs(st->treelist->tree[i]);
787  Update_Dirs(st->tree);
788
789  each                = 4;
790  step                = 0;
791  lambda              = .75;
792  n_tested            = 0;
793  n_without_swap      = 0;
794  old_loglk           = UNLIKELY; 
795  it_lim_without_swap = 2;
796  st_b                = NULL; 
797  do
798    {
799      For(i,st->n_part) Check_Dirs(st->treelist->tree[i]);
800
801      ++step;     
802      each--;
803
804/*       PART_Print_Bl(st); */
805
806      /* Compute the likelihood of the supertreee */
807      st->tree->c_lnL  = PART_Lk(st);
808      st->tree->c_pars = PART_Pars(st);
809/*       For(i,st->n_part) PhyML_Printf("\n. %s",Write_Tree(st->treelist->tree[i],NO)); */
810/*       PhyML_Printf("\n"); */
811
812      time(&(st->tree->t_current));
813      PhyML_Printf("\n. (%5d sec) [tot lnL=%15.5f] [# swaps=%3d]",
814             (int)(st->tree->t_current-st->tree->t_beg),
815             st->tree->c_lnL,n_tested);
816/*       For(i,st->n_part) PhyML_Printf("\n[gt %3d lnL=%15.5f]",i,st->treelist->tree[i]->c_lnL); */
817     
818      if((FABS(old_loglk-st->tree->c_lnL) < st->tree->mod->s_opt->min_diff_lk_global) || 
819         (n_without_swap > it_lim_without_swap)) break;
820
821      if(st->tree->c_lnL < old_loglk)
822        {
823          PhyML_Printf("\n. Moving backward (topology + branch lengths) \n");
824         
825          if(!PART_Mov_Backward_Topo_Bl(st,old_loglk,tested_b,n_tested))
826            Warn_And_Exit("\n. Err: mov_back failed\n");
827
828          if(!st->tree->n_swap) n_neg = 0;
829           
830          PART_Record_Br_Len(st);
831          For(i,st->n_part) Optimiz_All_Free_Param(st->treelist->tree[i],0);
832        }
833      else 
834        {
835          if(!each)
836            {
837              each = 4;
838              /* Markov model parameters are free to vary across data partitions */
839              For(i,st->n_part) Optimiz_All_Free_Param(st->treelist->tree[i],0);             
840              For(i,st->n_part) Set_Both_Sides(YES,st->treelist->tree[i]);
841              st->tree->c_lnL  = PART_Lk(st);
842              st->tree->c_pars = PART_Pars(st);
843            }
844         
845          old_loglk = st->tree->c_lnL;
846         
847
848          For(i,2*st->tree->n_otu-3) Init_NNI(st->tree->a_edges[i]->nni);
849
850          /* Test NNIs */
851          For(i,2*st->tree->n_otu-3)
852            {
853              st_b = st->tree->a_edges[i];
854              if((!st_b->left->tax) && (!st_b->rght->tax)) PART_NNI(st_b,st);
855            }
856         
857          /* Optimise external branch lengths */
858          For(i,2*st->tree->n_otu-3)
859            {
860              st_b = st->tree->a_edges[i];
861              if((st_b->left->tax) || (st_b->rght->tax))
862                {
863                  PART_Record_Br_Len(st);
864                  PART_Br_Len_Brent(st_b,0,st);
865                  For(j,st->n_part) st->bl0[st->bl_partition[j]][st_b->num] = st->bl[st->bl_partition[j]][st_b->num];
866                  st_b->nni->score     = .0;
867                  st_b->nni->best_conf =  0;
868                  PART_Restore_Br_Len(st);
869                  PART_Lk_At_Given_Edge(st_b,st);
870                }
871            }
872
873          /* Select and sort swaps */
874          n_neg = 0;
875          Select_Edges_To_Swap(st->tree,sorted_b,&n_neg);
876          Sort_Edges_NNI_Score(st->tree,sorted_b,n_neg);         
877
878          n_tested = 0;
879          For(i,(int)CEIL((phydbl)n_neg*(lambda)))
880            tested_b[n_tested++] = sorted_b[i];
881
882          if(n_tested > 0) n_without_swap = 0;
883          else             n_without_swap++;
884
885          PART_Record_Br_Len(st);
886         
887          /* Apply swaps */
888          PART_Make_N_Swap(tested_b,0,n_tested,st);
889
890          /* Update branch lengths (all edges first and then swaped edges) */
891          PART_Update_Bl(lambda,st);
892          PART_Update_Bl_Swaped(tested_b,n_tested,st);
893
894           
895        }
896    }
897  while(1);
898 
899  PhyML_Printf("\n\n. End of PART_Simu \n");
900  Free(sorted_b);
901  Free(tested_b);
902
903  if((n_without_swap > it_lim_without_swap))
904    {
905      PhyML_Printf("\n. Last optimization step...\n");
906      For(i,st->n_part) Round_Optimize(st->treelist->tree[i],st->treelist->tree[i]->data,ROUND_MAX);
907      st->tree->c_lnL = PART_Lk(st);
908      PART_Simu(st);
909    }
910}
911
912
913//////////////////////////////////////////////////////////////
914//////////////////////////////////////////////////////////////
915
916
917int PART_Mov_Backward_Topo_Bl(supert_tree *st, phydbl lk_old, t_edge **tested_b, int n_tested)
918{
919  int i,j,step,beg,end;
920  t_edge *st_b;
921  phydbl **l_init;
922  int dim;
923
924  l_init = (phydbl **)mCalloc(st->n_part,sizeof(phydbl *));
925  For(i,st->n_part) l_init[i] = (phydbl *)mCalloc(2*st->tree->n_otu-3,sizeof(phydbl ));
926
927  For(i,2*st->tree->n_otu-3) For(j,st->n_part) l_init[st->bl_partition[j]][i] = st->bl[st->bl_partition[j]][i];
928
929  step = 2;
930  do
931    {
932      For(i,2*st->tree->n_otu-3)
933        {
934          For(j,st->n_part)
935            {
936              st->bl[st->bl_partition[j]][i] = st->bl_cpy[st->bl_partition[j]][i] + 
937                (1./step) * (l_init[st->bl_partition[j]][i] - st->bl_cpy[st->bl_partition[j]][i]);
938/*            st->bl[st->bl_partition[j]][i] = st->bl_cpy[st->bl_partition[j]][i]; */
939            }
940        }
941     
942      beg = (int)FLOOR((phydbl)n_tested/(step-1));
943      end = 0;
944      st_b = NULL;
945      dim = 2*st->tree->n_otu-2;
946
947      for(i=beg-1;i>=end;i--)
948        {
949          st_b = tested_b[i];
950
951          PART_Swap(st_b->nni->swap_node_v2->v[st->tree->t_dir[st_b->nni->swap_node_v2->num*dim+st_b->nni->swap_node_v1->num]],
952                  st_b->nni->swap_node_v2,
953                  st_b->nni->swap_node_v3,
954                  st_b->nni->swap_node_v3->v[st->tree->t_dir[st_b->nni->swap_node_v3->num*dim+st_b->nni->swap_node_v4->num]],
955                  st);
956
957          Swap(st_b->nni->swap_node_v2->v[st->tree->t_dir[st_b->nni->swap_node_v2->num*dim+st_b->nni->swap_node_v1->num]],
958               st_b->nni->swap_node_v2,
959               st_b->nni->swap_node_v3,
960               st_b->nni->swap_node_v3->v[st->tree->t_dir[st_b->nni->swap_node_v3->num*dim+st_b->nni->swap_node_v4->num]],
961               st->tree);
962
963          PART_Do_Mapping(st);
964
965        }
966
967      beg = 0;
968      end = (int)FLOOR((phydbl)n_tested/step);
969      st_b = NULL;
970      dim = 2*st->tree->n_otu-2;
971
972      for(i=beg;i<end;i++)
973        {
974          st_b = tested_b[i];
975         
976          PART_Swap(st_b->nni->swap_node_v2->v[st->tree->t_dir[st_b->nni->swap_node_v2->num*dim+st_b->nni->swap_node_v1->num]],
977                  st_b->nni->swap_node_v2,
978                  st_b->nni->swap_node_v3,
979                  st_b->nni->swap_node_v3->v[st->tree->t_dir[st_b->nni->swap_node_v3->num*dim+st_b->nni->swap_node_v4->num]],
980                  st);
981
982          Swap(st_b->nni->swap_node_v2->v[st->tree->t_dir[st_b->nni->swap_node_v2->num*dim+st_b->nni->swap_node_v1->num]],
983               st_b->nni->swap_node_v2,
984               st_b->nni->swap_node_v3,
985               st_b->nni->swap_node_v3->v[st->tree->t_dir[st_b->nni->swap_node_v3->num*dim+st_b->nni->swap_node_v4->num]],
986               st->tree);
987
988          PART_Do_Mapping(st);
989
990        }
991     
992      if(!end) st->tree->n_swap = 0;     
993
994      PART_Lk(st);
995
996      PhyML_Printf("\n. lnL = %15.5f",st->tree->c_lnL);
997      step++;
998    }
999  while((st->tree->c_lnL < lk_old) && (step < 100));
1000
1001  if(step == 100)
1002    {
1003      For(i,2*st->tree->n_otu-3) For(j,st->n_part)
1004        st->bl[st->bl_partition[j]][i] = st->bl_cpy[st->bl_partition[j]][i];
1005    }
1006
1007  st->tree->n_swap = 0;
1008  For(i,2*st->tree->n_otu-3) 
1009    {
1010      if(st->tree->a_edges[i]->nni->score < 0.0) st->tree->n_swap++;
1011      st->tree->a_edges[i]->nni->score = +1.0;
1012    }
1013
1014  PART_Lk(st);
1015
1016  if(st->tree->c_lnL > lk_old)                                                    return  1;
1017  else if(FABS(st->tree->c_lnL-lk_old) < st->tree->mod->s_opt->min_diff_lk_local) return -1;
1018  else                                                                            return  0;
1019
1020  For(i,st->n_part) Free(l_init[i]);
1021  Free(l_init);
1022}
1023
1024//////////////////////////////////////////////////////////////
1025//////////////////////////////////////////////////////////////
1026
1027
1028void PART_Check_Extra_Taxa(supert_tree *st)
1029{
1030  int i,j,k;
1031  int sum;
1032  int *st_taxa;
1033
1034  st_taxa = (int *)mCalloc(st->tree->n_otu,sizeof(int));
1035
1036  For(i,st->tree->n_otu)
1037    {
1038      For(j,st->n_part)
1039        {
1040          For(k,st->treelist->tree[j]->n_otu) 
1041            if(!strcmp(st->treelist->tree[j]->a_nodes[k]->name,st->tree->a_nodes[i]->name)) break;
1042          if(k != st->treelist->tree[j]->n_otu) { st_taxa[i] = 1; break; }
1043        }
1044    }
1045
1046  sum = 0;
1047  For(i,st->tree->n_otu) if(st_taxa[i]) sum++;
1048  if(sum != st->tree->n_otu) 
1049    {
1050      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1051      Warn_And_Exit("");
1052    }
1053  Free(st_taxa);
1054}
1055
1056//////////////////////////////////////////////////////////////
1057//////////////////////////////////////////////////////////////
1058
1059
1060int PART_Get_Species_Found_In_St(supert_tree *st, calign *data)
1061{
1062  int i,j;
1063 
1064  For(i,data->n_otu)
1065    {
1066      For(j,st->tree->n_otu)
1067        {
1068          if(!strcmp(data->c_seq[i]->name,st->tree->a_nodes[j]->name))
1069            {
1070              break;
1071            }
1072        }
1073      if(j == st->tree->n_otu) return 0;
1074    }
1075  return 1;
1076
1077}
1078
1079//////////////////////////////////////////////////////////////
1080//////////////////////////////////////////////////////////////
1081
1082
1083void PART_Map_St_Nodes_In_Gt(t_tree *gt, supert_tree *st)
1084{
1085  int i;
1086
1087  For(i,2*st->tree->n_otu-2)
1088    {
1089      st->map_st_node_in_gt[gt->dp][i][0][0] = NULL;
1090      st->map_st_node_in_gt[gt->dp][i][1][0] = NULL;
1091      st->map_st_node_in_gt[gt->dp][i][2][0] = NULL;
1092
1093      st->map_st_node_in_gt[gt->dp][i][0][1] = NULL;
1094      st->map_st_node_in_gt[gt->dp][i][1][1] = NULL;
1095      st->map_st_node_in_gt[gt->dp][i][2][1] = NULL;
1096    }
1097
1098 
1099  /* Root */
1100  PART_Map_St_Nodes_In_Gt_One_Edge(st->tree->a_nodes[0]->v[0],
1101                                 st->tree->a_nodes[0],
1102                                 st->tree->a_nodes[0]->b[0],
1103                                 gt,st);
1104
1105  /* Internal nodes */
1106  PART_Map_St_Nodes_In_Gt_Post(st->tree->a_nodes[0],st->tree->a_nodes[0]->v[0],gt,st);
1107  PART_Map_St_Nodes_In_Gt_Pre (st->tree->a_nodes[0],st->tree->a_nodes[0]->v[0],gt,st);
1108 
1109  /* Root */
1110  PART_Map_St_Nodes_In_Gt_One_Edge(st->tree->a_nodes[0],
1111                                 st->tree->a_nodes[0]->v[0],
1112                                 st->tree->a_nodes[0]->b[0],
1113                                 gt,st);
1114 
1115}
1116
1117//////////////////////////////////////////////////////////////
1118//////////////////////////////////////////////////////////////
1119
1120
1121void PART_Map_St_Nodes_In_Gt_Post(t_node *a_st, t_node *d_st, t_tree *gt, supert_tree *st)
1122{
1123  int i;
1124
1125  if(d_st->tax) return;
1126  else
1127    {
1128      For(i,3)
1129        if(d_st->v[i] != a_st)
1130          PART_Map_St_Nodes_In_Gt_Post(d_st,d_st->v[i],gt,st);
1131     
1132      For(i,3)
1133        if(d_st->v[i] != a_st)
1134          {
1135            PART_Map_St_Nodes_In_Gt_One_Edge(d_st,d_st->v[i],d_st->b[i],gt,st);
1136          }
1137    }
1138}
1139
1140//////////////////////////////////////////////////////////////
1141//////////////////////////////////////////////////////////////
1142
1143
1144void PART_Map_St_Nodes_In_Gt_Pre(t_node *a_st, t_node *d_st, t_tree *gt, supert_tree *st)
1145{
1146  int i;
1147
1148  if(d_st->tax) return;
1149  else
1150    {
1151      For(i,3)
1152        {
1153          if(d_st->v[i] != a_st)
1154            {         
1155              PART_Map_St_Nodes_In_Gt_One_Edge(d_st->v[i],d_st,d_st->b[i],gt,st);
1156              PART_Map_St_Nodes_In_Gt_Pre(d_st,d_st->v[i],gt,st);
1157            }
1158        }
1159    }
1160}
1161
1162//////////////////////////////////////////////////////////////
1163//////////////////////////////////////////////////////////////
1164
1165
1166void PART_Map_St_Nodes_In_Gt_One_Edge(t_node *a_st, t_node *d_st, t_edge *b_st, t_tree *gt, supert_tree *st)
1167{
1168  if(d_st->tax)
1169    {
1170#ifdef DEBUG
1171      if(b_st->rght != d_st)
1172        {
1173          PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1174          Warn_And_Exit("");
1175        }
1176#endif
1177     
1178      st->map_st_node_in_gt[gt->dp][d_st->num][0][0] = st->match_st_node_in_gt[gt->dp][d_st->num];
1179    }
1180  else
1181    {
1182      t_node **list_of_nodes_d, **list_of_nodes_v1, **list_of_nodes_v2;
1183      int dir1, dir2;
1184      int i;
1185
1186      list_of_nodes_d  = NULL;
1187      list_of_nodes_v1 = NULL;
1188      list_of_nodes_v2 = NULL;
1189
1190      dir1 = dir2 = -1;
1191      For(i,3) 
1192        {
1193          if(d_st->v[i] != a_st) (dir1 < 0)?(dir1 = i):(dir2 = i);
1194          else list_of_nodes_d = st->map_st_node_in_gt[gt->dp][d_st->num][i];
1195        }
1196
1197      For(i,3) 
1198        if((d_st->v[dir1]->v[i]) && (d_st->v[dir1]->v[i] == d_st)) 
1199          {
1200            list_of_nodes_v1 = st->map_st_node_in_gt[gt->dp][d_st->v[dir1]->num][i];
1201            break;
1202          }
1203
1204      For(i,3) 
1205        if((d_st->v[dir2]->v[i]) && (d_st->v[dir2]->v[i] == d_st)) 
1206          {
1207            list_of_nodes_v2 = st->map_st_node_in_gt[gt->dp][d_st->v[dir2]->num][i];
1208            break;
1209          }
1210
1211      /* d_st matches one t_node in gt */
1212      if(st->match_st_node_in_gt[gt->dp][d_st->num])
1213        {
1214          list_of_nodes_d[0] = st->match_st_node_in_gt[gt->dp][d_st->num];
1215          list_of_nodes_d[1] = NULL;
1216        }
1217      else
1218        {
1219          /* list_of_nodes = union of  list_of_nodes_v1  &  list_of_nodes_v2 */
1220         
1221          if(!list_of_nodes_v1[0])
1222            {
1223              list_of_nodes_d[0] = list_of_nodes_v2[0];
1224              list_of_nodes_d[1] = list_of_nodes_v2[1];
1225            }
1226          else if(!list_of_nodes_v2[0])
1227            {
1228              list_of_nodes_d[0] = list_of_nodes_v1[0];
1229              list_of_nodes_d[1] = list_of_nodes_v1[1];
1230            }
1231          else
1232            {
1233              list_of_nodes_d[0] = list_of_nodes_v1[0];
1234              list_of_nodes_d[1] = list_of_nodes_v2[0];
1235
1236              if(list_of_nodes_v1[1] || list_of_nodes_v2[1])
1237                {
1238                  Print_Node(st->tree->a_nodes[0],
1239                             st->tree->a_nodes[0]->v[0],
1240                             st->tree);
1241                 
1242                  PhyML_Printf("\n\n--------------------------\n\n");
1243                  Print_Node(gt->a_nodes[0],
1244                             gt->a_nodes[0]->v[0],
1245                             gt);
1246
1247                  PhyML_Printf("\n\n--------------------------\n\n");
1248                 
1249                  PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1250                  Warn_And_Exit("");
1251                }
1252            }
1253        }
1254    }
1255}
1256
1257//////////////////////////////////////////////////////////////
1258//////////////////////////////////////////////////////////////
1259
1260
1261void PART_Map_St_Edges_In_Gt(t_tree *gt, supert_tree *st)
1262{
1263  int i,j;
1264  t_edge *st_b;
1265  t_node *gt_a, *gt_d;
1266 
1267  gt_a = NULL;
1268  gt_d = NULL;
1269
1270  For(i,2*st->tree->n_otu-3) st->map_st_edge_in_gt[gt->dp][i] = NULL;
1271
1272  For(i,2*st->tree->n_otu-3)
1273    {
1274      st_b = st->tree->a_edges[i];
1275
1276      if(!st->map_st_node_in_gt[gt->dp][st_b->left->num][st_b->l_r][0])
1277        {
1278          gt_a = st->map_st_node_in_gt[gt->dp][st_b->rght->num][st_b->r_l][0];
1279          gt_d = st->map_st_node_in_gt[gt->dp][st_b->rght->num][st_b->r_l][1];
1280
1281          For(j,3)
1282            {
1283              if((gt_a->v[j]) && (gt_a->v[j] == gt_d))
1284                {
1285                  st->map_st_edge_in_gt[gt->dp][st_b->num] = gt_a->b[j];
1286                  break;
1287                }
1288            }
1289        }
1290      else if(!st->map_st_node_in_gt[gt->dp][st_b->rght->num][st_b->r_l][0])
1291        {
1292          gt_a = st->map_st_node_in_gt[gt->dp][st_b->left->num][st_b->l_r][0];
1293          gt_d = st->map_st_node_in_gt[gt->dp][st_b->left->num][st_b->l_r][1];
1294
1295          For(j,3)
1296            {
1297              if((gt_a->v[j]) && (gt_a->v[j] == gt_d))
1298                {
1299                  st->map_st_edge_in_gt[gt->dp][st_b->num] = gt_a->b[j];
1300                  break;
1301                }
1302            }
1303        }
1304      else
1305        {         
1306          gt_a = st->map_st_node_in_gt[gt->dp][st_b->left->num][st_b->l_r][0];
1307          gt_d = st->map_st_node_in_gt[gt->dp][st_b->rght->num][st_b->r_l][0];
1308
1309          #ifdef DEBUG
1310          if((!gt_a) || (!gt_d))
1311            {
1312              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1313              Warn_And_Exit("");
1314            }
1315          #endif
1316
1317          For(j,3)
1318            {
1319              if((gt_a->v[j]) && (gt_a->v[j] == gt_d))
1320                {
1321                  st->map_st_edge_in_gt[gt->dp][st_b->num] = gt_a->b[j];
1322                  break;
1323                }
1324            }
1325        }
1326    }
1327}
1328
1329//////////////////////////////////////////////////////////////
1330//////////////////////////////////////////////////////////////
1331
1332
1333void PART_Map_Gt_Edges_In_St(t_tree *gt, supert_tree *st)
1334{
1335  int i;
1336  t_edge *st_b, *gt_b;
1337 
1338  For(i,2*st->tree->n_otu-3) st->size_map_gt_edge_in_st[gt->dp][i] = 0;
1339
1340  st_b = NULL;
1341  gt_b = NULL;
1342  For(i,2*st->tree->n_otu-3)
1343    {
1344      st_b = st->tree->a_edges[i];
1345      gt_b = st->map_st_edge_in_gt[gt->dp][st_b->num];
1346
1347      if(gt_b)
1348        {
1349          st->map_gt_edge_in_st[gt->dp][gt_b->num][st->size_map_gt_edge_in_st[gt->dp][gt_b->num]] = st_b;
1350          st->size_map_gt_edge_in_st[gt->dp][gt_b->num]++;
1351        }
1352    }
1353}
1354
1355//////////////////////////////////////////////////////////////
1356//////////////////////////////////////////////////////////////
1357
1358
1359int PART_Pars(supert_tree *st)
1360{
1361  int i;
1362 
1363  st->tree->c_pars = 0;
1364  For(i,st->n_part) 
1365    {
1366      Set_Both_Sides(YES,st->treelist->tree[i]);         
1367      Pars(NULL,st->treelist->tree[i]);
1368      st->tree->c_pars += st->treelist->tree[i]->c_pars;
1369    }
1370
1371  return st->tree->c_pars;
1372}
1373
1374//////////////////////////////////////////////////////////////
1375//////////////////////////////////////////////////////////////
1376
1377
1378int PART_Spr(phydbl init_lnL, supert_tree *st)
1379{
1380  int gt;
1381  int i;
1382  t_edge *pruned;
1383  int best_move;
1384  t_node *gt_a, *gt_d;
1385
1386  st->tree->n_root     = st->tree->a_nodes[0];
1387  pruned               = NULL;
1388  gt_a                 = NULL;
1389  gt_d                 = NULL;
1390 
1391  Set_Both_Sides(YES,st->tree);
1392
1393  For(i,2*st->tree->n_otu-3)
1394    {
1395      pruned            = st->tree->a_edges[i];
1396      st->tree->n_moves = 0;
1397
1398      Reset_Spr_List(st->tree);
1399      For(gt,st->n_part) Reset_Spr_List(st->treelist->tree[gt]);
1400     
1401      if(!pruned->rght->tax)
1402        {
1403          For(gt,st->n_part)
1404            {
1405              /* Check constraints at prune site on gt tree */
1406              gt_a = st->map_st_node_in_gt[gt][pruned->rght->num][pruned->r_l][0];
1407              gt_d = st->map_st_node_in_gt[gt][pruned->left->num][pruned->l_r][0];
1408
1409              if((gt_a) && (gt_d) && (!st->map_st_edge_in_gt[gt][pruned->num]->rght->tax))
1410                {
1411                  Test_All_Spr_Targets(st->map_st_edge_in_gt[gt][pruned->num],
1412                                       st->map_st_edge_in_gt[gt][pruned->num]->rght,
1413                                       st->treelist->tree[gt]);
1414                }
1415            }
1416        }
1417
1418      if(!pruned->left->tax)
1419        {
1420          For(gt,st->n_part)
1421            {
1422              /* Check constraints at prune site on gt tree */       
1423              gt_a = st->map_st_node_in_gt[gt][pruned->rght->num][pruned->r_l][0];
1424              gt_d = st->map_st_node_in_gt[gt][pruned->left->num][pruned->l_r][0];
1425
1426              if((gt_a) && (gt_d) && (!st->map_st_edge_in_gt[gt][pruned->num]->left->tax))
1427                {
1428                  Test_All_Spr_Targets(st->map_st_edge_in_gt[gt][pruned->num],
1429                                       st->map_st_edge_in_gt[gt][pruned->num]->left,
1430                                       st->treelist->tree[gt]);
1431                }
1432            }
1433        }
1434
1435     
1436      if(!pruned->left->tax)
1437        {
1438          PART_Test_All_Spr_Targets(st->tree->a_edges[i],
1439                                  st->tree->a_edges[i]->left,
1440                                  st);     
1441        }
1442     
1443      if(!pruned->rght->tax)
1444        {
1445          PART_Test_All_Spr_Targets(st->tree->a_edges[i],
1446                                  st->tree->a_edges[i]->rght,
1447                                  st);     
1448        }
1449
1450
1451      if(st->tree->n_moves)
1452        {
1453          best_move = PART_Test_List_Of_Regraft_Pos(st->tree->spr_list,
1454                                                  (int)CEIL(0.1*(st->tree->n_moves)),
1455                                                  st);   
1456         
1457          if(st->tree->spr_list[best_move]->lnL > init_lnL)
1458            {
1459              PART_Try_One_Spr_Move(st->tree->spr_list[best_move],st);
1460            }
1461          else
1462            {
1463              Set_Both_Sides(YES,st->tree);
1464              st->tree->c_lnL  = PART_Lk(st);
1465              st->tree->c_pars = PART_Pars(st);       
1466            }
1467        }
1468    }
1469  return 1;
1470}
1471
1472//////////////////////////////////////////////////////////////
1473//////////////////////////////////////////////////////////////
1474
1475
1476void PART_Speed_Spr(supert_tree *st)
1477{
1478  int step;
1479  int gt;
1480  phydbl old_lnL;
1481
1482  Make_Spr_List(st->tree);
1483  For(gt,st->n_part) Make_Spr_List(st->treelist->tree[gt]);
1484
1485  Set_Both_Sides(YES,st->tree); 
1486  For(gt,st->n_part) 
1487    {
1488      Set_Both_Sides(YES,st->treelist->tree[gt]);
1489      Record_Br_Len(st->treelist->tree[gt]);
1490    }
1491 
1492  st->tree->c_pars = PART_Pars(st);
1493  st->tree->c_lnL  = PART_Lk(st);
1494 
1495
1496  st->tree->best_lnL = st->tree->c_lnL;
1497  old_lnL            = st->tree->c_lnL;
1498  step               = 0;
1499  do
1500    {
1501      ++step;
1502
1503      PhyML_Printf("\n. Starting a SPR cycle... \n");
1504
1505      old_lnL = st->tree->c_lnL;
1506
1507      st->tree->n_improvements         = 0;
1508      st->tree->perform_spr_right_away = 1;
1509      PART_Spr(UNLIKELY,st);
1510
1511 
1512      time(&(st->tree->t_current));     
1513      PhyML_Printf("\n. (%5d sec) [00] [%10.2f] [%5d]\n",
1514             (int)(st->tree->t_current-st->tree->t_beg),
1515             PART_Lk(st),PART_Pars(st));
1516
1517      /* Optimise parameters of the Markov t_mod */
1518      For(gt,st->n_part) Optimiz_All_Free_Param(st->treelist->tree[gt],
1519                                              st->treelist->tree[gt]->mod->s_opt->print);
1520
1521      time(&(st->tree->t_current));     
1522      PhyML_Printf("\n. (%5d sec) [ 0] [%10.2f] [%5d]\n",
1523             (int)(st->tree->t_current-st->tree->t_beg),
1524             PART_Lk(st),PART_Pars(st));
1525
1526      /* Optimise branch lengths */
1527      For(gt,st->n_part)
1528        {
1529          Optimize_Br_Len_Serie(st->treelist->tree[gt]);
1530        }
1531
1532
1533      /* Update partial likelihoods & parsimony */
1534      Set_Both_Sides(YES,st->tree); 
1535      st->tree->c_pars = PART_Pars(st);
1536      st->tree->c_lnL  = PART_Lk(st);
1537     
1538     
1539      time(&(st->tree->t_current));     
1540      PhyML_Printf("\n. (%5d sec) [**] [%10.2f] [%5d]\n",
1541             (int)(st->tree->t_current-st->tree->t_beg),
1542             st->tree->c_lnL,st->tree->c_pars);
1543
1544      /* Record the current best log-likleihood  */
1545      st->tree->best_lnL = st->tree->c_lnL;
1546
1547      if(st->tree->c_lnL < old_lnL)
1548        {
1549          PhyML_Printf("\n. old_lnL = %f c_lnL = %f\n",old_lnL,st->tree->c_lnL); 
1550          PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1551          Warn_And_Exit("");
1552        }
1553
1554      /* Record the current best branch lengths  */
1555      For(gt,st->n_part) Record_Br_Len(st->treelist->tree[gt]);
1556
1557      /* Exit if no improvements after complete optimization */
1558      if((!st->tree->n_improvements) || 
1559         (FABS(old_lnL-st->tree->c_lnL) < st->tree->mod->s_opt->min_diff_lk_global)) break;
1560           
1561    }while(1);
1562 
1563}
1564
1565//////////////////////////////////////////////////////////////
1566//////////////////////////////////////////////////////////////
1567
1568
1569
1570void PART_Test_All_Spr_Targets(t_edge *pruned, t_node *n_link, supert_tree *st)
1571{
1572  int i,j;
1573
1574  For(i,3)
1575    {
1576      if(n_link->b[i] != pruned)
1577        {
1578          For(j,3)
1579            {
1580              if((n_link->v[i]->v[j]) && (n_link->v[i]->v[j] != n_link))
1581                {
1582                  PART_Test_One_Spr_Target_Recur(n_link->v[i],n_link->v[i]->v[j],n_link->v[i]->b[j],pruned,n_link,st);
1583                }
1584            }
1585        }
1586    }
1587}
1588
1589//////////////////////////////////////////////////////////////
1590//////////////////////////////////////////////////////////////
1591
1592
1593void PART_Test_One_Spr_Target_Recur(t_node *a, t_node *d, t_edge *target, t_edge *pruned, t_node *n_link, supert_tree *st)
1594{
1595
1596  PART_Test_One_Spr_Target(pruned,target,n_link,st);
1597
1598  if(d->tax) return;
1599  else
1600    {
1601      int i;
1602
1603      For(i,3)
1604        if(d->v[i] != a)
1605          {
1606            PART_Test_One_Spr_Target_Recur(d,d->v[i],d->b[i],pruned,n_link,st);
1607          }
1608    }
1609}
1610
1611//////////////////////////////////////////////////////////////
1612//////////////////////////////////////////////////////////////
1613
1614
1615void PART_Test_One_Spr_Target(t_edge *st_p, t_edge *st_t, t_node *n_link, supert_tree *st) 
1616{
1617  int gt, move;
1618
1619  st->tree->n_moves++;
1620  st->tree->spr_list[st->tree->size_spr_list]->b_target      = st_t;
1621  st->tree->spr_list[st->tree->size_spr_list]->n_link        = n_link;
1622  st->tree->spr_list[st->tree->size_spr_list]->n_opp_to_link = (n_link == st_p->left)?(st_p->rght):(st_p->left);
1623  st->tree->spr_list[st->tree->size_spr_list]->b_opp_to_link = st_p;
1624  st->tree->spr_list[st->tree->size_spr_list]->pars          = 0;
1625
1626  For(gt,st->n_part)
1627    {
1628      move = Map_Spr_Move(st_p,st_t,n_link,st->treelist->tree[gt],st);
1629     
1630      if(move > -1)
1631        st->tree->spr_list[st->tree->size_spr_list]->pars += st->treelist->tree[gt]->spr_list[move]->pars;
1632      else if(move == -1 || move == -2)
1633        st->tree->spr_list[st->tree->size_spr_list]->pars += st->treelist->tree[gt]->c_pars;
1634    }
1635
1636  Include_One_Spr_To_List_Of_Spr(st->tree->spr_list[st->tree->size_spr_list],st->tree);
1637}
1638
1639//////////////////////////////////////////////////////////////
1640//////////////////////////////////////////////////////////////
1641
1642
1643int Map_Spr_Move(t_edge *st_pruned, t_edge *st_target, t_node *st_link, t_tree *gt, supert_tree *st)
1644{
1645  int i;
1646  t_edge *gt_pruned, *gt_target;
1647  t_node *gt_link, *gt_a, *gt_d;
1648
1649  gt_pruned = NULL;
1650  gt_target = NULL;
1651  gt_link   = NULL;
1652  gt_a      = NULL;
1653  gt_d      = NULL;
1654
1655  /* Check contraints at prune and regraft sites on the gt tree */
1656
1657  /* st_pruned is not on a path that matches a branch in gt */
1658  gt_a = st->map_st_node_in_gt[gt->dp][st_pruned->left->num][st_pruned->l_r][0];
1659  gt_d = st->map_st_node_in_gt[gt->dp][st_pruned->rght->num][st_pruned->r_l][0];
1660
1661  if((!gt_a) || (!gt_d)) return -1;
1662  else
1663    {
1664      /* which gt nodes matches st_link ? */
1665      gt_link = (st_pruned->left == st_link)?(gt_a):(gt_d);
1666     
1667      if(gt_link->tax) return -1;
1668      else
1669        {
1670          gt_pruned = st->map_st_edge_in_gt[gt->dp][st_pruned->num];
1671          gt_target = st->map_st_edge_in_gt[gt->dp][st_target->num];
1672           
1673          if((gt_pruned->left == gt_target->left) ||
1674             (gt_pruned->left == gt_target->rght) ||
1675             (gt_pruned->rght == gt_target->left) ||
1676             (gt_pruned->rght == gt_target->rght)) return -1;
1677          else
1678            {
1679              For(i,gt->size_spr_list)
1680                {
1681                  if((gt_pruned == gt->spr_list[i]->b_opp_to_link) && (gt_target == gt->spr_list[i]->b_target))
1682                    return i;
1683                }
1684            }
1685        }
1686    }
1687  return -2;
1688}
1689
1690//////////////////////////////////////////////////////////////
1691//////////////////////////////////////////////////////////////
1692
1693
1694int PART_Test_List_Of_Regraft_Pos(t_spr **st_spr_list, int list_size, supert_tree *st)
1695{
1696
1697  int i,j,best_move;
1698  t_spr *move;
1699  t_edge *init_target, *b_residual;
1700  phydbl best_lnL, init_lnL;
1701  int dir_v0, dir_v1, dir_v2;
1702  int gt;
1703  int move_num;
1704 
1705
1706  best_lnL = UNLIKELY;
1707  init_target = b_residual = NULL;
1708  best_move = -1;
1709
1710#ifdef DEBUG
1711  if(!list_size)
1712    {
1713      PhyML_Printf("\n\n. List size is 0 !");
1714      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1715      Warn_And_Exit(""); 
1716    }
1717#endif
1718 
1719  init_lnL = UNLIKELY;
1720
1721  For(i,list_size)
1722    {
1723      st->tree->spr_list[i]->lnL = .0;
1724
1725      For(gt,st->n_part) 
1726        {
1727          move_num = Map_Spr_Move(st->tree->spr_list[i]->b_opp_to_link,
1728                                  st->tree->spr_list[i]->b_target,
1729                                  st->tree->spr_list[i]->n_link,
1730                                  st->treelist->tree[gt],st);
1731
1732          if(move_num > -1)
1733            {
1734              move = st->treelist->tree[gt]->spr_list[move_num];
1735
1736              if(move->b_target)
1737                {
1738                  init_lnL = st->treelist->tree[gt]->c_lnL;
1739
1740                  /* Record t_edge lengths */
1741                  Record_Br_Len(st->treelist->tree[gt]);
1742                 
1743                  /* Prune subtree */
1744                  Prune_Subtree(move->n_link,
1745                                move->n_opp_to_link,                       
1746                                &init_target,
1747                                &b_residual,
1748                                st->treelist->tree[gt]);
1749
1750                  /* Rough optimisation of the branch length */
1751                  Fast_Br_Len(init_target,st->treelist->tree[gt],0);
1752                 
1753                  /* Update the change proba matrix at prune position */
1754                  Update_PMat_At_Given_Edge(init_target,st->treelist->tree[gt]);
1755             
1756                  /* Update partial likelihood along the path from the prune to
1757                     the regraft position */
1758                  Update_P_Lk_Along_A_Path(move->path,move->depth_path,st->treelist->tree[gt]);
1759
1760                  /* Regraft subtree */
1761                  Graft_Subtree(move->b_target,move->n_link,b_residual,st->treelist->tree[gt]);
1762             
1763                  /* Estimate the three t_edge lengths at the regraft site */
1764                  Triple_Dist(move->n_link,st->treelist->tree[gt],-1);
1765             
1766                  /* Update the transition proba matrices along edges defining
1767                     the regraft site */
1768                  For(j,3)
1769                    if(move->n_link->v[j] != move->n_opp_to_link)
1770                      Update_PMat_At_Given_Edge(move->n_link->b[j],st->treelist->tree[gt]);
1771             
1772                  /* Compute the likelihood */
1773                  Update_P_Lk(st->treelist->tree[gt],
1774                              move->b_opp_to_link,
1775                              move->n_link);
1776                 
1777                  move->lnL = Lk(move->b_opp_to_link,st->treelist->tree[gt]);
1778
1779                 
1780                  st->tree->spr_list[i]->lnL += move->lnL;
1781
1782                  /* Record branch lengths */
1783                  dir_v1 = dir_v2 = dir_v0 = -1;
1784                  For(j,3)
1785                    {
1786                      if(move->n_link->v[j] == move->n_opp_to_link) dir_v0 = j;
1787                      else if(dir_v1 < 0)                           dir_v1 = j;
1788                      else                                          dir_v2 = j;
1789                    }
1790                 
1791                  move->l0 = move->n_link->b[dir_v0]->l->v;
1792                 
1793                  if(move->n_link->v[dir_v1]->num > move->n_link->v[dir_v2]->num)
1794                    {
1795                      move->l1 = move->n_link->b[dir_v2]->l->v;
1796                      move->l2 = move->n_link->b[dir_v1]->l->v;
1797                    }
1798                  else
1799                    {
1800                      move->l1 = move->n_link->b[dir_v1]->l->v;
1801                      move->l2 = move->n_link->b[dir_v2]->l->v;
1802                    }
1803                         
1804                  /* Regraft the subtree at its original position */
1805                  Prune_Subtree(move->n_link,
1806                                move->n_opp_to_link,
1807                                &move->b_target,
1808                                &b_residual,
1809                                st->treelist->tree[gt]);
1810                 
1811                  Graft_Subtree(init_target,
1812                                move->n_link,
1813                                b_residual,
1814                                st->treelist->tree[gt]);
1815                 
1816                  /* Restore branch lengths */
1817                  Restore_Br_Len(st->treelist->tree[gt]);
1818             
1819                  /* Update relevant change proba matrices */
1820                  Update_PMat_At_Given_Edge(move->b_target,st->treelist->tree[gt]);
1821                  For(j,3) Update_PMat_At_Given_Edge(move->n_link->b[j],st->treelist->tree[gt]);
1822                 
1823                  /* Update relevant partial likelihoods */
1824                  For(j,3) Update_P_Lk(st->treelist->tree[gt],move->n_link->b[j],move->n_link);
1825                 
1826                  st->treelist->tree[gt]->c_lnL = init_lnL;
1827                }
1828            }
1829          else
1830            {
1831              st->tree->spr_list[i]->lnL += st->treelist->tree[gt]->c_lnL;
1832            }
1833        }
1834
1835      if(st->tree->spr_list[i]->lnL > best_lnL)
1836        {
1837          best_lnL  = st->tree->spr_list[i]->lnL;
1838          best_move = i;
1839        }
1840    }
1841
1842  return best_move; 
1843
1844}
1845
1846
1847//////////////////////////////////////////////////////////////
1848//////////////////////////////////////////////////////////////
1849
1850
1851int PART_Try_One_Spr_Move(t_spr *st_move, supert_tree *st)
1852{
1853  int j;
1854  t_spr **gt_move;
1855  t_edge **init_target, **b_residual;
1856  int dir_v0, dir_v1, dir_v2;
1857  int gt;
1858  int gt_move_num;
1859  int n_moves;
1860
1861 
1862  init_target = (t_edge **)mCalloc(st->n_part,sizeof(t_edge *));
1863  b_residual  = (t_edge **)mCalloc(st->n_part,sizeof(t_edge *));
1864  gt_move     = (t_spr **)mCalloc(st->n_part,sizeof(t_spr *));
1865
1866
1867  n_moves = 0;
1868  For(gt,st->n_part) 
1869    {
1870      gt_move_num = Map_Spr_Move(st_move->b_opp_to_link,
1871                                 st_move->b_target,
1872                                 st_move->n_link,
1873                                 st->treelist->tree[gt],st);
1874     
1875      if(gt_move_num > -1)
1876        {
1877          n_moves++;
1878
1879          gt_move[gt] = st->treelist->tree[gt]->spr_list[gt_move_num];
1880         
1881          if(gt_move[gt]->b_target)
1882            {
1883              /* Record t_edge lengths */
1884              Record_Br_Len(st->treelist->tree[gt]);
1885
1886              /* Prune subtree */
1887              Prune_Subtree(gt_move[gt]->n_link,
1888                            gt_move[gt]->n_opp_to_link,
1889                            &(init_target[gt]),
1890                            &(b_residual[gt]),
1891                            st->treelist->tree[gt]);
1892             
1893              /* Rough optimisation of the branch length */
1894              Fast_Br_Len(init_target[gt],st->treelist->tree[gt],0);
1895             
1896              /* Update the change proba matrix at prune position */
1897              Update_PMat_At_Given_Edge(init_target[gt],st->treelist->tree[gt]); /* TO DO : NECESSARY ?? */
1898             
1899              /* Update partial likelihood along the path from the prune to
1900                 the regraft position */
1901              Update_P_Lk_Along_A_Path(gt_move[gt]->path,gt_move[gt]->depth_path,st->treelist->tree[gt]); /* TO DO : NECESSARY ?? */
1902             
1903              /* Regraft subtree */
1904              Graft_Subtree(gt_move[gt]->b_target,gt_move[gt]->n_link,b_residual[gt],st->treelist->tree[gt]);
1905             
1906              dir_v1 = dir_v2 = dir_v0 = -1;
1907              For(j,3)
1908                {
1909                  if(gt_move[gt]->n_link->v[j] == gt_move[gt]->n_opp_to_link) dir_v0 = j;
1910                  else if(dir_v1 < 0)                                         dir_v1 = j;
1911                  else                                                        dir_v2 = j;
1912                }
1913             
1914              gt_move[gt]->n_link->b[dir_v0]->l->v = gt_move[gt]->l0;
1915                 
1916              if(gt_move[gt]->n_link->v[dir_v1]->num > gt_move[gt]->n_link->v[dir_v2]->num)
1917                {
1918                  gt_move[gt]->n_link->b[dir_v2]->l->v = gt_move[gt]->l1;
1919                  gt_move[gt]->n_link->b[dir_v1]->l->v = gt_move[gt]->l2;
1920                }
1921              else
1922                {
1923                  gt_move[gt]->n_link->b[dir_v1]->l->v = gt_move[gt]->l1;
1924                  gt_move[gt]->n_link->b[dir_v2]->l->v = gt_move[gt]->l2;
1925                }
1926            }
1927        }
1928    }
1929 
1930
1931  if(n_moves)
1932    {
1933      if(st_move->lnL > st->tree->best_lnL)
1934        {
1935          t_edge *st_target, *st_residual;
1936
1937          /* Apply the move on the super-tree */
1938          Prune_Subtree(st_move->n_link,
1939                        st_move->n_opp_to_link,
1940                        &st_target,
1941                        &st_residual,
1942                        st->tree);
1943         
1944          Graft_Subtree(st_move->b_target,
1945                        st_move->n_link,
1946                        st_residual,
1947                        st->tree);
1948         
1949                 
1950          /* Map gt and st nodes and edges */
1951          PART_Do_Mapping(st);
1952
1953          time(&(st->tree->t_current));
1954
1955          Set_Both_Sides(YES,st->tree);   
1956          st->tree->c_lnL      = PART_Lk(st);
1957          st->tree->c_pars     = PART_Pars(st);
1958         
1959         
1960          if(FABS(st->tree->c_lnL - st_move->lnL) > st->tree->mod->s_opt->min_diff_lk_local)
1961            {
1962              PhyML_Printf("\n. st->tree->c_lnL = %f st_move->lnL = %f\n",
1963                     st->tree->c_lnL,st_move->lnL);
1964
1965              For(gt,st->n_part)
1966                {
1967                  PhyML_Printf("\n. truth -> %f ; move -> %f",
1968                               Lk(NULL,st->treelist->tree[gt]),
1969                               gt_move[gt] ? gt_move[gt]->lnL : -1.);
1970                }
1971            }
1972         
1973          PhyML_Printf("\n. (%5d sec) [+ ] [%10.2f] [%5d] -- ",
1974                 (int)(st->tree->t_current - st->tree->t_beg),
1975                 st->tree->c_lnL,st->tree->c_pars);       
1976         
1977          For(gt,st->n_part)
1978            PhyML_Printf("[%10.2f] ",st->treelist->tree[gt]->c_lnL);
1979         
1980         
1981          st->tree->n_improvements++;
1982          st->tree->best_lnL = st->tree->c_lnL;
1983          For(gt,st->n_part) Record_Br_Len(st->treelist->tree[gt]);
1984         
1985          Free(init_target);
1986          Free(b_residual);
1987          Free(gt_move);
1988         
1989          return 1;
1990        }
1991/*       else */
1992/*      { */
1993/*        For(gt,st->n_part)  */
1994/*          { */
1995/*            if(gt_move[gt]) */
1996/*              { */
1997/*                Lk(st->treelist->tree[gt]); */
1998/*                Fast_Br_Len_Recur(st->treelist->tree[gt]->a_nodes[0], */
1999/*                                  st->treelist->tree[gt]->a_nodes[0]->v[0], */
2000/*                                  st->treelist->tree[gt]->a_nodes[0]->b[0], */
2001/*                                  st->treelist->tree[gt]); */
2002/*              } */
2003/*          } */
2004         
2005/*        time(&(st->tree->t_current)); */
2006/*        st->tree->both_sides = 1; */
2007/*        st->tree->c_lnL      = PART_Lk(st); */
2008         
2009/*        if(st->tree->c_lnL > st->tree->best_lnL) */
2010/*          { */
2011/*            t_edge *st_target, *st_residual; */
2012             
2013/*            /\* Apply the move on the super-tree *\/ */
2014/*            Prune_Subtree(st_move->n_link, */
2015/*                          st_move->n_opp_to_link,                          */
2016/*                          &st_target, */
2017/*                          &st_residual, */
2018/*                          st->tree); */
2019             
2020/*            Graft_Subtree(st_move->b_target, */
2021/*                          st_move->n_link, */
2022/*                          st_residual, */
2023/*                          st->tree); */
2024             
2025             
2026/*            /\* Map gt and st nodes and edges *\/ */
2027/*            PART_Do_Mapping(st); */
2028
2029
2030/*            st->tree->c_pars = PART_Pars(st); */
2031/*            PhyML_Printf("\n. (%5d sec) [++] [%10.2f] [%5d] -- ", */
2032/*                   (int)(st->tree->t_current-st->tree->t_beg), */
2033/*                   st->tree->c_lnL, */
2034/*                   st->tree->c_pars); */
2035/*            For(gt,st->n_part) */
2036/*              PhyML_Printf("[%10.2f] ",st->treelist->tree[gt]->c_lnL); */
2037             
2038/*            st->tree->n_improvements++; */
2039/*            st->tree->best_lnL = st->tree->c_lnL; */
2040/*            For(gt,st->n_part) Record_Br_Len(st->treelist->tree[gt]); */
2041
2042/*            Free(init_target); */
2043/*            Free(b_residual); */
2044/*            Free(gt_move); */
2045
2046/*            return 1; */
2047/*          } */
2048/*      } */
2049    }
2050 
2051  For(gt,st->n_part) 
2052    {
2053      if(gt_move[gt])
2054        {         
2055          /* Regraft the subtree at its original position */
2056          Prune_Subtree(gt_move[gt]->n_link,
2057                        gt_move[gt]->n_opp_to_link,
2058                        &(gt_move[gt]->b_target),
2059                        &(b_residual[gt]),
2060                        st->treelist->tree[gt]);
2061
2062          Graft_Subtree(init_target[gt],
2063                        gt_move[gt]->n_link,
2064                        b_residual[gt],
2065                        st->treelist->tree[gt]);         
2066
2067          /* Restore branch lengths */
2068          Restore_Br_Len(st->treelist->tree[gt]);
2069        }
2070    }
2071 
2072  Set_Both_Sides(YES,st->tree);
2073  st->tree->c_lnL      = PART_Lk(st);
2074  st->tree->c_pars     = PART_Pars(st);
2075
2076  time(&(st->tree->t_current));
2077 
2078  PhyML_Printf("\n. (%5d sec) [--] [%10.2f] [%5d] -- ",
2079         (int)(st->tree->t_current - st->tree->t_beg),
2080         st->tree->c_lnL,st->tree->c_pars);       
2081 
2082  For(gt,st->n_part) PhyML_Printf("[%10.2f] ",st->treelist->tree[gt]->c_lnL);
2083
2084  Free(init_target);
2085  Free(b_residual);
2086  Free(gt_move);
2087
2088  return 0;
2089
2090}
2091
2092//////////////////////////////////////////////////////////////
2093//////////////////////////////////////////////////////////////
2094
2095
2096void PART_NNI(t_edge *st_b, supert_tree *st)
2097{ 
2098  t_node *v1, *v2, *v3, *v4;
2099  phydbl lk0_opt, lk1_opt, lk2_opt;
2100  int i,j;
2101  phydbl *init_bl;
2102  t_edge **map_edge_bef_swap, **map_edge_aft_swap;
2103
2104
2105  init_bl = (phydbl *)mCalloc(st->n_bl_part,sizeof(phydbl));
2106  map_edge_bef_swap = (t_edge **)mCalloc(st->n_part,sizeof(t_edge *));
2107  map_edge_aft_swap = (t_edge **)mCalloc(st->n_part,sizeof(t_edge *));
2108
2109
2110  v1 = st_b->left->v[st_b->l_v1];
2111  v2 = st_b->left->v[st_b->l_v2];
2112  v3 = st_b->rght->v[st_b->r_v1];
2113  v4 = st_b->rght->v[st_b->r_v2];
2114
2115  lk0_opt  = lk1_opt  = lk2_opt  = UNLIKELY;
2116
2117  if(v1->num < v2->num)
2118    {
2119      Check_Dirs(st->tree);
2120      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2121      Warn_And_Exit("");
2122    }
2123
2124  if(v3->num < v4->num)
2125    {
2126      Check_Dirs(st->tree);
2127      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2128      Warn_And_Exit("");
2129    }
2130
2131 
2132/*   PhyML_Printf("oooooooo\n"); */
2133/*   Print_Node(st->tree->a_nodes[0], */
2134/*           st->tree->a_nodes[0]->v[0], */
2135/*           st->tree); */
2136/*   PhyML_Printf(">>>>>>>\n"); */
2137/*   For(i,st->n_part) */
2138/*     { */
2139/*       Print_Node(st->treelist->tree[i]->a_nodes[0], */
2140/*               st->treelist->tree[i]->a_nodes[0]->v[0], */
2141/*               st->treelist->tree[i]); */
2142/*       PhyML_Printf("<<<<<<<\n"); */
2143/*     } */
2144
2145 
2146  PART_Record_Br_Len(st);
2147
2148  For(i,st->n_part) map_edge_bef_swap[i] = NULL;
2149  For(i,st->n_part) if(st->map_st_edge_in_gt[i][st_b->num]) map_edge_bef_swap[i] = st->map_st_edge_in_gt[i][st_b->num];
2150
2151  /* First alternative topological configuration */
2152  /* Swap */
2153  PART_Swap(v2,st_b->left,st_b->rght,v3,st);
2154  Swap(v2,st_b->left,st_b->rght,v3,st->tree);
2155  PART_Do_Mapping(st);
2156  PART_Set_Bl(st->bl,st);
2157  For(i,st->n_part) map_edge_aft_swap[i] = NULL;
2158  For(i,st->n_part) if(st->map_st_edge_in_gt[i][st_b->num]) map_edge_aft_swap[i] = st->map_st_edge_in_gt[i][st_b->num];
2159  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2160  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2161  For(i,st->n_part) if(map_edge_bef_swap[i] && map_edge_aft_swap[i]) 
2162    {
2163      For(j,3) if((map_edge_aft_swap[i]->left->v[j]) && 
2164                  (map_edge_aft_swap[i]->left->b[j] == map_edge_bef_swap[i]) &&
2165                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2166        Update_P_Lk(st->treelist->tree[i],
2167                    map_edge_aft_swap[i],
2168                    map_edge_aft_swap[i]->left);
2169      For(j,3) if((map_edge_aft_swap[i]->rght->v[j]) && 
2170                  (map_edge_aft_swap[i]->rght->b[j] == map_edge_bef_swap[i]) &&
2171                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2172        Update_P_Lk(st->treelist->tree[i],
2173                    map_edge_aft_swap[i],
2174                    map_edge_aft_swap[i]->rght);
2175    }
2176  PART_Update_Lk_At_Given_Edge(st_b,st);
2177  lk1_opt  = PART_Br_Len_Brent(st_b,0,st);
2178  For(i,st->n_part) st->bl1[st->bl_partition[i]][st_b->num] = st->bl[st->bl_partition[i]][st_b->num];
2179  /* Unswap */
2180  PART_Swap(v3,st_b->left,st_b->rght,v2,st);
2181  Swap(v3,st_b->left,st_b->rght,v2,st->tree);
2182  PART_Do_Mapping(st);
2183  PART_Restore_Br_Len(st);
2184  PART_Set_Bl(st->bl,st);
2185  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2186  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2187  For(i,st->n_part) if(map_edge_bef_swap[i] && map_edge_aft_swap[i]) 
2188    {
2189      For(j,3) if((map_edge_aft_swap[i]->left->v[j]) && 
2190                  (map_edge_aft_swap[i]->left->b[j] == map_edge_bef_swap[i]) &&
2191                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2192        Update_P_Lk(st->treelist->tree[i],
2193                    map_edge_aft_swap[i],
2194                    map_edge_aft_swap[i]->left);
2195      For(j,3) if((map_edge_aft_swap[i]->rght->v[j]) && 
2196                  (map_edge_aft_swap[i]->rght->b[j] == map_edge_bef_swap[i]) &&
2197                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2198        Update_P_Lk(st->treelist->tree[i],
2199                    map_edge_aft_swap[i],
2200                    map_edge_aft_swap[i]->rght);
2201    }
2202
2203
2204
2205
2206  /* Second alternative topological configuration */
2207  /* Swap */
2208  PART_Swap(v2,st_b->left,st_b->rght,v4,st);
2209  Swap(v2,st_b->left,st_b->rght,v4,st->tree);
2210  PART_Do_Mapping(st);
2211  PART_Set_Bl(st->bl,st);
2212  For(i,st->n_part) map_edge_aft_swap[i] = NULL;
2213  For(i,st->n_part) if(st->map_st_edge_in_gt[i][st_b->num]) map_edge_aft_swap[i] = st->map_st_edge_in_gt[i][st_b->num];
2214  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2215  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2216  For(i,st->n_part) if(map_edge_bef_swap[i] && map_edge_aft_swap[i]) 
2217    {
2218      For(j,3) if((map_edge_aft_swap[i]->left->v[j]) && 
2219                  (map_edge_aft_swap[i]->left->b[j] == map_edge_bef_swap[i]) &&
2220                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2221        Update_P_Lk(st->treelist->tree[i],
2222                    map_edge_aft_swap[i],
2223                    map_edge_aft_swap[i]->left);
2224      For(j,3) if((map_edge_aft_swap[i]->rght->v[j]) && 
2225                  (map_edge_aft_swap[i]->rght->b[j] == map_edge_bef_swap[i]) &&
2226                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2227        Update_P_Lk(st->treelist->tree[i],
2228                    map_edge_aft_swap[i],
2229                    map_edge_aft_swap[i]->rght);
2230    }
2231
2232  PART_Update_Lk_At_Given_Edge(st_b,st);
2233  lk2_opt  = PART_Br_Len_Brent(st_b,0,st);
2234  For(i,st->n_part) st->bl2[st->bl_partition[i]][st_b->num] = st->bl[st->bl_partition[i]][st_b->num];
2235  /*   PhyML_Printf("\n. lk2_init = %f lk2_opt = %f",lk2_init,lk2_opt); */
2236  /* Unswap */
2237  PART_Swap(v4,st_b->left,st_b->rght,v2,st);
2238  Swap(v4,st_b->left,st_b->rght,v2,st->tree);
2239  PART_Do_Mapping(st);
2240  PART_Restore_Br_Len(st);
2241  PART_Set_Bl(st->bl,st);
2242  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2243  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2244  For(i,st->n_part) if(map_edge_bef_swap[i] && map_edge_aft_swap[i]) 
2245    {
2246      For(j,3) if((map_edge_aft_swap[i]->left->v[j]) && 
2247                  (map_edge_aft_swap[i]->left->b[j] == map_edge_bef_swap[i]) &&
2248                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2249        Update_P_Lk(st->treelist->tree[i],
2250                    map_edge_aft_swap[i],
2251                    map_edge_aft_swap[i]->left);
2252      For(j,3) if((map_edge_aft_swap[i]->rght->v[j]) && 
2253                  (map_edge_aft_swap[i]->rght->b[j] == map_edge_bef_swap[i]) &&
2254                  (map_edge_aft_swap[i] != map_edge_bef_swap[i])) 
2255        Update_P_Lk(st->treelist->tree[i],
2256                    map_edge_aft_swap[i],
2257                    map_edge_aft_swap[i]->rght);
2258    }
2259
2260
2261  /* Back to the initial topological configuration
2262   * and branch lengths.
2263   */
2264  PART_Do_Mapping(st);
2265  PART_Set_Bl(st->bl,st);
2266  PART_Restore_Br_Len(st);
2267  For(i,st->n_part) map_edge_aft_swap[i] = NULL;
2268  For(i,st->n_part) if(st->map_st_edge_in_gt[i][st_b->num]) map_edge_aft_swap[i] = st->map_st_edge_in_gt[i][st_b->num];
2269  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2270  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2271  PART_Update_Lk_At_Given_Edge(st_b,st);
2272  lk0_opt  = PART_Br_Len_Brent(st_b,0,st);
2273  For(i,st->n_part) st->bl0[st->bl_partition[i]][st_b->num] = st->bl[st->bl_partition[i]][st_b->num];
2274
2275  PART_Restore_Br_Len(st);
2276  PART_Set_Bl(st->bl,st);
2277  For(i,st->n_part) if(map_edge_bef_swap[i]) Update_PMat_At_Given_Edge(map_edge_bef_swap[i],st->treelist->tree[i]);
2278  For(i,st->n_part) if(map_edge_aft_swap[i]) Update_PMat_At_Given_Edge(map_edge_aft_swap[i],st->treelist->tree[i]);
2279  PART_Update_Lk_At_Given_Edge(st_b,st);
2280
2281
2282
2283/*   For(i,2*st->tree->n_otu-3) */
2284/*     PhyML_Printf("\n. 3 Edge %3d --> lnL=%f",i,PART_Lk_At_Given_Edge(st->tree->a_edges[i],st)); */
2285
2286
2287
2288  st_b->nni->lk0 = lk0_opt;
2289  st_b->nni->lk1 = lk1_opt;
2290  st_b->nni->lk2 = lk2_opt;
2291
2292  st_b->nni->score = lk0_opt - MAX(lk1_opt,lk2_opt);
2293
2294  if((st_b->nni->score <  st->tree->mod->s_opt->min_diff_lk_local) &&
2295     (st_b->nni->score > -st->tree->mod->s_opt->min_diff_lk_local))
2296    {
2297      st_b->nni->score = .0;
2298      st_b->nni->lk1   = st_b->nni->lk0;
2299      st_b->nni->lk2   = st_b->nni->lk0;
2300     }
2301
2302  PART_Restore_Br_Len(st);
2303  PART_Update_Lk_At_Given_Edge(st_b,st); /* to replace by PART_Update_PMat_At_Given_Edge(st_b,st); */
2304/*   PhyML_Printf("\n. lk_end = %f",st->tree->c_lnL); */
2305/*   For(i,2*st->tree->n_otu-3) PhyML_Printf("\n. %f",PART_Lk_At_Given_Edge(st->tree->a_edges[i],st)); */
2306/*   PhyML_Printf("\n. lk_end = %f",PART_Lk(st)); */
2307/*   PhyML_Printf("\n"); */
2308
2309/*   PhyML_Printf("\n. Edge %3d, score = %20f",st_b->num,st_b->nni->score); */
2310
2311  if(st_b->num == 90)
2312    PhyML_Printf("\n. v1=%d v2=%d v3=%d v4=%d left-%d right-%d",
2313           v1->num,
2314           v2->num,
2315           v3->num,
2316           v4->num,
2317           st_b->left->num,
2318           st_b->rght->num);
2319
2320
2321
2322  if(lk0_opt > MAX(lk1_opt,lk2_opt))
2323    {
2324      st_b->nni->best_conf = 0;
2325      st_b->nni->swap_node_v1 = NULL;
2326      st_b->nni->swap_node_v2 = NULL;
2327      st_b->nni->swap_node_v3 = NULL;
2328      st_b->nni->swap_node_v4 = NULL;
2329    }
2330  else if(lk1_opt > MAX(lk0_opt,lk2_opt))
2331    {
2332      st_b->nni->best_conf    = 1;
2333      st_b->nni->swap_node_v1 = v2;
2334      st_b->nni->swap_node_v2 = st_b->left;
2335      st_b->nni->swap_node_v3 = st_b->rght;
2336      st_b->nni->swap_node_v4 = v3;
2337    }
2338  else if(lk2_opt > MAX(lk0_opt,lk1_opt))
2339    {
2340      st_b->nni->best_conf    = 2;
2341      st_b->nni->swap_node_v1 = v2;
2342      st_b->nni->swap_node_v2 = st_b->left;
2343      st_b->nni->swap_node_v3 = st_b->rght;
2344      st_b->nni->swap_node_v4 = v4;
2345    }
2346  else
2347    {
2348      st_b->nni->score        = +1.0;
2349      st_b->nni->best_conf    = 0;
2350      st_b->nni->swap_node_v1 = NULL;
2351      st_b->nni->swap_node_v2 = NULL;
2352      st_b->nni->swap_node_v3 = NULL;
2353      st_b->nni->swap_node_v4 = NULL;
2354    }
2355 
2356  Free(init_bl);
2357  Free(map_edge_aft_swap);
2358  Free(map_edge_bef_swap);
2359}
2360
2361//////////////////////////////////////////////////////////////
2362//////////////////////////////////////////////////////////////
2363
2364
2365void PART_Swap(t_node *st_a, t_node *st_b, t_node *st_c, t_node *st_d, supert_tree *st)
2366{
2367  int i,j;
2368  t_node *gt_a, *gt_b, *gt_c, *gt_d;
2369  int ab, ba, cd, dc, bc;
2370
2371  ab = ba = cd = dc = bc = -1;
2372
2373  For(i,3) if(st_a->v[i] == st_b) { ab = i; break; }
2374  For(i,3) if(st_b->v[i] == st_a) { ba = i; break; }
2375  For(i,3) if(st_c->v[i] == st_d) { cd = i; break; }
2376  For(i,3) if(st_d->v[i] == st_c) { dc = i; break; }
2377  For(i,3) if(st_b->v[i] == st_c) { bc = i; break; }
2378
2379  if(ab < 0 || ba < 0 || cd < 0 || dc < 0)
2380    {
2381      PhyML_Printf("\n. Nodes %d %d %d %d\n",st_a->num,st_b->num,st_c->num,st_d->num);
2382      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2383      Warn_And_Exit("");
2384    }
2385
2386  gt_a = gt_b = gt_c = gt_d = NULL;
2387 
2388  For(i,st->n_part)
2389    {
2390      gt_b = st->match_st_node_in_gt[i][st_b->num];
2391      gt_c = st->match_st_node_in_gt[i][st_c->num];
2392     
2393      if(gt_b && gt_c) /* The st t_edge with st_b and st_c at its extremities
2394                        * matches an t_edge in gt
2395                        */
2396        {
2397#ifdef DEBUG
2398          For(j,3) if((gt_b->v[j]) && (gt_b->v[j] == gt_c)) break;
2399          if(j == 3)
2400            {
2401              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2402              Warn_And_Exit("");
2403            }
2404#endif
2405          gt_a = st->map_st_node_in_gt[i][st_a->num][ab][0];
2406          gt_d = st->map_st_node_in_gt[i][st_d->num][dc][0];
2407          Swap(gt_a,gt_b,gt_c,gt_d,st->treelist->tree[i]);
2408        }
2409    }
2410}
2411
2412//////////////////////////////////////////////////////////////
2413//////////////////////////////////////////////////////////////
2414
2415
2416void PART_Set_Bl(phydbl **bl, supert_tree *st)
2417{
2418  int i,j;
2419  t_edge *gt_b;
2420
2421  gt_b = NULL;
2422                                               
2423  /* Set all the actual branch lengths to 0.0
2424   */
2425  For(i,st->n_part)
2426    {
2427      For(j,2*st->treelist->tree[i]->n_otu-3)
2428        {
2429          gt_b = st->treelist->tree[i]->a_edges[j];
2430          gt_b->l->v = .0;
2431        }
2432    }
2433
2434  /* Update every branch length
2435   */ 
2436  For(i,2*st->tree->n_otu-3)
2437    {
2438      For(j,st->n_part)
2439        {
2440          gt_b = st->map_st_edge_in_gt[j][i];   
2441         
2442          /* Need to make sure that st->tree->a_edges[i] is on an existing path in gt */
2443          if((st->map_st_node_in_gt[j][st->tree->a_edges[i]->left->num][st->tree->a_edges[i]->l_r][0]) &&
2444             (st->map_st_node_in_gt[j][st->tree->a_edges[i]->rght->num][st->tree->a_edges[i]->r_l][0]))
2445            {
2446              gt_b->l->v += bl[st->bl_partition[j]][i];
2447            }
2448        }
2449    }
2450}
2451
2452//////////////////////////////////////////////////////////////
2453//////////////////////////////////////////////////////////////
2454
2455
2456void PART_Record_Br_Len(supert_tree *st)
2457{
2458  int i,j;
2459  For(i,st->n_part) For(j,2*st->tree->n_otu-3) st->bl_cpy[i][j] = st->bl[i][j];
2460}
2461
2462//////////////////////////////////////////////////////////////
2463//////////////////////////////////////////////////////////////
2464
2465
2466void PART_Restore_Br_Len(supert_tree *st)
2467{
2468  int i,j;
2469  For(i,st->n_part) For(j,2*st->tree->n_otu-3) st->bl[i][j] = st->bl_cpy[i][j];
2470}
2471
2472//////////////////////////////////////////////////////////////
2473//////////////////////////////////////////////////////////////
2474
2475
2476phydbl PART_Lk(supert_tree *st)
2477{
2478  int i;
2479
2480  PART_Do_Mapping(st);
2481  PART_Set_Bl(st->bl,st); 
2482
2483  st->tree->c_lnL = .0;
2484  For(i,st->n_part) 
2485    {
2486      Set_Both_Sides(YES,st->treelist->tree[i]);         
2487      Lk(NULL,st->treelist->tree[i]);
2488/*       PhyML_Printf("\n. Tree %3d lnL = %f",i+1,st->treelist->tree[i]->c_lnL); */
2489      st->tree->c_lnL += st->treelist->tree[i]->c_lnL;
2490    }
2491  return st->tree->c_lnL;
2492}
2493
2494//////////////////////////////////////////////////////////////
2495//////////////////////////////////////////////////////////////
2496
2497
2498phydbl PART_Lk_At_Given_Edge(t_edge *st_b, supert_tree *st)
2499{
2500  int i;
2501  t_edge *gt_b;
2502  phydbl lnL;
2503
2504  PART_Set_Bl(st->bl,st);
2505
2506  gt_b = NULL;
2507  st->tree->c_lnL = .0;
2508  lnL = .0;
2509  For(i,st->n_part)
2510    {     
2511      gt_b = st->map_st_edge_in_gt[i][st_b->num];
2512      lnL = Lk(gt_b,st->treelist->tree[i]);
2513      st->tree->c_lnL += lnL;
2514/*       PhyML_Printf("\n. gt %d st t_edge %d gt t_edge %d lnL=%f l=%f ",i,st_b->num,gt_b->num,lnL,gt_b->l->v); */
2515    }
2516  return st->tree->c_lnL;
2517}
2518
2519//////////////////////////////////////////////////////////////
2520//////////////////////////////////////////////////////////////
2521
2522
2523phydbl PART_Update_Lk_At_Given_Edge(t_edge *st_b, supert_tree *st)
2524{
2525  int i;
2526  t_edge *gt_b;
2527
2528  PART_Set_Bl(st->bl,st);
2529 
2530  gt_b = NULL;
2531  st->tree->c_lnL = .0;
2532  For(i,st->n_part)
2533    {
2534      gt_b = st->map_st_edge_in_gt[i][st_b->num];
2535      if(gt_b) st->tree->c_lnL += Update_Lk_At_Given_Edge(gt_b,st->treelist->tree[i]);
2536      else     st->tree->c_lnL += st->treelist->tree[i]->c_lnL;
2537    }
2538  return st->tree->c_lnL;
2539}
2540
2541
2542//////////////////////////////////////////////////////////////
2543//////////////////////////////////////////////////////////////
2544
2545
2546void PART_Fill_Model_Partitions_Table(supert_tree *st)
2547{
2548  int i,j;
2549  char *c;
2550  char *abc;
2551  int lig, col;
2552  int n_groups;
2553  int *encountered_vals;
2554
2555  c = (char *)mCalloc(10,sizeof(char));
2556  abc = (char *)mCalloc(20,sizeof(char));
2557  encountered_vals = (int *)mCalloc(st->n_part,sizeof(int));
2558
2559  strcpy(abc,"ABCDEFGHIJKLMNOP\0");
2560
2561  PhyML_Printf("\n\n\n");
2562  lig = col = 0;
2563  while(1)
2564    {
2565      PhyML_Printf("\n\n");
2566      For(i,st->n_part)
2567        PhyML_Printf(". Data set %3d : %s\n",i+1,st->optionlist[i]->in_align_file);
2568
2569      PhyML_Printf("\n. Data set             ");
2570      For(i,st->n_part) PhyML_Printf("%3d ",i+1);
2571      PhyML_Printf("\n. -A- t_edge lengths     ");
2572      For(i,st->n_part) PhyML_Printf("%3d ",st->bl_partition[i]);
2573
2574      if(lig == 1) break;
2575
2576      PhyML_Printf("\n. (%c-%2d)> ",abc[lig],col+1);
2577      Getstring_Stdin(c);
2578     
2579      switch(lig)
2580        {
2581        case 0 :
2582          {
2583            st->bl_partition[col] = atoi(c);
2584            break;
2585          }
2586        default :
2587          {
2588            break;
2589          }
2590        }
2591
2592      col++;
2593
2594      if(col == st->n_part)
2595        {
2596          col = 0;
2597          lig++;
2598        }
2599    }
2600   
2601  n_groups = 0;
2602  For(i,st->n_part) 
2603    {
2604      For(j,n_groups)
2605        if(encountered_vals[j] == st->bl_partition[i])
2606          break;
2607
2608      if(j == n_groups) 
2609        {
2610          encountered_vals[n_groups] = st->bl_partition[i];
2611          n_groups++;
2612        }
2613      st->bl_partition[i] = j;
2614    }
2615 
2616  st->n_bl_part = n_groups;
2617
2618  Free(encountered_vals);
2619  Free(c);
2620  Free(abc);
2621}
2622
2623//////////////////////////////////////////////////////////////
2624//////////////////////////////////////////////////////////////
2625
2626
2627phydbl PART_Br_Len_Brent(t_edge *st_b, int quickdirty, supert_tree *st)
2628{
2629  phydbl ax,  cx;
2630  int part;
2631  phydbl cur_l;
2632
2633
2634  For(part,st->n_bl_part)
2635    {
2636      cur_l = st->bl[part][st_b->num];
2637     
2638      ax = 10.*cur_l;
2639      cx = st->tree->mod->l_min;
2640
2641      Generic_Brent_Lk(&(st->bl[part][st_b->num]),
2642                       ax,cx,
2643                       st->tree->mod->s_opt->min_diff_lk_local,
2644                       st->tree->mod->s_opt->brent_it_max,
2645                       st->tree->mod->s_opt->quickdirty,
2646                       Wrap_Part_Lk_At_Given_Edge,st_b,NULL,st,NO);
2647     
2648    }
2649  return st->tree->c_lnL;
2650}
2651
2652//////////////////////////////////////////////////////////////
2653//////////////////////////////////////////////////////////////
2654
2655
2656void PART_Initialise_Bl_Partition(supert_tree *st)
2657{
2658  int i,j;
2659 
2660  For(i,st->n_bl_part)
2661    {
2662      For(j,2*st->tree->n_otu-3)
2663        {
2664          st->bl[i][j] = .1;
2665        }
2666    }
2667}
2668
2669//////////////////////////////////////////////////////////////
2670//////////////////////////////////////////////////////////////
2671
2672
2673void PART_Optimize_Br_Len_Serie(t_node *st_a, t_node *st_d, t_edge *st_b, supert_tree *st)
2674{
2675  phydbl lk_init;
2676  int i;
2677
2678  lk_init = st->tree->c_lnL;
2679 
2680  PART_Br_Len_Brent(st_b,0,st);
2681
2682  if(st->tree->c_lnL < lk_init - st->tree->mod->s_opt->min_diff_lk_local)
2683    { 
2684      PhyML_Printf("\n== %f -- %f",lk_init,st->tree->c_lnL);
2685      PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
2686      Warn_And_Exit("");
2687    }
2688   
2689  if(st_d->tax) return;
2690  else For(i,3) if(st_d->v[i] != st_a)
2691    {
2692      PART_Update_P_Lk(st_d->b[i],st_d,st);
2693      PART_Optimize_Br_Len_Serie(st_d,st_d->v[i],st_d->b[i],st);
2694    }
2695
2696  For(i,3) if((st_d->v[i] == st_a) && (!st_d->v[i]->tax)) PART_Update_P_Lk(st_d->b[i],st_d,st);
2697
2698}
2699
2700//////////////////////////////////////////////////////////////
2701//////////////////////////////////////////////////////////////
2702
2703
2704void PART_Update_P_Lk(t_edge *st_b, t_node *st_n, supert_tree *st)
2705{
2706  int i,dir;
2707
2708  dir = -1;
2709  For(i,3) if((st_n->b[i]) && (st_n->b[i] == st_b)) {dir = i; break;}
2710 
2711  if(dir < 0)
2712    {
2713      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2714      Warn_And_Exit("");
2715    }
2716
2717  For(i,st->n_part)
2718    {
2719      if((st->map_st_node_in_gt[i][st_n->num][dir][0]) && (!st->map_st_node_in_gt[i][st_n->num][dir][0]->tax))
2720        {         
2721          Update_P_Lk(st->treelist->tree[i],
2722                      st->map_st_edge_in_gt[i][st_b->num],
2723                      st->map_st_node_in_gt[i][st_n->num][dir][0]);
2724        }
2725    }
2726}
2727
2728//////////////////////////////////////////////////////////////
2729//////////////////////////////////////////////////////////////
2730
2731
2732void PART_Make_N_Swap(t_edge **st_b, int beg, int end, supert_tree *st)
2733{
2734  int i;
2735  int dim;
2736
2737  dim = 2*st->tree->n_otu-2;
2738
2739  st->tree->n_swap = 0;
2740  for(i=beg;i<end;i++)
2741    {
2742      if(st_b[i]->left->tax || st_b[i]->rght->tax)
2743        {
2744          PhyML_Printf("\n. Edge %d is external.",st_b[i]->num);
2745          PhyML_Printf("\n. Err in file %s at line %d\n\n.",__FILE__,__LINE__);
2746          Warn_And_Exit("");
2747        }
2748     
2749      PART_Swap(st_b[i]->nni->swap_node_v2->v[st->tree->t_dir[st_b[i]->nni->swap_node_v2->num*dim+st_b[i]->nni->swap_node_v1->num]],
2750              st_b[i]->nni->swap_node_v2,
2751              st_b[i]->nni->swap_node_v3,
2752              st_b[i]->nni->swap_node_v3->v[st->tree->t_dir[st_b[i]->nni->swap_node_v3->num*dim+st_b[i]->nni->swap_node_v4->num]],
2753              st);
2754
2755      Swap(st_b[i]->nni->swap_node_v2->v[st->tree->t_dir[st_b[i]->nni->swap_node_v2->num*dim+st_b[i]->nni->swap_node_v1->num]],
2756           st_b[i]->nni->swap_node_v2,
2757           st_b[i]->nni->swap_node_v3,
2758           st_b[i]->nni->swap_node_v3->v[st->tree->t_dir[st_b[i]->nni->swap_node_v3->num*dim+st_b[i]->nni->swap_node_v4->num]],
2759           st->tree);
2760
2761      PART_Do_Mapping(st);
2762
2763      st->tree->n_swap++;
2764    }
2765}
2766
2767//////////////////////////////////////////////////////////////
2768//////////////////////////////////////////////////////////////
2769
2770
2771void PART_Update_Bl(phydbl fact, supert_tree *st)
2772{
2773  int i,j;
2774 
2775  For(i,2*st->tree->n_otu-3)
2776    {
2777      For(j,st->n_part)
2778        st->bl[st->bl_partition[j]][i] = 
2779        st->bl_cpy[st->bl_partition[j]][i] + 
2780        (st->bl0[st->bl_partition[j]][i] - st->bl_cpy[st->bl_partition[j]][i]) * fact;
2781    }
2782  PART_Set_Bl(st->bl,st);
2783}
2784
2785//////////////////////////////////////////////////////////////
2786//////////////////////////////////////////////////////////////
2787
2788
2789void PART_Update_Bl_Swaped(t_edge **st_b, int n, supert_tree *st)
2790{
2791  int i,j;
2792 
2793  For(i,n)
2794    {
2795      For(j,st->n_part)
2796        {
2797          st->bl[st->bl_partition[j]][st_b[i]->num] = 
2798            (st_b[i]->nni->best_conf == 1)?
2799            (st->bl1[st->bl_partition[j]][st_b[i]->num]):
2800            (st->bl2[st->bl_partition[j]][st_b[i]->num]);
2801        }
2802    }
2803  PART_Set_Bl(st->bl,st);
2804}
2805
2806//////////////////////////////////////////////////////////////
2807//////////////////////////////////////////////////////////////
2808
2809
2810void PART_Do_Mapping(supert_tree *st)
2811{
2812  int k;
2813
2814  Fill_Dir_Table(st->tree);
2815  For(k,st->n_part)
2816    {
2817      Fill_Dir_Table(st->treelist->tree[k]);
2818      PART_Match_St_Nodes_In_Gt(st->treelist->tree[k],st);           
2819      PART_Match_St_Edges_In_Gt(st->treelist->tree[k],st);           
2820      PART_Map_St_Nodes_In_Gt(st->treelist->tree[k],st);
2821      PART_Map_St_Edges_In_Gt(st->treelist->tree[k],st);
2822      PART_Map_Gt_Edges_In_St(st->treelist->tree[k],st);
2823    }
2824}
2825
2826//////////////////////////////////////////////////////////////
2827//////////////////////////////////////////////////////////////
2828
2829
2830void PART_Print_Bl(supert_tree *st)
2831{
2832  int i,j;
2833 
2834  For(j,2*st->tree->n_otu-3)
2835    { 
2836      PhyML_Printf("\n. t_edge %4d ",j);
2837      For(i,st->n_bl_part)
2838        {
2839          PhyML_Printf("%f ",st->bl[i][j]);
2840        }
2841    }
2842}
2843
2844//////////////////////////////////////////////////////////////
2845//////////////////////////////////////////////////////////////
2846
2847//////////////////////////////////////////////////////////////
2848//////////////////////////////////////////////////////////////
2849
Note: See TracBrowser for help on using the repository browser.