source: branches/lib/GDE/PHYML20130708/phyml/src/m4.c

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

added most recent version of phyml
File size: 47.3 KB
Line 
1/*
2
3PhyML:  a program that  computes maximum likelihood phylogenies from
4DNA or AA homoLOGous sequences.
5
6Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8All parts of the source except where indicated are distributed under
9the GNU public licence. See http://www.opensource.org for details.
10
11*/
12
13
14/* Routines for Markov-Modulated Markov Models (M4) */
15
16#include "m4.h"
17
18int M4_main(int argc, char **argv)
19{
20 
21  calign *cdata;
22  option *io;
23  t_tree *tree;
24  int num_data_set;
25  int num_tree,num_rand_tree;
26  t_mod *mod;
27  time_t t_beg,t_end;
28  phydbl best_lnL;
29  int r_seed;
30  char *most_likely_tree=NULL;
31
32 
33#ifdef MPI
34  int rc;
35  rc = MPI_Init(&argc,&argv);
36  if (rc != MPI_SUCCESS) {
37    PhyML_Printf("\n. Error starting MPI program. Terminating.\n");
38    MPI_Abort(MPI_COMM_WORLD, rc);
39  }
40  MPI_Comm_size(MPI_COMM_WORLD,&Global_numTask);
41  MPI_Comm_rank(MPI_COMM_WORLD,&Global_myRank);
42#endif
43
44#ifdef QUIET
45  setvbuf(stdout,NULL,_IOFBF,2048);
46#endif
47
48
49  tree             = NULL;
50  mod              = NULL;
51  best_lnL         = UNLIKELY;
52     
53  io = (option *)Get_Input(argc,argv);
54  r_seed = (io->r_seed < 0)?(time(NULL)):(io->r_seed);
55  srand(r_seed);
56  io->r_seed = r_seed;
57
58
59  if(io->in_tree == 2) Test_Multiple_Data_Set_Format(io);
60  else io->n_trees = 1;
61
62
63  if((io->n_data_sets > 1) && (io->n_trees > 1))
64    {
65      io->n_data_sets = MIN(io->n_trees,io->n_data_sets);
66      io->n_trees     = MIN(io->n_trees,io->n_data_sets);
67    }
68
69  For(num_data_set,io->n_data_sets)
70    {
71      best_lnL = UNLIKELY;
72      Get_Seq(io);
73      Make_Model_Complete(io->mod);
74      Set_Model_Name(io->mod);
75      Print_Settings(io);
76      mod = io->mod;
77     
78      if(io->data)
79        {
80          if(io->n_data_sets > 1) PhyML_Printf("\n. Data set [#%d]\n",num_data_set+1);
81          cdata = Compact_Data(io->data,io);
82
83          Free_Seq(io->data,cdata->n_otu);
84         
85          if(cdata) Check_Ambiguities(cdata,io->datatype,io->state_len);
86          else
87            {
88              PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
89              Warn_And_Exit("");
90            }
91
92          for(num_tree=(io->n_trees == 1)?(0):(num_data_set);num_tree < io->n_trees;num_tree++)
93            {
94              if(!io->mod->s_opt->random_input_tree) io->mod->s_opt->n_rand_starts = 1;
95
96              For(num_rand_tree,io->mod->s_opt->n_rand_starts)
97                {
98                  if((io->mod->s_opt->random_input_tree) && (io->mod->s_opt->topo_search != NNI_MOVE))
99                    if(!io->quiet) PhyML_Printf("\n. [Random start %3d/%3d]\n",num_rand_tree+1,io->mod->s_opt->n_rand_starts);
100
101                  Init_Model(cdata,mod,io);
102
103                  if(io->mod->use_m4mod) M4_Init_Model(mod->m4mod,cdata,mod);
104
105                  switch(io->in_tree)
106                    {
107                    case 0 : case 1 : { tree = Dist_And_BioNJ(cdata,mod,io); break; }
108                    case 2 :          { tree = Read_User_Tree(cdata,mod,io); break; }
109                    }
110
111
112                  if(!tree) continue;
113
114                  time(&t_beg);
115                  time(&(tree->t_beg));
116     
117
118                  tree->mod         = mod;
119                  tree->io          = io;
120                  tree->data        = cdata;
121                  tree->n_pattern   = tree->data->crunch_len;
122
123                  Set_Both_Sides(YES,tree);
124   
125                  if(mod->s_opt->random_input_tree) Random_Tree(tree);
126
127                  if((!num_data_set) && (!num_tree) && (!num_rand_tree)) Check_Memory_Amount(tree);
128
129                  Prepare_Tree_For_Lk(tree);
130
131                  if(io->in_tree == 1) Spr_Pars(tree);
132                 
133                  if(io->do_alias_subpatt)
134                    {
135                      MIXT_Set_Alias_Subpatt(YES,tree);
136                      Lk(NULL,tree);
137                      MIXT_Set_Alias_Subpatt(NO,tree);
138                    }
139
140                  if(tree->mod->s_opt->opt_topo)
141                    {
142                      if(tree->mod->s_opt->topo_search      == NNI_MOVE) Simu_Loop(tree);
143                      else if(tree->mod->s_opt->topo_search == SPR_MOVE) Speed_Spr_Loop(tree);
144                      else                                               Best_Of_NNI_And_SPR(tree);
145                    }
146                  else
147                    {
148                      if(tree->mod->s_opt->opt_subst_param || 
149                         tree->mod->s_opt->opt_bl)                       Round_Optimize(tree,tree->data,ROUND_MAX);
150                      else                                               Lk(NULL,tree);
151                    }
152                 
153                 
154                  Set_Both_Sides(YES,tree);     
155                  Lk(NULL,tree);
156                  Pars(NULL,tree);
157                  Get_Tree_Size(tree);
158                  PhyML_Printf("\n. Log likelihood of the current tree: %f.\n",tree->c_lnL);
159
160                  Exit("\n");
161
162                  /* */
163                  M4_Compute_Proba_Hidden_States_On_Edges(tree);
164                  /* */
165
166                  Get_Best_Root_Position(tree);
167
168                  /* Print the tree estimated using the current random (or BioNJ) starting tree */
169                  if(io->mod->s_opt->n_rand_starts > 1)
170                    {
171                      Br_Len_Involving_Invar(tree);
172                      Print_Tree(io->fp_out_trees,tree);
173                      fflush(NULL);
174                    }
175
176                  /* Record the most likely tree in a string of characters */
177                  if(tree->c_lnL > best_lnL)
178                    {
179                      best_lnL = tree->c_lnL;
180                      Br_Len_Involving_Invar(tree);
181                      if(most_likely_tree) Free(most_likely_tree);
182                      most_likely_tree = Write_Tree(tree,NO);
183                      Get_Tree_Size(tree);
184                    }
185
186/*                JF(tree); */
187
188                  time(&t_end);
189                  Print_Fp_Out(io->fp_out_stats,t_beg,t_end,tree,
190                               io,num_data_set+1,
191                               (tree->mod->s_opt->n_rand_starts > 1)?
192                               (num_rand_tree):(num_tree),YES);
193                 
194                  if(tree->io->print_site_lnl) Print_Site_Lk(tree,io->fp_out_lk);
195
196                  /* Start from BioNJ tree */
197                  if((num_rand_tree == io->mod->s_opt->n_rand_starts-1) && (tree->mod->s_opt->random_input_tree))
198                    {
199                      /* Do one more iteration in the loop, but don't randomize the tree */
200                      num_rand_tree--;
201                      tree->mod->s_opt->random_input_tree = 0;
202                    }
203                 
204                  Free_Spr_List(tree);
205                  Free_One_Spr(tree->best_spr);
206                  if(tree->mat) Free_Mat(tree->mat);
207                  Free_Triplet(tree->triplet_struct);
208                  Free_Tree_Pars(tree);
209                  Free_Tree_Lk(tree);
210                  Free_Tree(tree);
211                }
212
213
214              /* Launch bootstrap analysis */
215              if(mod->bootstrap) 
216                {
217                  if(!io->quiet) PhyML_Printf("\n. Launch bootstrap analysis on the most likely tree...\n");
218
219                  #ifdef MPI
220                  MPI_Bcast (most_likely_tree, strlen(most_likely_tree)+1, MPI_CHAR, 0, MPI_COMM_WORLD);
221                  if(!io->quiet)  PhyML_Printf("\n. The bootstrap analysis will use %d CPUs.\n",Global_numTask);
222                  #endif
223
224                  most_likely_tree = Bootstrap_From_String(most_likely_tree,cdata,mod,io);
225                }
226              else if(io->ratio_test) 
227                {
228                  /* Launch aLRT */
229                  if(!io->quiet) PhyML_Printf("\n. Compute aLRT branch supports on the most likely tree...\n");
230                  most_likely_tree = aLRT_From_String(most_likely_tree,cdata,mod,io);
231                }
232
233              /* Print the most likely tree in the output file */
234              if(!io->quiet) PhyML_Printf("\n. Printing the most likely tree in file '%s'...\n", Basename(io->out_tree_file));
235              if(io->n_data_sets == 1) rewind(io->fp_out_tree);
236              PhyML_Fprintf(io->fp_out_tree,"%s\n",most_likely_tree);
237             
238
239              if(io->n_trees > 1 && io->n_data_sets > 1) break;
240            }
241          Free_Cseq(cdata);
242        }
243      else
244        {
245          PhyML_Printf("\n. No data was found.\n");
246          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
247          Warn_And_Exit("");
248        }
249      Free_Model_Complete(mod);
250    }
251 
252  if(most_likely_tree) Free(most_likely_tree);
253
254  if(mod->s_opt->n_rand_starts > 1) PhyML_Printf("\n. Best log likelihood: %f\n",best_lnL);
255
256  Free_Optimiz(mod->s_opt);
257  Free_Model_Basic(mod);
258
259  if(io->fp_in_align)  fclose(io->fp_in_align);
260  if(io->fp_in_tree)   fclose(io->fp_in_tree);
261  if(io->fp_out_lk)    fclose(io->fp_out_lk);
262  if(io->fp_out_tree)  fclose(io->fp_out_tree);
263  if(io->fp_out_trees) fclose(io->fp_out_trees);
264  if(io->fp_out_stats) fclose(io->fp_out_stats);
265
266  Free_Input(io);
267
268  time(&t_end);
269  Print_Time_Info(t_beg,t_end);
270
271#ifdef MPI
272  MPI_Finalize();
273#endif
274
275  return 0;
276}
277
278//////////////////////////////////////////////////////////////
279//////////////////////////////////////////////////////////////
280
281
282/* Allocate memory */
283m4 *M4_Make_Light()
284{
285  m4 *m4mod;
286
287  m4mod = (m4 *)mCalloc(1,sizeof(m4));
288  M4_Set_M4mod_Default(m4mod);
289  return m4mod;
290}
291
292//////////////////////////////////////////////////////////////
293//////////////////////////////////////////////////////////////
294
295
296void M4_Set_M4mod_Default(m4 *m4mod)
297{
298  m4mod->use_cov_alpha = 1;
299  m4mod->use_cov_alpha = 0;
300  m4mod->n_h           = 3;
301  m4mod->n_o           = 4;
302}
303
304//////////////////////////////////////////////////////////////
305//////////////////////////////////////////////////////////////
306
307/* Allocate memory */
308void M4_Make_Complete(int n_h, int n_o, m4 *m4mod)
309{
310  int i;
311
312  m4mod->n_h = n_h;
313  m4mod->n_o = n_o;
314  m4mod->n_o = n_o;
315  m4mod->o_rr = (phydbl *)mCalloc(n_o*n_o,sizeof(phydbl));
316  m4mod->o_fq = (phydbl *)mCalloc(n_o,sizeof(phydbl));
317  m4mod->o_mats = (phydbl **)mCalloc(n_h,sizeof(phydbl *));
318  For(i,n_h) m4mod->o_mats[i] = (phydbl *)mCalloc(n_o*n_o,sizeof(phydbl));
319  m4mod->h_mat = (phydbl *)mCalloc(n_h*n_h,sizeof(phydbl));
320  m4mod->h_rr = (phydbl *)mCalloc(n_h*n_h,sizeof(phydbl));
321  m4mod->h_fq = (phydbl *)mCalloc(n_h,sizeof(phydbl));
322  m4mod->multipl = (phydbl *)mCalloc(n_h,sizeof(phydbl));
323  m4mod->multipl_unscaled = (phydbl *)mCalloc(n_h,sizeof(phydbl));
324  m4mod->h_fq_unscaled = (phydbl *)mCalloc(n_h,sizeof(phydbl));
325}
326
327//////////////////////////////////////////////////////////////
328//////////////////////////////////////////////////////////////
329
330
331/* Fill in the (big) rate matrix of the M4 t_mod */ 
332void M4_Update_Qmat(m4 *m4mod, t_mod *mod)
333{
334  int i,j;
335  int n_s, n_o, n_h;
336  phydbl mr, sum;
337
338  /* The number of states in M4 models is the product
339     of the number of hidden states (or classes) by the
340     number of observable states
341   */
342  n_s = mod->ns;
343  n_o = m4mod->n_o;
344  n_h = m4mod->n_h;
345 
346  /* Set the relative substitution rates */
347  if(mod->m4mod->use_cov_alpha)
348    {
349      DiscreteGamma(m4mod->h_fq,m4mod->multipl,m4mod->alpha,m4mod->alpha,m4mod->n_h,mod->ras->gamma_median);
350    }
351  else if(mod->m4mod->use_cov_free)
352    {
353      sum = .0;
354      For(i,mod->m4mod->n_h) sum += FABS(mod->m4mod->h_fq_unscaled[i]);
355      For(i,mod->m4mod->n_h) mod->m4mod->h_fq[i] = FABS(mod->m4mod->h_fq_unscaled[i])/sum;
356     
357      do
358        {
359          sum = .0;
360          For(i,mod->m4mod->n_h)
361            {
362              if(mod->m4mod->h_fq[i] < 0.01) mod->m4mod->h_fq[i]=0.01;
363              if(mod->m4mod->h_fq[i] > 0.99) mod->m4mod->h_fq[i]=0.99;
364              sum += mod->m4mod->h_fq[i];
365            }
366
367          For(i,mod->m4mod->n_h) mod->m4mod->h_fq[i]/=sum;
368        }
369      while((sum > 1.01) || (sum < 0.99));
370
371
372      /* Make sure the multipliers are centered on 1.0 */
373      sum = .0;
374      For(i,mod->m4mod->n_h) sum += FABS(mod->m4mod->multipl_unscaled[i]) * mod->m4mod->h_fq[i];
375      For(i,mod->m4mod->n_h) mod->m4mod->multipl[i] = mod->m4mod->multipl_unscaled[i] / sum;
376     
377      /* printf("\n. WARNING\n"); */
378      /* mod->m4mod->h_fq[0] = 1./3; */
379      /* mod->m4mod->h_fq[1] = 1./3; */
380      /* mod->m4mod->h_fq[2] = 1./3; */
381
382      /* mod->m4mod->multipl[0] = 1.0; */
383      /* mod->m4mod->multipl[1] = 1.0; */
384      /* mod->m4mod->multipl[2] = 1.0; */
385
386      sum = 0;
387      For(i,mod->m4mod->n_h) sum += mod->m4mod->multipl[i] * mod->m4mod->h_fq[i];
388      if(sum < 0.99 || sum > 1.01)
389        {
390          PhyML_Printf("\n. sum = %f",sum);
391          PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
392          Warn_And_Exit("\n");
393        }
394
395      /* PhyML_Printf("\n__ "); */
396      /* For(i,mod->m4mod->n_h) PhyML_Printf("\n.%f %f %f", */
397      /*                                    mod->m4mod->h_fq[i], */
398      /*                                    mod->m4mod->h_fq_unscaled[i], */
399      /*                                    mod->m4mod->multipl[i]); */
400    }
401
402
403  /* PhyML_Printf("\n."); */
404  /* PhyML_Printf("\n. M4 model parameters"); */
405  /* m4mod->delta=.0; */
406  /* PhyML_Printf("\n. Delta = %f",m4mod->delta); */
407  /* For(i,mod->m4mod->n_h) PhyML_Printf("\n. multipl %d = %f",i,m4mod->multipl[i]); */
408  /* For(i,mod->m4mod->n_h) PhyML_Printf("\n. fq %d = %f",i,m4mod->h_fq[i]); */
409
410
411  /* Set up the stationary frequency vector */
412  For(i,n_s) mod->e_frq->pi->v[i] = m4mod->o_fq[i%n_o] * m4mod->h_fq[i/n_o];
413
414
415  if(mod->whichmodel != CUSTOM &&
416     mod->whichmodel != GTR    &&
417     mod->io->datatype == NT)   
418    {
419      phydbl kappa1,kappa2;
420
421      if((mod->whichmodel != F84) && (mod->whichmodel != TN93)) mod->lambda->v = 1.; 
422      else if(mod->whichmodel == F84)
423        {
424          mod->lambda->v = Get_Lambda_F84(mod->e_frq->pi->v,&(mod->kappa->v));
425        }
426
427      kappa2 = mod->kappa->v*2./(1.+mod->lambda->v);
428      kappa1 = kappa2 * mod->lambda->v;
429
430      /* A <-> C */ m4mod->o_rr[0] = 1.0;
431      /* A <-> G */ m4mod->o_rr[1] = kappa2;
432      /* A <-> T */ m4mod->o_rr[2] = 1.0;
433      /* C <-> G */ m4mod->o_rr[3] = 1.0;
434      /* C <-> T */ m4mod->o_rr[4] = kappa1;
435    }
436
437  /* Fill in the matrices of nucleotide or amino-acid substitution rates here */
438  Update_Qmat_Generic(m4mod->o_rr, m4mod->o_fq, m4mod->n_o, m4mod->o_mats[0]);
439
440  /* Print_Square_Matrix_Generic(n_o,m4mod->o_mats[0]); */
441
442  /* Multiply each of these matrices by a relative substitution rate */
443  for(i=1;i<m4mod->n_h;i++) For(j,n_o*n_o) m4mod->o_mats[i][j] = m4mod->o_mats[0][j]*m4mod->multipl[i];
444  For(j,n_o*n_o) m4mod->o_mats[0][j] *= m4mod->multipl[0];
445
446  For(i,n_s*n_s) mod->r_mat->qmat->v[i] = .0;
447
448  /* Diagonal blocks (i.e, nucleotide substitutions), symmetric */
449  For(i,n_s)
450    {
451      for(j=i+1;j<n_s;j++)
452        {
453          if((int)(j/n_o) == (int)(i/n_o))
454            {
455              mod->r_mat->qmat->v[i*n_s+j] = m4mod->o_mats[(int)(i/n_o)][(i%n_o)*n_o+j%n_o];
456              mod->r_mat->qmat->v[j*n_s+i] = mod->r_mat->qmat->v[i*n_s+j] * m4mod->o_fq[i%n_o] / m4mod->o_fq[j%n_o];
457            }
458        }
459    }
460
461  /* Work out scaling factor such that the  expected number of observed state substitution
462     along a branch of length 1 is 1.*/
463  mr = .0;
464  For(i,n_s)
465    {
466      sum = .0;
467      For(j,n_s) sum += mod->r_mat->qmat->v[i*n_s+j];
468      mr += sum * m4mod->o_fq[i%n_o] * m4mod->h_fq[(int)(i/n_o)]; 
469    }
470 
471  /* Scale the diagonal blocks */
472  For(i,n_s*n_s) mod->r_mat->qmat->v[i] /= mr;
473 
474  /* We are done with the diagonal blocks. Let's fill the non-diagonal ones now. */
475
476  /* Fill the matrix of substitution rate across classes (switches) here */
477  Update_Qmat_Generic(m4mod->h_rr, m4mod->h_fq, m4mod->n_h, m4mod->h_mat);
478
479/*   Print_Square_Matrix_Generic(m4mod->n_h,m4mod->h_mat); */
480
481  /* Multiply this matrix by the switching rate */
482  For(i,n_h*n_h) m4mod->h_mat[i] *= m4mod->delta;
483
484  /* Fill the non diagonal blocks */
485  For(i,n_s)
486    {
487      for(j=i+1;j<n_s;j++)
488        {
489          if((int)(j/n_o) != (int)(i/n_o))
490            {
491              if(i%n_o == j%n_o)
492                {
493                  mod->r_mat->qmat->v[i*n_s+j] = m4mod->h_mat[(int)(i/n_o)*n_h+(int)(j/n_o)];
494                  mod->r_mat->qmat->v[j*n_s+i] = mod->r_mat->qmat->v[i*n_s+j] * m4mod->h_fq[(int)(i/n_o)] / m4mod->h_fq[(int)(j/n_o)]; 
495                }
496            }
497        }
498    }
499
500  /* Note: class equilibrium frequencies are already built in the h_mat matrix.
501     No need to 'add' these frequencies later on. */
502
503  /* We are done with the non diagonal blocks */
504
505
506  /* Diagonal cells */
507  For(i,n_s)
508    {
509      sum = .0;
510      For(j,n_s)
511        {
512          if(j != i)
513            sum += mod->r_mat->qmat->v[i*n_s+j];
514        }
515      mod->r_mat->qmat->v[i*n_s+i] = -sum;
516    }
517
518  For(i,n_s*n_s) mod->eigen->q[i] = mod->r_mat->qmat->v[i];
519}
520
521//////////////////////////////////////////////////////////////
522//////////////////////////////////////////////////////////////
523
524
525void M4_Init_P_Lk_Tips_Double(t_tree *tree)
526{
527  int curr_site,i,j,k,l,dim1,dim2,dim3;
528 
529  dim1 = tree->mod->ras->n_catg * tree->mod->m4mod->n_h * tree->mod->m4mod->n_o;
530  dim2 = tree->mod->m4mod->n_h * tree->mod->m4mod->n_o;
531  dim3 = tree->mod->m4mod->n_o;
532
533
534  Fors(curr_site,tree->data->crunch_len,tree->mod->io->state_len)
535    {
536      For(i,tree->n_otu)
537        {
538          for(j=1;j<tree->mod->m4mod->n_h;j++)
539            {
540              For(k,tree->mod->m4mod->n_o)
541                {
542                  tree->a_nodes[i]->b[0]->p_lk_rght[curr_site*dim1 + 0*dim2 + j*dim3+k] = 
543                    tree->a_nodes[i]->b[0]->p_lk_rght[curr_site*dim1 + 0*dim2 + 0*dim3+k];
544                 
545                  printf("\n() i=%d plk=%f",
546                         curr_site*dim1 + 0*dim2 + j*dim3+k,
547                         tree->a_nodes[i]->b[0]->p_lk_rght[curr_site*dim1 + 0*dim2 + j*dim3+k]);
548
549                  /* tree->a_nodes[i]->b[0]->p_lk_rght[curr_site][0][j*tree->mod->m4mod->n_o+k] =  */
550                  /* tree->a_nodes[i]->b[0]->p_lk_rght[curr_site][0][k]; */
551                }
552
553
554              For(k,tree->mod->m4mod->n_o)
555                for(l=1;l<tree->mod->ras->n_catg;l++)
556                  tree->a_nodes[i]->b[0]->p_lk_rght[curr_site*dim1 + l*dim2 + j*dim3+k] = 
557                  tree->a_nodes[i]->b[0]->p_lk_rght[curr_site*dim1 + 0*dim2 + j*dim3+k];
558                  /* tree->a_nodes[i]->b[0]->p_lk_rght[curr_site][l][j*tree->mod->m4mod->n_o+k] =  */
559                  /* tree->a_nodes[i]->b[0]->p_lk_rght[curr_site][0][j*tree->mod->m4mod->n_o+k]; */
560            }
561        }
562    }
563}
564
565//////////////////////////////////////////////////////////////
566//////////////////////////////////////////////////////////////
567
568
569void M4_Init_P_Lk_Tips_Int(t_tree *tree)
570{
571  int curr_site,i,j,k,dim2,dim3;
572
573  dim2 = tree->mod->m4mod->n_h * tree->mod->m4mod->n_o;
574  dim3 = tree->mod->m4mod->n_o;
575
576  Fors(curr_site,tree->data->crunch_len,tree->mod->io->state_len)
577    {
578      For(i,tree->n_otu)
579        {
580          for(j=1;j<tree->mod->m4mod->n_h;j++)
581            {
582              For(k,tree->mod->m4mod->n_o)
583                {
584                  tree->a_nodes[i]->b[0]->p_lk_tip_r[curr_site*dim2 + j*dim3+k] = 
585                    tree->a_nodes[i]->b[0]->p_lk_tip_r[curr_site*dim2 + 0*dim3+k];
586                  /* tree->a_nodes[i]->b[0]->p_lk_tip_r[curr_site][j*tree->mod->m4mod->n_o+k] =  */
587                  /* tree->a_nodes[i]->b[0]->p_lk_tip_r[curr_site][k]; */
588                }
589            }
590        }
591    }
592}
593
594//////////////////////////////////////////////////////////////
595//////////////////////////////////////////////////////////////
596
597
598phydbl ****M4_Integral_Term_On_One_Edge(t_edge *b, t_tree *tree)
599{
600  phydbl ****integral,*P1,*P2; 
601  int ns;
602  int g,i,j,k,l;
603  int step;
604 
605
606  ns = tree->mod->ns;
607
608
609  P1 = (phydbl *)mCalloc(tree->mod->ras->n_catg*ns*ns,sizeof(phydbl));
610  P2 = (phydbl *)mCalloc(tree->mod->ras->n_catg*ns*ns,sizeof(phydbl));
611
612
613  integral = (phydbl ****)mCalloc(tree->mod->ras->n_catg,sizeof(phydbl ***));
614  For(g,tree->mod->ras->n_catg)
615    {
616      integral[g] = (phydbl ***)mCalloc(ns,sizeof(phydbl **));
617      For(j,ns)
618        {
619          integral[g][j] = (phydbl **)mCalloc(ns,sizeof(phydbl *));
620          For(k,ns) integral[g][j][k] = (phydbl *)mCalloc(ns,sizeof(phydbl));
621        }
622    }
623
624  /* Integral calculation */
625  step = 100;
626
627  PhyML_Printf("\n. [");
628  For(i,step)
629    {
630      For(g,tree->mod->ras->n_catg)
631        {
632          PMat(((phydbl)(i+0.5)/step)*b->l->v*tree->mod->ras->gamma_rr->v[g],tree->mod,g*ns*ns,P1);
633          PMat(((phydbl)(step-i-0.5)/step)*b->l->v*tree->mod->ras->gamma_rr->v[g],tree->mod,g*ns*ns,P2);
634
635          For(j,ns)
636            {
637              For(k,ns)
638                {
639                  For(l,ns)
640                    {
641                      /* integral[g][j][k][l] += P1[g][j][k] * P2[g][j][l]  / ((phydbl)(step)); */
642                      integral[g][j][k][l] += P1[g*ns*ns + j*ns+k] * P2[g*ns*ns + j*ns+l] / ((phydbl)(step));
643                    }
644                }
645            }     
646        }
647      PhyML_Printf("."); fflush(NULL);
648    }
649  PhyML_Printf("]\n");
650
651  Free(P1);
652  Free(P2);
653
654  return integral;
655}
656
657//////////////////////////////////////////////////////////////
658//////////////////////////////////////////////////////////////
659
660
661void M4_Post_Prob_H_Class_Edge_Site(t_edge *b, phydbl ****integral, phydbl *postprob, t_tree *tree)
662{
663  /* Calculation of the expected frequencies of each hidden
664     class at a given site. */
665
666  phydbl site_lk;
667  int g,i,j,k,l;
668  int n_h;
669  phydbl sum;
670  int dim1,dim2;
671
672  dim1 = tree->mod->ras->n_catg * tree->mod->ns;
673  dim2 = tree->mod->ns;
674
675  n_h = tree->mod->m4mod->n_h; /* number of classes, i.e., number of hidden states */
676
677  site_lk = (phydbl)EXP(tree->cur_site_lk[tree->curr_site]);
678
679  if(b->rght->tax)
680    {
681      sum = .0;
682      For(i,n_h)
683        {
684          postprob[i] = .0;
685          For(j,tree->mod->m4mod->n_o)
686            {
687              For(g,tree->mod->ras->n_catg)
688                {
689                  For(k,tree->mod->ns)
690                    {
691                      For(l,tree->mod->ns)
692                        {
693                          postprob[i] +=
694
695                            (1./site_lk) *
696                            tree->mod->ras->gamma_r_proba->v[g] *
697                            tree->mod->m4mod->h_fq[i] *
698                            tree->mod->m4mod->o_fq[j] *
699                            b->p_lk_left[tree->curr_site*dim1 + g*dim2 + k] *
700                            b->p_lk_tip_r[tree->curr_site*dim2 + l] *
701                            integral[g][i*tree->mod->m4mod->n_o+j][k][l];
702
703                            /* (1./site_lk) * */
704                            /* tree->mod->ras->gamma_r_proba[g] * */
705                            /* tree->mod->m4mod->h_fq[i] * */
706                            /* tree->mod->m4mod->o_fq[j] * */
707                            /* b->p_lk_left[tree->curr_site][g][k] * */
708                            /* b->p_lk_tip_r[tree->curr_site][l] * */
709                            /* /\*                      b->p_lk_rght[tree->curr_site][0][l] * *\/ */
710                            /* integral[g][i*tree->mod->m4mod->n_o+j][k][l]; */
711                        }
712                    }
713                }
714            }
715          sum += postprob[i];
716        }
717
718      /* TO DO */
719      For(i,n_h) postprob[i] *= EXP(b->sum_scale_left[tree->curr_site]); 
720
721    }
722  else
723    {
724      sum = .0;
725      For(i,n_h)
726        {
727          postprob[i] = .0;
728          For(j,tree->mod->m4mod->n_o)
729            {
730              For(g,tree->mod->ras->n_catg)
731                {
732                  For(k,tree->mod->ns)
733                    {
734                      For(l,tree->mod->ns)
735                        {
736                          postprob[i] +=
737
738                            (1./site_lk) *
739                            tree->mod->ras->gamma_r_proba->v[g] *
740                            tree->mod->m4mod->h_fq[i] *
741                            tree->mod->m4mod->o_fq[j] *
742                            b->p_lk_left[tree->curr_site*dim1 + g*dim2 + k] *
743                            b->p_lk_rght[tree->curr_site*dim1 + g*dim2 + l] *
744                            integral[g][i*tree->mod->m4mod->n_o+j][k][l];
745
746                            /* (1./site_lk) * */
747                            /* tree->mod->ras->gamma_r_proba[g] * */
748                            /* tree->mod->m4mod->h_fq[i] * */
749                            /* tree->mod->m4mod->o_fq[j] * */
750                            /* b->p_lk_left[tree->curr_site][g][k] * */
751                            /* b->p_lk_rght[tree->curr_site][g][l] * */
752                            /* integral[g][i*tree->mod->m4mod->n_o+j][k][l]; */
753                        }
754                    }
755                }
756            }
757          sum += postprob[i];
758        }
759
760      /* TO DO */
761      For(i,n_h) postprob[i] *= EXP(b->sum_scale_left[tree->curr_site] + b->sum_scale_rght[tree->curr_site]); 
762
763    }
764
765  For(i,n_h) 
766    if((postprob[i] < -1.E-5) || (postprob[i] > 1.0+1.E-5))
767      {
768        PhyML_Printf("\n. Cat : %d Prob : %f\n",i,postprob[i]);
769        PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
770        Warn_And_Exit("\n");
771      }
772
773  sum = 0.0;
774  For(i,n_h) sum += postprob[i];
775
776  if((sum > 1.0+1.E-2) || (sum < 1.0-1.E-2))
777    {
778      PhyML_Printf("\n. Sum = %f\n",sum);
779      PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
780      Exit("\n");
781    }
782
783  return;
784}
785
786//////////////////////////////////////////////////////////////
787//////////////////////////////////////////////////////////////
788
789
790phydbl ***M4_Compute_Proba_Hidden_States_On_Edges(t_tree *tree)
791{
792  int i;
793  phydbl ***post_probs;
794  phydbl ****integral;
795
796
797  post_probs = (phydbl ***)mCalloc(2*tree->n_otu-3,sizeof(phydbl **));
798
799  For(i,2*tree->n_otu-3)
800    {
801      post_probs[i] = (phydbl **)mCalloc(tree->n_pattern,sizeof(phydbl *));
802      For(tree->curr_site,tree->n_pattern) 
803        post_probs[i][tree->curr_site] = (phydbl *)mCalloc(tree->mod->m4mod->n_h,sizeof(phydbl));
804    }
805
806
807  /* Compute posterior probabilities of each hidden class (usually a rate class)
808     on each edge, at each site.
809  */
810  For(i,2*tree->n_otu-3) 
811    {
812      PhyML_Printf("\n. Edge %4d/%4d",i+1,2*tree->n_otu-3);
813
814      integral = M4_Integral_Term_On_One_Edge(tree->a_edges[i],tree);
815
816      For(tree->curr_site,tree->n_pattern)
817        M4_Post_Prob_H_Class_Edge_Site(tree->a_edges[i],
818                                       integral,
819                                       post_probs[i][tree->curr_site],
820                                       tree);
821     
822      M4_Free_Integral_Term_On_One_Edge(integral,tree);
823    }
824  return post_probs;
825}
826
827//////////////////////////////////////////////////////////////
828//////////////////////////////////////////////////////////////
829
830
831/* Estimate the (posterior) mean relative rate of substitution on each branch
832   at each site. The posterior mean rates averaged over sites is also estimated
833   for each edge. The corresponding trees are printed in a postscript file. Tree 0
834   is the tree with posterior mean rates averaged over the sites. The following trees
835   have posterior mean rates computed for each site.
836*/
837void M4_Compute_Posterior_Mean_Rates(phydbl ***post_probs, t_tree *tree)
838{
839  char *s;
840  int i;
841  phydbl **mean_post_probs;
842  phydbl *mrr;
843  phydbl sum;
844  int patt,br,rcat;
845  phydbl *mean_br_len;
846  int best_r,len_var;
847  phydbl max_prob;
848
849  mean_br_len = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
850  mean_post_probs = (phydbl **)mCalloc(2*tree->n_otu-3,sizeof(phydbl *));
851  For(i,2*tree->n_otu-3) mean_post_probs[i] = (phydbl *)mCalloc(tree->mod->m4mod->n_h,sizeof(phydbl ));
852  mrr = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
853
854  Record_Br_Len(tree);
855  M4_Scale_Br_Len(tree);
856
857  /* Compute the posterior mean relative rate on each branch averaged over the
858     whole set of patterns (sites) */
859  len_var = 0;
860  For(patt,tree->n_pattern) 
861    {
862      if(!Is_Invar(patt,1,NT,tree->data))
863        {
864          For(br,2*tree->n_otu-3)
865            {
866              max_prob = -1.;
867              best_r = -1;
868              For(rcat,tree->mod->m4mod->n_h)
869                {
870                  if(post_probs[br][patt][rcat] > max_prob)
871                    {
872                      max_prob = post_probs[br][patt][rcat];
873                      best_r = rcat;
874                    }
875                }
876
877/* /\*        Add weight on each category, weight is proportional to the corresponding posterior probability *\/ */
878/*            For(rcat,tree->mod->m4mod->n_h) */
879/*              { */
880/*                mean_post_probs[br][rcat] += post_probs[br][patt][rcat] * tree->data->wght[patt]; */
881/*              } */
882
883              /* Add weight on the most probable rate category only */
884              mean_post_probs[br][best_r] += tree->data->wght[patt];
885            }
886          len_var += tree->data->wght[patt];
887        }
888    }
889
890  For(br,2*tree->n_otu-3) 
891    {
892      For(rcat,tree->mod->m4mod->n_h)
893        {
894          mean_post_probs[br][rcat] /= (phydbl)len_var;
895        }
896    }
897
898  /* Compute the posterior mean relative rate and scale
899     each branch length using this factor */
900  For(br,2*tree->n_otu-3)
901    {
902      For(rcat,tree->mod->m4mod->n_h)
903        {
904          mrr[br] += mean_post_probs[br][rcat] * tree->mod->m4mod->multipl[rcat];
905        }
906      tree->a_edges[br]->l->v *= mrr[br];
907    }
908
909  PhyML_Fprintf(tree->io->fp_out_stats,"\n. Mean posterior probabilities & rates\n");
910  For(rcat,tree->mod->m4mod->n_h) PhyML_Fprintf(tree->io->fp_out_stats,"%2.4f ",tree->mod->m4mod->multipl[rcat]);
911  PhyML_Fprintf(tree->io->fp_out_stats,"\n");
912  For(br,2*tree->n_otu-3) 
913    {
914      For(rcat,tree->mod->m4mod->n_h)
915        {
916          PhyML_Fprintf(tree->io->fp_out_stats,"%2.4f ",mean_post_probs[br][rcat]);
917        }
918/*       PhyML_Fprintf(tree->io->fp_out_stats," -- %f -> %f x %f = %f",mrr[br],tree->a_edges[br]->l->v,mrr[br],tree->a_edges[br]->l->v*mrr[br]); */
919
920      PhyML_Fprintf(tree->io->fp_out_stats," mrr=%f ",mrr[br]);
921
922      if(mrr[br] > 1.) PhyML_Fprintf(tree->io->fp_out_stats,"FAST ");
923      else             PhyML_Fprintf(tree->io->fp_out_stats,"SLOW ");
924     
925      PhyML_Fprintf(tree->io->fp_out_stats,"%s",tree->a_edges[br]->labels[0]);
926
927      PhyML_Fprintf(tree->io->fp_out_stats,"\n");
928    }
929
930  /* Print the tree */
931  PhyML_Fprintf(tree->io->fp_out_tree,"Constrained tree with corrected branch lengths = ");
932  s = Write_Tree(tree,NO);
933  PhyML_Fprintf(tree->io->fp_out_tree,"%s\n",s);
934  Free(s);
935  tree->ps_tree = DR_Make_Tdraw_Struct(tree);
936  DR_Print_Postscript_Header(tree->n_pattern,tree->io->fp_out_ps);
937  tree->ps_page_number = 0;
938  DR_Print_Tree_Postscript(tree->ps_page_number++,YES,tree->io->fp_out_ps,tree);
939
940  /* Go back to the initial scaled branch lengths */
941  For(br,2*tree->n_otu-3) tree->a_edges[br]->l->v /= mrr[br];
942
943  /* Compute the posterior mean relative rate at each site, for each branch
944     and each rate category. Scale branch lengths using these factors and
945     print each tree (i.e., on tree per site pattern) */
946  For(patt,tree->n_pattern) 
947    {
948      For(br,2*tree->n_otu-3) 
949        {
950          mrr[br] = .0;
951          max_prob = -1.;
952          best_r = -1;
953          For(rcat,tree->mod->m4mod->n_h) /* For each rate class */
954            {
955              mrr[br] += post_probs[br][patt][rcat] * tree->mod->m4mod->multipl[rcat];
956              if(post_probs[br][patt][rcat] > max_prob)
957                {
958                  max_prob = post_probs[br][patt][rcat];
959                  best_r = rcat;
960                }
961            }
962/*        mrr[br] = tree->mod->m4mod->multipl[best_r]; /\* Use the most probable relative rate instead of mean *\/ */
963          tree->a_edges[br]->l->v *= mrr[br];
964        }
965
966      For(br,2*tree->n_otu-3) mean_br_len[br] += tree->a_edges[br]->l->v * tree->data->wght[patt];
967
968      PhyML_Fprintf(tree->io->fp_out_stats,"\n. Posterior probabilities site %4d (weight=%d, is_inv=%d)\n",
969             patt,
970             tree->data->wght[patt],
971             Is_Invar(patt,1,NT,tree->data));
972
973      For(rcat,tree->mod->m4mod->n_h) PhyML_Fprintf(tree->io->fp_out_stats,"%2.4f ",tree->mod->m4mod->multipl[rcat]);
974      PhyML_Fprintf(tree->io->fp_out_stats,"\n");
975      For(br,2*tree->n_otu-3)
976        {
977          PhyML_Fprintf(tree->io->fp_out_stats,"Edge %3d ",br);
978          max_prob = -1.0;
979          best_r = -1;
980          For(rcat,tree->mod->m4mod->n_h)
981            {
982              if(post_probs[br][patt][rcat] > max_prob)
983                {
984                  max_prob = post_probs[br][patt][rcat];
985                  best_r = rcat;
986                }
987            }
988
989          For(rcat,tree->mod->m4mod->n_h)
990            {
991              PhyML_Fprintf(tree->io->fp_out_stats,"%2.4f",post_probs[br][patt][rcat]);
992              if(rcat == best_r) PhyML_Fprintf(tree->io->fp_out_stats,"* ");
993              else               PhyML_Fprintf(tree->io->fp_out_stats,"  ");
994            }
995
996/*        PhyML_Fprintf(tree->io->fp_out_stats," -- %f -> %f x %f = %f",mrr[br],tree->a_edges[br]->l->v,mrr[br],tree->a_edges[br]->l->v*mrr[br]); */
997         
998          if(mrr[br] > 1.01)      PhyML_Fprintf(tree->io->fp_out_stats," %s ","FAST");
999          else if(mrr[br] < 0.99) PhyML_Fprintf(tree->io->fp_out_stats," %s ","SLOW");
1000          else                    PhyML_Fprintf(tree->io->fp_out_stats," %s ","MEDIUM");
1001          PhyML_Fprintf(tree->io->fp_out_stats,"%s ",tree->a_edges[br]->labels[0]);
1002          PhyML_Fprintf(tree->io->fp_out_stats,"\n");
1003        }
1004
1005      PhyML_Fprintf(tree->io->fp_out_tree,"tree %d = ",patt+1);
1006      s = Write_Tree(tree,NO);
1007      PhyML_Fprintf(tree->io->fp_out_tree,"%s\n",s);
1008      Free(s);
1009      DR_Print_Tree_Postscript(tree->ps_page_number++,YES,tree->io->fp_out_ps,tree);
1010
1011      /* Go back to the initial scaled branch lengths */
1012      For(br,2*tree->n_otu-3) tree->a_edges[br]->l->v /= mrr[br];
1013
1014      For(br,2*tree->n_otu-3) 
1015        {
1016          sum = .0;
1017          For(rcat,tree->mod->m4mod->n_h)
1018            {
1019              sum += post_probs[br][patt][rcat];
1020            }
1021         
1022          if((sum < 0.99) || (sum > 1.01))
1023            {
1024              PhyML_Fprintf(tree->io->fp_out_stats,"\n. sum = %f\n",sum);
1025              PhyML_Fprintf(tree->io->fp_out_stats,"\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1026              Warn_And_Exit("\n");
1027            }
1028        }
1029    }
1030 
1031  /* Mean branch lengths */
1032  For(br,2*tree->n_otu-3)
1033    {
1034      mean_br_len[br] /= (phydbl)tree->data->init_len;
1035      tree->a_edges[br]->l->v = mean_br_len[br];
1036    }
1037  PhyML_Fprintf(tree->io->fp_out_tree,"Mean branch lengths=");
1038  s = Write_Tree(tree,NO);
1039  PhyML_Fprintf(tree->io->fp_out_tree,"%s\n",s);
1040  Free(s);
1041/*   DR_Print_Tree_Postscript(tree->ps_page_number++,tree->io->fp_out_ps,tree); */
1042
1043  Restore_Br_Len(tree);
1044
1045  DR_Print_Postscript_EOF(tree->io->fp_out_ps);
1046
1047  For(br,2*tree->n_otu-3)
1048    {
1049      For(tree->curr_site,tree->n_pattern)
1050        Free(post_probs[br][tree->curr_site]);
1051      Free(post_probs[br]);
1052    }
1053  Free(post_probs);
1054  For(i,2*tree->n_otu-3) Free(mean_post_probs[i]);
1055  Free(mean_post_probs);
1056  Free(mrr);
1057  Free(mean_br_len);
1058}
1059
1060//////////////////////////////////////////////////////////////
1061//////////////////////////////////////////////////////////////
1062
1063
1064/* Classifiy each branch, at each site, among one of the rate classes */
1065phydbl **M4_Site_Branch_Classification(phydbl ***post_probs, t_tree *tree)
1066{
1067  int patt, br, rcat, i;
1068  phydbl **best_probs;
1069  phydbl post_prob_fast, post_prob_slow;
1070
1071  best_probs = (phydbl **)mCalloc(tree->n_pattern,sizeof(phydbl *));
1072  For(i,tree->n_pattern) best_probs[i] = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
1073
1074  tree->write_labels = YES;
1075
1076  For(patt,tree->n_pattern) 
1077    {
1078      For(br,2*tree->n_otu-3) 
1079        {
1080          post_prob_fast = .0;
1081          post_prob_slow = .0;
1082
1083          For(rcat,tree->mod->m4mod->n_h) /* For each rate class */
1084            {         
1085              if(tree->mod->m4mod->multipl[rcat] > 1.0) 
1086                post_prob_fast += post_probs[br][patt][rcat];
1087              else
1088                post_prob_slow += post_probs[br][patt][rcat];
1089            }
1090
1091          best_probs[patt][br] = (post_prob_fast > post_prob_slow)?(post_prob_fast):(post_prob_slow);
1092
1093          if(!(tree->a_edges[br]->n_labels%BLOCK_LABELS)) Make_New_Edge_Label(tree->a_edges[br]);
1094
1095/*        if((post_prob_fast > post_prob_slow) && (best_probs[patt][br] > 0.95)) */
1096/*          strcpy(tree->a_edges[br]->labels[tree->a_edges[br]->n_labels],"FASTER"); */
1097/*        else if((post_prob_fast < post_prob_slow) && (best_probs[patt][br] > 0.95)) */
1098/*          strcpy(tree->a_edges[br]->labels[tree->a_edges[br]->n_labels],"SLOWER"); */
1099/*        else */
1100/*          strcpy(tree->a_edges[br]->labels[tree->a_edges[br]->n_labels],"UNKNOWN"); */
1101
1102          if(post_prob_fast > post_prob_slow)
1103            strcpy(tree->a_edges[br]->labels[tree->a_edges[br]->n_labels],"FASTER");
1104          else 
1105            strcpy(tree->a_edges[br]->labels[tree->a_edges[br]->n_labels],"SLOWER");
1106
1107          tree->a_edges[br]->n_labels++;
1108        }
1109    }
1110  return best_probs;
1111}
1112
1113//////////////////////////////////////////////////////////////
1114//////////////////////////////////////////////////////////////
1115
1116
1117void M4_Site_Branch_Classification_Experiment(t_tree *tree)
1118{
1119  calign *cpy_data;
1120  short int **true_rclass, **est_rclass;
1121  phydbl **best_probs;
1122  int i,j;
1123  phydbl correct_class, mis_class, unknown;
1124 
1125  true_rclass = (short int **)mCalloc(tree->data->init_len, sizeof(short int *));
1126  est_rclass  = (short int **)mCalloc(tree->data->init_len, sizeof(short int *));
1127 
1128  For(i,tree->data->init_len)
1129    {
1130      true_rclass[i] = (short int *)mCalloc(2*tree->n_otu-3,sizeof(short int));
1131      est_rclass[i]  = (short int *)mCalloc(2*tree->n_otu-3,sizeof(short int));
1132    }
1133
1134  if(tree->io->datatype != NT && tree->io->datatype != AA)
1135    {
1136      PhyML_Printf("\n. Not implemented yet.");
1137      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1138      Warn_And_Exit("");
1139    }
1140 
1141  cpy_data = Copy_Cseq(tree->data,tree->io);
1142
1143  /* Generate a simulated data set under H0, with the right sequence length. */
1144  PhyML_Printf("\n. Evolving sequences (delta=%f, alpha=%f) ...\n",tree->mod->m4mod->delta,tree->mod->m4mod->alpha);
1145  Evolve(cpy_data,tree->mod,tree);
1146
1147  For(i,cpy_data->init_len)
1148    {
1149      For(j,2*tree->n_otu-3)
1150        {
1151          if(!strcmp(tree->a_edges[j]->labels[i],"FASTER"))
1152            {
1153              true_rclass[i][j] = 1;
1154            }
1155          else if(!strcmp(tree->a_edges[j]->labels[i],"SLOWER"))
1156            {
1157              true_rclass[i][j] = 0;
1158            }
1159          else
1160            {
1161              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1162              Warn_And_Exit("\n");
1163            }
1164        }
1165    }
1166 
1167  For(j,2*tree->n_otu-3) 
1168    {
1169      Free_Edge_Labels(tree->a_edges[j]);
1170      tree->a_edges[j]->n_labels = 0;
1171    }
1172
1173  /* Generate the memory needed for likelihood calculation because
1174     we will need bigger arrays
1175  */
1176  Free_Tree_Lk(tree);
1177  Free_Tree_Pars(tree);
1178
1179  tree->data      = cpy_data;
1180  tree->n_pattern = tree->data->crunch_len;
1181
1182  /* Allocate memory and initialize likelihood structure with
1183     simulated sequences (i.e., columns are not compressed)
1184  */
1185  Make_Tree_4_Pars(tree,cpy_data,cpy_data->init_len);
1186  Make_Tree_4_Lk(tree,cpy_data,cpy_data->init_len);
1187
1188  /* Estimate model parameters */
1189  PhyML_Printf("\n. Estimating model parameters...\n");
1190  tree->mod->s_opt->opt_cov_alpha = 1;
1191  tree->mod->s_opt->opt_cov_delta = 1;
1192  Round_Optimize(tree,tree->data,ROUND_MAX);
1193
1194  tree->both_sides = 1;
1195  Lk(NULL,tree);
1196
1197  /* Classify branches */
1198  best_probs = M4_Site_Branch_Classification(M4_Compute_Proba_Hidden_States_On_Edges(tree),tree);
1199
1200  For(i,tree->data->init_len)
1201    {
1202      For(j,2*tree->n_otu-3)
1203        {
1204          if(!strcmp(tree->a_edges[j]->labels[i],"FASTER"))
1205            {
1206              est_rclass[i][j] = 1;
1207            }
1208          else if(!strcmp(tree->a_edges[j]->labels[i],"SLOWER"))
1209            {
1210              est_rclass[i][j] = 0;
1211            }
1212          else if(!strcmp(tree->a_edges[j]->labels[i],"UNKNOWN"))
1213            {
1214              est_rclass[i][j] = -1;
1215            }
1216          else
1217            {
1218              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1219              Warn_And_Exit("\n");
1220            }
1221        }
1222    }
1223
1224  unknown       = .0;
1225  correct_class = .0;
1226  mis_class     = .0;
1227  For(i,tree->data->init_len)
1228    {
1229      For(j,2*tree->n_otu-3)
1230        {
1231/*        PhyML_Printf("\n. Edge %3d %4d %4d - %f",j,true_rclass[i][j],est_rclass[i][j],best_probs[i][j]); */
1232          if(est_rclass[i][j] == -1)
1233            {
1234              unknown += 1.;
1235            }
1236          else if(est_rclass[i][j] == true_rclass[i][j])
1237            {
1238              correct_class += 1.;
1239            }
1240          else if(est_rclass[i][j] != true_rclass[i][j])
1241            {
1242              mis_class += 1.;
1243            }
1244          else
1245            {
1246              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
1247              Warn_And_Exit("\n");
1248            }
1249        }
1250/*       PhyML_Printf("\n"); */
1251    }
1252
1253  correct_class /= ((tree->data->init_len * (2*tree->n_otu-3)) - unknown);
1254  mis_class     /= ((tree->data->init_len * (2*tree->n_otu-3)) - unknown);
1255  unknown       /= (tree->data->init_len  * (2*tree->n_otu-3));
1256 
1257  PhyML_Printf("\n. correct_class = %f mis_class = %f unknown = %f",
1258         correct_class,
1259         mis_class,
1260         unknown);
1261
1262
1263  For(i,tree->data->init_len)
1264    {
1265      Free(true_rclass[i]);
1266      Free(est_rclass[i]);
1267      Free(best_probs[i]);
1268    }
1269  Free(true_rclass);
1270  Free(est_rclass);
1271  Free(best_probs);
1272
1273}
1274
1275//////////////////////////////////////////////////////////////
1276//////////////////////////////////////////////////////////////
1277
1278
1279  /* Scale branch lengths such that they express expected number
1280     of nucleotide or amino-acid substitutions */
1281
1282void M4_Scale_Br_Len(t_tree *tree)
1283{
1284  phydbl scale_fact,mrs;
1285  int i,j;
1286
1287  /* (1) Work out the relative mean rate of switches */
1288  mrs = .0;
1289  For(i,tree->mod->m4mod->n_h)
1290    {
1291      For(j,tree->mod->m4mod->n_h)
1292        {
1293          if(j != i)
1294            mrs += tree->mod->m4mod->h_fq[i] * tree->mod->m4mod->h_mat[i*tree->mod->m4mod->n_h+j];
1295        }
1296    }
1297
1298  /* (2) scale_fact = (1 + delta x mrs) */
1299  scale_fact = 1.0 + tree->mod->m4mod->delta * mrs;
1300
1301  /* (3) Scale branch lengths */
1302  For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v /= scale_fact;
1303}
1304
1305//////////////////////////////////////////////////////////////
1306//////////////////////////////////////////////////////////////
1307
1308
1309void M4_Free_Integral_Term_On_One_Edge(phydbl ****integral, t_tree *tree)
1310{
1311  int g,i,j;
1312
1313  For(g,tree->mod->ras->n_catg)
1314    {
1315      For(i,tree->mod->m4mod->n_h)
1316        {
1317          For(j,tree->mod->m4mod->n_h)
1318            {
1319              Free(integral[g][i][j]);
1320            }
1321          Free(integral[g][i]);
1322        }     
1323      Free(integral[g]);
1324    }
1325  Free(integral);
1326}
1327
1328//////////////////////////////////////////////////////////////
1329//////////////////////////////////////////////////////////////
1330
1331
1332void M4_Detect_Site_Switches_Experiment(t_tree *tree)
1333{
1334  t_mod *nocov_mod,*cov_mod,*ori_mod;
1335  calign *ori_data,*cpy_data;
1336  int i,n_iter;
1337  phydbl *nocov_bl,*cov_bl;
1338  phydbl *site_lnl_nocov, *site_lnl_cov;
1339
1340  nocov_bl       = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
1341  cov_bl         = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
1342  site_lnl_nocov = (phydbl *)mCalloc(tree->data->init_len,sizeof(phydbl));
1343  site_lnl_cov   = (phydbl *)mCalloc(tree->data->init_len,sizeof(phydbl));
1344
1345  ori_data = tree->data;
1346  ori_mod  = tree->mod;
1347
1348  if(tree->io->datatype != NT && tree->io->datatype != AA)
1349    {
1350      PhyML_Printf("\n== Not implemented yet.");
1351      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1352      Warn_And_Exit("");
1353    }
1354
1355  cpy_data = Copy_Cseq(tree->data,tree->io);
1356
1357  PhyML_Printf("\n. Estimate model parameters under non-switching substitution model.\n");
1358  Switch_From_M4mod_To_Mod(tree->mod);
1359  Simu_Loop(tree);
1360  nocov_mod = (t_mod *)Copy_Model(tree->mod); /* Record model parameters */
1361  For(i,2*tree->n_otu-3) nocov_bl[i] = tree->a_edges[i]->l->v; /* Record branch lengths */
1362  For(i,tree->data->crunch_len) site_lnl_nocov[i] = tree->cur_site_lk[i];
1363  Print_Lk(tree,"[LnL under non-switching substitution model]");
1364 
1365  PhyML_Printf("\n. Estimate model parameters under switching substitution model.\n");
1366  Switch_From_Mod_To_M4mod(tree->mod);
1367  Simu_Loop(tree);
1368  cov_mod = (t_mod *)Copy_Model(tree->mod); /* Record model parameters */
1369  For(i,2*tree->n_otu-3) cov_bl[i] = tree->a_edges[i]->l->v; /* Record branch lengths */
1370  For(i,tree->data->crunch_len) site_lnl_cov[i] = tree->cur_site_lk[i];
1371  Print_Lk(tree,"[LnL under switching substitution model]");
1372 
1373
1374  PhyML_Printf("\n");
1375  For(i,tree->data->crunch_len) PhyML_Printf("TRUTH %f %f\n",site_lnl_nocov[i],site_lnl_cov[i]);
1376
1377  /* Generate a simulated data set under H0, with the right sequence length. */
1378  tree->mod = nocov_mod;
1379  Evolve(cpy_data, nocov_mod, tree);
1380
1381  /* Generate the memory needed for likelihood calculation because
1382     we will need bigger arrays
1383  */
1384  tree->mod = cov_mod;
1385  Free_Tree_Lk(tree);
1386  Free_Tree_Pars(tree);
1387
1388  tree->data      = cpy_data;
1389  tree->n_pattern = tree->data->crunch_len;
1390  tree->mod       = cov_mod;
1391
1392  /* Allocate memory and initialize likelihood structure with
1393     simulated sequences (i.e., columns are not compressed)
1394  */
1395  Make_Tree_4_Pars(tree,cpy_data,cpy_data->init_len);
1396  Make_Tree_4_Lk(tree,cpy_data,cpy_data->init_len);
1397
1398 
1399  n_iter = 0;
1400  do
1401    {
1402      /* Get the transition proba right to generate sequences */
1403      tree->mod = nocov_mod;
1404      For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = nocov_bl[i];
1405      For(i,2*tree->n_otu-3) Update_PMat_At_Given_Edge(tree->a_edges[i],tree);
1406     
1407      /* Generate sequences */
1408      Evolve(cpy_data, nocov_mod, tree);
1409      tree->data = cpy_data;
1410
1411      if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(tree);
1412      else Init_P_Lk_Tips_Int(tree);
1413     
1414      tree->mod = nocov_mod;
1415      For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = nocov_bl[i];
1416      Lk(NULL,tree);
1417      For(i,tree->data->crunch_len) site_lnl_nocov[i] = tree->cur_site_lk[i];
1418      Print_Lk(tree,"[CPY LnL under non-switching substitution model]");
1419
1420      tree->mod = cov_mod;
1421      For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = cov_bl[i];
1422      Lk(NULL,tree);
1423      For(i,tree->data->crunch_len) site_lnl_cov[i] = tree->cur_site_lk[i];
1424      Print_Lk(tree,"[CPY LnL under switching substitution model]");
1425
1426      PhyML_Printf("\n");
1427      For(i,tree->data->crunch_len) PhyML_Printf("SYNTH %f %f\n",site_lnl_nocov[i],site_lnl_cov[i]);
1428    }
1429  while(++n_iter < 200);
1430
1431  Free_Tree_Lk(tree);
1432  Free_Tree_Pars(tree);
1433
1434  /* Back to the original data set to check that everything is ok */
1435  tree->data      = ori_data;
1436  tree->n_pattern = tree->data->crunch_len;
1437
1438  Make_Tree_4_Pars(tree,ori_data,ori_data->init_len);
1439  Make_Tree_4_Lk(tree,ori_data,ori_data->init_len);
1440
1441  tree->mod = nocov_mod;
1442  For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = nocov_bl[i]; 
1443  Lk(NULL,tree);
1444  Print_Lk(tree,"[FINAL LnL under non-switching substitution model]");
1445 
1446  tree->mod = cov_mod;
1447  For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = cov_bl[i]; 
1448  Lk(NULL,tree);
1449  Print_Lk(tree,"[FINAL LnL under switching substitution model]");
1450
1451  tree->mod = ori_mod;
1452
1453  Free_Model(cov_mod);
1454  Free_Model(nocov_mod);
1455  Free(site_lnl_cov);
1456  Free(site_lnl_nocov);
1457
1458  Free_Cseq(cpy_data);
1459  Free(nocov_bl);
1460  Free(cov_bl);
1461}
1462
1463//////////////////////////////////////////////////////////////
1464//////////////////////////////////////////////////////////////
1465
1466
1467void M4_Posterior_Prediction_Experiment(t_tree *tree)
1468{
1469  calign *ori_data,*cpy_data;
1470  int i,n_iter,n_simul;
1471  FILE *fp_nocov,*fp_cov,*fp_obs;
1472  char *s;
1473  t_edge *best_edge;
1474
1475  s = (char *)mCalloc(100,sizeof(char));
1476 
1477  best_edge = NULL;
1478
1479  strcpy(s,tree->io->in_align_file);
1480  fp_nocov = Openfile(strcat(s,"_nocov"),1);
1481  strcpy(s,tree->io->in_align_file);
1482  fp_cov   = Openfile(strcat(s,"_cov"),1);
1483  strcpy(s,tree->io->in_align_file);
1484  fp_obs = Openfile(strcat(s,"_obs"),1);
1485 
1486  Free(s);
1487
1488  Print_Diversity_Header(fp_nocov, tree);
1489  Print_Diversity_Header(fp_cov, tree);
1490  Print_Diversity_Header(fp_obs, tree);
1491
1492  ori_data = tree->data;
1493
1494  if(tree->io->datatype != NT && tree->io->datatype != AA)
1495   {
1496      PhyML_Printf("\n. Not implemented yet.");
1497      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1498      Warn_And_Exit("");
1499    }
1500
1501  cpy_data = Copy_Cseq(tree->data,tree->io);
1502
1503  /* Generate a simulated data set under H0, with the right sequence length. */
1504  Set_Model_Parameters(tree->mod);
1505  For(i,2*tree->n_otu-3) Update_PMat_At_Given_Edge(tree->a_edges[i],tree);
1506  Evolve(cpy_data,tree->mod,tree);
1507
1508  /* Generate the memory needed for likelihood calculation because
1509     we will need bigger arrays
1510  */
1511  Free_Tree_Lk(tree);
1512  Free_Tree_Pars(tree);
1513
1514  tree->data      = cpy_data;
1515  tree->n_pattern = tree->data->crunch_len;
1516
1517  /* Allocate memory and initialize likelihood structure with
1518     simulated sequences (i.e., columns are not compressed)
1519  */
1520  Make_Tree_4_Pars(tree,cpy_data,cpy_data->init_len);
1521  Make_Tree_4_Lk(tree,cpy_data,cpy_data->init_len);
1522
1523  /* Go back to the original data set */
1524  tree->data      = ori_data;
1525  tree->n_pattern = ori_data->crunch_len;
1526 
1527  if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(tree);
1528  else Init_P_Lk_Tips_Int(tree);
1529
1530  PhyML_Printf("\n. Estimate model parameters under non-switching substitution model.\n");
1531  Switch_From_M4mod_To_Mod(tree->mod);
1532
1533  tree->bl_from_node_stamps = 1;
1534  /* best_edge = TIMES_Find_Best_Root_Position(tree); */
1535  PhyML_Printf("\n. Put root on t_edge %3d",i);
1536  TIMES_Least_Square_Node_Times(best_edge,tree);
1537  TIMES_Adjust_Node_Times(tree);
1538  /* TIMES_Round_Optimize(tree); */
1539
1540/*   Round_Optimize(tree,tree->data); */
1541/*   Simu_Loop(tree); */
1542  Print_Lk(tree,"[LnL under non-switching substitution model]");
1543  Print_Tree(tree->io->fp_out_tree,tree);
1544 
1545  /* Print observed diversities */
1546  Init_Ui_Tips(tree);
1547  Site_Diversity(tree);
1548  Print_Diversity(fp_obs,tree);
1549
1550  n_simul = 100;
1551  n_iter = 0;
1552  do
1553    {
1554      Evolve(cpy_data,tree->mod,tree);
1555      tree->data      = cpy_data;
1556      tree->n_pattern = cpy_data->init_len;
1557
1558      if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(tree);
1559      else Init_P_Lk_Tips_Int(tree);
1560
1561      Lk(NULL,tree);
1562
1563      Init_Ui_Tips(tree);
1564      Site_Diversity(tree);
1565      Print_Diversity(fp_nocov,tree);
1566
1567      Print_Lk(tree,"[CPY under non-switching substitution model]");
1568    }while(++n_iter < n_simul);
1569
1570
1571  /* Go back to the original data set */
1572  tree->data      = ori_data;
1573  tree->n_pattern = ori_data->crunch_len;
1574 
1575  if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(tree);
1576  else Init_P_Lk_Tips_Int(tree);
1577
1578  PhyML_Printf("\n. Estimate model parameters under switching substitution model.\n");
1579  Switch_From_Mod_To_M4mod(tree->mod);
1580  /* TIME_Round_Optimize(tree); */
1581/*   Round_Optimize(tree,tree->data); */
1582  /*   Simu_Loop(tree); */
1583  Print_Lk(tree,"[LnL under switching substitution model]");
1584  Print_Tree(tree->io->fp_out_tree,tree);
1585
1586  n_iter = 0;
1587  do
1588    {
1589      Evolve(cpy_data,tree->mod,tree);
1590      tree->data      = cpy_data;
1591      tree->n_pattern = cpy_data->init_len;
1592      if(tree->mod->s_opt->greedy) Init_P_Lk_Tips_Double(tree);
1593      else Init_P_Lk_Tips_Int(tree);
1594
1595      Lk(NULL,tree);
1596
1597      Init_Ui_Tips(tree);
1598      Site_Diversity(tree);
1599      Print_Diversity(fp_cov,tree);
1600
1601      Print_Lk(tree,"[LnL under switching substitution model]");
1602    }while(++n_iter < n_simul);
1603
1604  fclose(fp_obs);
1605  fclose(fp_nocov);
1606  fclose(fp_cov);
1607}
1608
1609//////////////////////////////////////////////////////////////
1610//////////////////////////////////////////////////////////////
1611
1612m4 *M4_Copy_M4_Model(t_mod *ori_mod, m4 *ori_m4mod)
1613{
1614  int i,j,n_h, n_o;
1615  m4 *cpy_m4mod;
1616 
1617  if(ori_mod->io->datatype != NT && ori_mod->io->datatype != AA)
1618    {
1619      PhyML_Printf("\n== Not implemented yet.");
1620      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1621      Exit("\n");
1622    }
1623
1624
1625  cpy_m4mod = (m4 *)M4_Make_Light();
1626  cpy_m4mod->n_o = ori_m4mod->n_o;
1627  cpy_m4mod->n_h = ori_m4mod->n_h;
1628
1629  if(ori_mod->use_m4mod)
1630    {
1631      M4_Make_Complete(cpy_m4mod->n_h,
1632                       cpy_m4mod->n_o,
1633                       cpy_m4mod);
1634     
1635      n_h = cpy_m4mod->n_h;
1636      n_o = cpy_m4mod->n_o;
1637     
1638      cpy_m4mod->n_h = ori_m4mod->n_h;
1639      cpy_m4mod->n_o = ori_m4mod->n_o;
1640      For(i,n_h) For(j,n_o*n_o) cpy_m4mod->o_mats[i][j] = ori_m4mod->o_mats[i][j];
1641      For(i,n_h) cpy_m4mod->multipl[i] = ori_m4mod->multipl[i];
1642      For(i,n_h) cpy_m4mod->multipl_unscaled[i] = ori_m4mod->multipl_unscaled[i]; 
1643      For(i,n_o*n_o) cpy_m4mod->o_rr[i] = ori_m4mod->o_rr[i];
1644      For(i,n_h*n_h) cpy_m4mod->h_rr[i] = ori_m4mod->h_rr[i];
1645      For(i,n_h*n_h) cpy_m4mod->h_mat[i] = ori_m4mod->h_mat[i];
1646      For(i,n_o) cpy_m4mod->o_fq[i] = ori_m4mod->o_fq[i];
1647      For(i,n_h) cpy_m4mod->h_fq[i] = ori_m4mod->h_fq[i];
1648      For(i,n_h) cpy_m4mod->h_fq_unscaled[i] = ori_m4mod->h_fq_unscaled[i];
1649      cpy_m4mod->delta = ori_m4mod->delta;
1650      cpy_m4mod->alpha = ori_m4mod->alpha;
1651    }
1652
1653  return cpy_m4mod;
1654}
1655
1656//////////////////////////////////////////////////////////////
1657//////////////////////////////////////////////////////////////
1658
Note: See TracBrowser for help on using the repository browser.