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

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

added most recent version of phyml
File size: 136.6 KB
Line 
1
2/*
3
4PhyML:  a program that  computes maximum likelihood phylogenies from
5DNA or AA homoLOGous sequences.
6
7Copyright (C) Stephane Guindon. Oct 2003 onward.
8
9All parts of the source except where indicated are distributed under
10the GNU public licence. See http://www.opensource.org for details.
11
12*/
13
14
15
16#include "mcmc.h"
17
18//////////////////////////////////////////////////////////////
19//////////////////////////////////////////////////////////////
20
21
22void MCMC(t_tree *tree)
23{
24  int move;
25  phydbl u;
26  int first,secod;
27  int i;
28
29  RATES_Set_Clock_And_Nu_Max(tree);
30  RATES_Set_Birth_Rate_Boundaries(tree);
31 
32  if(tree->mcmc->randomize == YES)
33    {
34      MCMC_Randomize_Birth(tree);
35      MCMC_Randomize_Nu(tree);
36      MCMC_Randomize_Node_Times(tree); 
37      MCMC_Sim_Rate(tree->n_root,tree->n_root->v[2],tree);
38      MCMC_Sim_Rate(tree->n_root,tree->n_root->v[1],tree);
39      MCMC_Randomize_Node_Rates(tree);
40      MCMC_Randomize_Clock_Rate(tree); /* Clock Rate must be the last parameter to be randomized */
41      MCMC_Randomize_Rate_Across_Sites(tree);
42      MCMC_Randomize_Kappa(tree);
43      MCMC_Randomize_Covarion_Rates(tree);
44      MCMC_Randomize_Covarion_Switch(tree);
45    }
46  else
47    {
48      MCMC_Read_Param_Vals(tree);
49    }
50
51  Switch_Eigen(YES,tree->mod);
52
53  MCMC_Initialize_Param_Val(tree->mcmc,tree);
54 
55  Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
56  Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
57
58  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES; 
59  RATES_Update_Cur_Bl(tree);
60  RATES_Lk_Rates(tree); 
61  TIMES_Lk_Times(tree); 
62  Set_Both_Sides(NO,tree); 
63  if(tree->mcmc->use_data) Lk(NULL,tree);
64  else tree->c_lnL = 0.0;
65  Switch_Eigen(NO,tree->mod);
66  MCMC_Print_Param(tree->mcmc,tree);
67 
68  //////////////////
69  if(tree->io->mutmap == YES)
70    {
71      int j;
72      char *s,*t;
73      FILE *fp;
74     
75      Make_Ancestral_Seq(tree);
76      Make_MutMap(tree);
77 
78      For(j,tree->n_otu)
79        {
80          s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
81          strcpy(s,tree->a_nodes[j]->name);
82          tree->a_nodes[j]->name = s;
83        }
84     
85      s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
86      t = (char *)mCalloc(T_MAX_NAME,sizeof(char));
87     
88      tree->write_tax_names = YES;
89      For(i,tree->n_pattern)
90        {
91          For(j,tree->n_otu)
92            {
93              strcpy(t,tree->a_nodes[j]->name);
94              s[0]=tree->data->c_seq[j]->state[i];
95              s[1]='\0';
96              strcat(s,"--");
97              sprintf(s+strlen(s),"%.0f",tree->rates->nd_t[j]);
98              /* strcat(s,tree->a_nodes[j]->name); */
99              strcpy(tree->a_nodes[j]->name,s);
100            }
101         
102          strcpy(s,"rosettatree.");
103          sprintf(s+strlen(s),"%d",i);
104          fp = fopen(s,"w");
105          s = Write_Tree(tree,NO);
106          PhyML_Fprintf(fp,"%s",s);
107          fclose(fp);
108         
109          For(j,tree->n_otu) strcpy(tree->a_nodes[j]->name,t);
110        }
111      Free(s);
112      Free(t);
113    }
114
115
116  first = 2;
117  secod = 1;
118  do
119    {
120      /* if(tree->mcmc->ess[tree->mcmc->num_move_tree_height] > 100 &&  */
121      /*         tree->mcmc->ess[tree->mcmc->num_move_nu] > 100          && */
122      /*         tree->mcmc->ess[tree->mcmc->num_move_clock_r] > 100     && */
123      /*         tree->mcmc->run > 1000) */
124      /*        { */
125      /*          FILE *fp; */
126      /*          char *s; */
127
128      /*          s = (char *)mCalloc(100,sizeof(char)); */
129
130      /*          sprintf(s,"simul_par.%d",getpid()); */
131      /*          fclose(tree->mcmc->out_fp_stats); */
132      /*          tree->mcmc->out_fp_stats = fopen(s,"w"); */
133      /*          tree->mcmc->run = 0; */
134      /*          tree->mcmc->nd_t_digits = 4; */
135      /*          MCMC_Print_Param(tree->mcmc,tree); */
136
137      /*          RATES_Update_Cur_Bl(tree); */
138      /*          printf("\n. %s",Write_Tree(tree,NO)); */
139      /*          Evolve(tree->data,tree->mod,tree); */
140
141      /*          sprintf(s,"simul_seq.%d",getpid()); */
142      /*          fp = fopen(s,"w"); */
143      /*          Print_CSeq(fp,NO,tree->data); */
144      /*          fflush(NULL); */
145      /*          fclose(fp); */
146      /*          Free(s); */
147
148      /*          Exit("\n"); */
149      /*        } */
150
151
152      /* if(tree->mcmc->ess[tree->mcmc->num_move_tree_height] > 100 && */
153      /*         tree->mcmc->ess[tree->mcmc->num_move_nu] > 100          && */
154      /*         tree->mcmc->ess[tree->mcmc->num_move_clock_r] > 100     && */
155      /*         tree->mcmc->run > 1000) */
156      /*        { */
157      /*          FILE *fp; */
158      /*          char *s,*t; */
159
160      /*          s = (char *)mCalloc(100,sizeof(char)); */
161         
162      /*          t = strrchr(tree->io->in_align_file,'.'); */
163      /*          sprintf(s,"res%s",t); */
164      /*          fp = fopen(s,"w"); */
165      /*          fclose(tree->mcmc->out_fp_stats); */
166      /*          tree->mcmc->out_fp_stats = fopen(s,"w"); */
167      /*          tree->mcmc->run = 0; */
168      /*          MCMC_Print_Param(tree->mcmc,tree); */
169      /*          fclose(fp); */
170      /*          Free(s); */
171      /*          Exit("\n"); */
172      /*        } */
173
174
175      u = Uni();
176
177      For(move,tree->mcmc->n_moves) if(tree->mcmc->move_weight[move] > u) break;
178     
179      if(u < .5) { first = 2; secod = 1; }
180      else       { first = 1; secod = 2; }
181 
182
183      /* printf("\n. [%15s] %15f ",tree->mcmc->move_name[move],tree->c_lnL); */
184
185
186      /* Clock rate */
187      if(!strcmp(tree->mcmc->move_name[move],"clock"))
188        {
189          For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;       
190          MCMC_Clock_R(tree); 
191        }
192
193      /* Nu */
194      else if(!strcmp(tree->mcmc->move_name[move],"nu"))
195        {
196          For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
197          MCMC_Nu(tree);
198        }
199
200      /* Tree height */
201      else if(!strcmp(tree->mcmc->move_name[move],"tree_height"))
202        { 
203          MCMC_Tree_Height(tree);
204        }
205
206      /* Subtree height */
207      else if(!strcmp(tree->mcmc->move_name[move],"subtree_height"))
208        { 
209          MCMC_Subtree_Height(tree); 
210        }
211
212      /* Subtree rates */
213      else if(!strcmp(tree->mcmc->move_name[move],"subtree_rates"))
214        {
215          MCMC_Subtree_Rates(tree);
216        }
217
218      /* Birth rate */
219      else if(!strcmp(tree->mcmc->move_name[move],"birth_rate"))
220        {
221          MCMC_Birth_Rate(tree);
222        }
223
224      /* Swing rates */
225      else if(!strcmp(tree->mcmc->move_name[move],"tree_rates"))
226        {
227          MCMC_Tree_Rates(tree);
228        }
229
230      else if(!strcmp(tree->mcmc->move_name[move],"updown_nu_cr"))
231        {
232          MCMC_Updown_Nu_Cr(tree);
233        }
234
235      else if(!strcmp(tree->mcmc->move_name[move],"updown_t_cr"))
236        {
237          MCMC_Updown_T_Cr(tree);
238        }
239
240      else if(!strcmp(tree->mcmc->move_name[move],"updown_t_br"))
241        {
242          MCMC_Updown_T_Br(tree);
243        }
244
245      /* Ts/tv ratio */
246      else if(!strcmp(tree->mcmc->move_name[move],"kappa"))
247        {
248          MCMC_Kappa(tree);
249        }
250
251      /* Gamma shape parameter */
252      else if(!strcmp(tree->mcmc->move_name[move],"ras"))
253        {
254          MCMC_Rate_Across_Sites(tree);
255        }
256
257      /* Covarion change calibration interval */
258      else if(!strcmp(tree->mcmc->move_name[move],"jump_calibration"))
259        {
260          MCMC_Jump_Calibration(tree);
261        }
262
263      /* Covarion model parameters */
264      else if(!strcmp(tree->mcmc->move_name[move],"cov_rates"))
265        {
266          MCMC_Covarion_Rates(tree);
267        }
268
269      /* Covarion model parameters */
270      else if(!strcmp(tree->mcmc->move_name[move],"cov_switch"))
271        {
272          MCMC_Covarion_Switch(tree);
273        }
274
275
276      /* Times */
277      else if(!strcmp(tree->mcmc->move_name[move],"time"))
278        {
279          Set_Both_Sides(YES,tree);     
280          if(tree->mcmc->use_data == YES) Lk(NULL,tree);
281          Set_Both_Sides(NO,tree);     
282
283          if(tree->mcmc->is == NO || tree->rates->model_log_rates == YES)
284            {
285              MCMC_Root_Time(tree);
286              MCMC_One_Time(tree->n_root,tree->n_root->v[first],YES,tree);
287              MCMC_One_Time(tree->n_root,tree->n_root->v[secod],YES,tree);
288            }
289          else
290            {
291              /* MCMC_One_Time(tree->n_root,tree->n_root->v[first],YES,tree); */
292              /* MCMC_One_Time(tree->n_root,tree->n_root->v[secod],YES,tree); */
293              RATES_Posterior_One_Time(tree->n_root,tree->n_root->v[first],YES,tree);
294              RATES_Posterior_One_Time(tree->n_root,tree->n_root->v[secod],YES,tree);
295            }
296        }
297     
298      /* Node Rates */
299
300      else if(!strcmp(tree->mcmc->move_name[move],"nd_rate"))
301        {
302          MCMC_One_Node_Rate(tree->n_root,tree->n_root->v[first],YES,tree);
303          MCMC_One_Node_Rate(tree->n_root,tree->n_root->v[secod],YES,tree);
304        }
305
306      /* Edge Rates */
307      else if(!strcmp(tree->mcmc->move_name[move],"br_rate"))
308        {
309
310          Set_Both_Sides(YES,tree);
311          if(tree->mcmc->use_data == YES) Lk(NULL,tree);
312          Set_Both_Sides(NO,tree);
313         
314          if(tree->mcmc->is == NO)
315            {
316              /* MCMC_Slice_One_Rate(tree->n_root,tree->n_root->v[first],YES,tree); */
317              /* MCMC_Slice_One_Rate(tree->n_root,tree->n_root->v[secod],YES,tree); */
318
319              MCMC_One_Rate(tree->n_root,tree->n_root->v[first],YES,tree);
320              MCMC_One_Rate(tree->n_root,tree->n_root->v[secod],YES,tree);
321            }
322          else
323            {
324              RATES_Posterior_One_Rate(tree->n_root,tree->n_root->v[first],YES,tree);
325              RATES_Posterior_One_Rate(tree->n_root,tree->n_root->v[secod],YES,tree);
326            }
327
328          /* MCMC_Sim_Rate(tree->n_root,tree->n_root->v[2],tree); */
329          /* MCMC_Sim_Rate(tree->n_root,tree->n_root->v[1],tree); */
330          /* if(tree->mcmc->use_data == YES) Lk(NULL,tree); */
331          /* RATES_Lk_Rates(tree); */
332        }
333
334      tree->mcmc->run++;
335      MCMC_Get_Acc_Rates(tree->mcmc);
336
337      MCMC_Print_Param(tree->mcmc,tree);
338      MCMC_Print_Param_Stdin(tree->mcmc,tree);
339
340      if(tree->io->mutmap == YES)
341        {
342          if(!(tree->mcmc->run%tree->mcmc->sample_interval)) 
343            {
344              Sample_Ancestral_Seq(YES,!tree->mcmc->use_data,tree);
345             
346              phydbl sum = 0.0;
347              int edge,site,mut;
348              char *s;
349              FILE *fp;
350             
351              s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
352             
353              strcpy(s,tree->mcmc->io->in_align_file);
354              strcat(s,"_");
355              strcat(s,tree->mcmc->out_filename);
356              strcat(s,".mutmap");
357              fp = fopen(s,"w");
358             
359              Free(s);
360             
361              For(i,(2*tree->n_otu-3)*(tree->n_pattern)*6) sum += tree->mutmap[i];
362              PhyML_Fprintf(fp,"edge\t site\t mut\t count");
363              For(i,(2*tree->n_otu-3)*(tree->n_pattern)*6) 
364                {
365                  Get_Mutmap_Coord(i,&edge,&site,&mut,tree);
366                  PhyML_Fprintf(fp,"\n%4d\t %4d\t %4d\t %10f",edge,site,mut,(phydbl)tree->mutmap[i]/sum);
367                }
368             
369              fclose(fp);             
370            }
371        }
372
373      (void)signal(SIGINT,MCMC_Terminate);
374    }
375  while(tree->mcmc->run < tree->mcmc->chain_len);
376
377
378}
379
380//////////////////////////////////////////////////////////////
381//////////////////////////////////////////////////////////////
382
383void MCMC_Single_Param_Generic(phydbl *val, 
384                               phydbl lim_inf, 
385                               phydbl lim_sup, 
386                               int move_num,
387                               phydbl *lnPrior,
388                               phydbl *lnLike,
389                               phydbl (*prior_func)(t_edge *,t_tree *,supert_tree *), 
390                               phydbl (*like_func)(t_edge *,t_tree *,supert_tree *),
391                               int move_type,
392                               int _log, /* _log == YES: the model describes the distribution of log(val) but the move applies to val. Need a correction factor */
393                               t_edge *branch, t_tree *tree, supert_tree *stree)
394{
395  phydbl cur_val,new_val,new_lnLike,new_lnPrior,cur_lnLike,cur_lnPrior;
396  phydbl u,alpha,ratio;
397  phydbl K;
398  phydbl new_lnval, cur_lnval;
399 
400  Record_Br_Len(tree); 
401
402  cur_val       = *val;
403  new_val       = -1.0;
404  ratio         =  0.0;
405  K             = tree->mcmc->tune_move[move_num];
406  cur_lnval     = LOG(*val);
407  new_lnval     = cur_lnval;
408
409  if(lnLike)
410    {
411      cur_lnLike = *lnLike;
412      new_lnLike = *lnLike;
413    }
414  else
415    {
416      cur_lnLike = 0.0;
417      new_lnLike = 0.0;
418    }
419 
420  if(lnPrior)
421    {
422      cur_lnPrior = *lnPrior;
423      new_lnPrior = *lnPrior;
424    }
425  else
426    {
427      cur_lnPrior = 0.0;
428      new_lnPrior = 0.0;
429    }
430
431  MCMC_Make_Move(&cur_val,&new_val,lim_inf,lim_sup,&ratio,K,move_type);
432
433  if(new_val < lim_sup && new_val > lim_inf)
434    {
435      *val = new_val;
436 
437      RATES_Update_Cur_Bl(tree);
438           
439      if(_log == YES) ratio += (cur_lnval - new_lnval);
440     
441      if(prior_func) /* Prior ratio */
442        { 
443          new_lnPrior = (*prior_func)(branch,tree,stree); 
444          ratio += (new_lnPrior - cur_lnPrior); 
445        }
446     
447      if(like_func && tree->mcmc->use_data == YES)  /* Likelihood ratio */
448        { 
449          new_lnLike  = (*like_func)(branch,tree,stree); 
450          ratio += (new_lnLike - cur_lnLike); 
451        }
452     
453      ratio = EXP(ratio);
454      alpha = MIN(1.,ratio);
455     
456      /* printf("\n. %s cur_val: %f new_val:%f cur_lnL: %f new_lnL: %f ratio: %f", */
457      /*        tree->mcmc->move_name[move_num], */
458      /*        cur_val, */
459      /*        new_val, */
460      /*        cur_lnLike, */
461      /*        new_lnLike, */
462      /*        ratio); */
463
464      u = Uni();
465      if(u > alpha) /* Reject */
466        {
467          *val    = cur_val;
468          new_val = cur_val;
469          if(lnPrior) *lnPrior = cur_lnPrior;
470          if(lnLike)  *lnLike  = cur_lnLike;
471          Restore_Br_Len(tree);
472          if(tree->mod && tree->mod->update_eigen) Update_Eigen(tree->mod);
473        }
474      else /* Accept */
475        {
476          tree->mcmc->acc_move[move_num]++;
477          if(lnPrior) *lnPrior = new_lnPrior;
478          if(lnLike)  *lnLike  = new_lnLike;
479        }
480    }
481     
482  tree->mcmc->run_move[move_num]++;
483}
484
485//////////////////////////////////////////////////////////////
486//////////////////////////////////////////////////////////////
487
488
489void MCMC_Update_Effective_Sample_Size(int move_num, t_mcmc *mcmc, t_tree *tree)
490{
491  phydbl new_val,cur_val;
492
493  mcmc->ess_run[move_num]++;
494
495  new_val = mcmc->new_param_val[move_num];
496  cur_val = mcmc->old_param_val[move_num];
497
498  if(mcmc->ess_run[move_num] == 1)
499    {
500      mcmc->first_val[move_num] = cur_val;
501      mcmc->sum_val[move_num]   = cur_val;
502      mcmc->sum_valsq[move_num] = POW(cur_val,2);
503      return;
504    }
505
506  mcmc->sum_val[move_num]            += new_val;
507  mcmc->sum_valsq[move_num]          += POW(new_val,2);
508  mcmc->sum_curval_nextval[move_num] += cur_val * new_val;
509
510
511  mcmc->ess[move_num] = 
512    Effective_Sample_Size(mcmc->first_val[move_num],
513                          new_val,
514                          mcmc->sum_val[move_num],
515                          mcmc->sum_valsq[move_num],
516                          mcmc->sum_curval_nextval[move_num],
517                          mcmc->ess_run[move_num]);
518       
519  mcmc->old_param_val[move_num] = new_val;
520
521  if(move_num == mcmc->num_move_nd_t+tree->n_root->num-tree->n_otu)
522    {
523      /* FILE *fp; */
524      /* fp = fopen("out","a"); */
525      /* fprintf(fp,"%f\n",new_val); */
526      /* fclose(fp); */
527      /* printf("\n. first=%G last=%G sum=%f sum_valsq=%f sum_cur_next=%f run=%d ess=%f", */
528      /*             mcmc->first_val[move_num],new_val, */
529      /*             mcmc->sum_val[move_num], */
530      /*             mcmc->sum_valsq[move_num], */
531      /*             mcmc->sum_curval_nextval[move_num], */
532      /*             mcmc->run_move[move_num]+1, */
533      /*             mcmc->ess[move_num] */
534      /*             ); */
535    }
536     
537}
538
539//////////////////////////////////////////////////////////////
540//////////////////////////////////////////////////////////////
541
542void MCMC_Clock_R(t_tree *mixt_tree)
543{
544  t_tree *tree;
545
546  tree = mixt_tree;
547  do
548    {
549      MCMC_Single_Param_Generic(&(tree->rates->clock_r),
550                                mixt_tree->rates->min_clock,
551                                mixt_tree->rates->max_clock,
552                                mixt_tree->mcmc->num_move_clock_r,
553                                NULL,&(mixt_tree->c_lnL),
554                                NULL,Wrap_Lk,
555                                mixt_tree->mcmc->move_type[mixt_tree->mcmc->num_move_clock_r],
556                                NO,NULL,mixt_tree,NULL);
557
558      tree = tree->next;
559    }
560  while(tree);
561}
562
563//////////////////////////////////////////////////////////////
564//////////////////////////////////////////////////////////////
565
566#if defined(GEO)
567
568// Sample dispersal parameter from a phylogeo model
569void MCMC_Geo_Sigma(t_tree *mixt_tree)
570{
571  t_tree *tree;
572
573  tree = mixt_tree;
574  do
575    {
576      MCMC_Single_Param_Generic(&(tree->geo->sigma),
577                                mixt_tree->geo->min_sigma,
578                                mixt_tree->geo->max_sigma,
579                                mixt_tree->mcmc->num_move_geo_sigma,
580                                NULL,&(mixt_tree->geo->c_lnL),
581                                NULL,GEO_Wrap_Lk,
582                                mixt_tree->mcmc->move_type[mixt_tree->mcmc->num_move_geo_sigma],
583                                NO,NULL,mixt_tree,NULL);
584
585      GEO_Update_Fmat(tree->geo);
586      tree = tree->next;
587    }
588  while(tree);
589}
590
591//////////////////////////////////////////////////////////////
592//////////////////////////////////////////////////////////////
593
594// Sample competition parameter from a phylogeo model
595void MCMC_Geo_Lbda(t_tree *mixt_tree)
596{
597  t_tree *tree;
598
599  tree = mixt_tree;
600  do
601    {
602      MCMC_Single_Param_Generic(&(tree->geo->lbda),
603                                mixt_tree->geo->min_lbda,
604                                mixt_tree->geo->max_lbda,
605                                mixt_tree->mcmc->num_move_geo_lambda,
606                                NULL,&(mixt_tree->geo->c_lnL),
607                                NULL,GEO_Wrap_Lk,
608                                mixt_tree->mcmc->move_type[mixt_tree->mcmc->num_move_geo_lambda],
609                                NO,NULL,mixt_tree,NULL);
610
611      tree = tree->next;
612    }
613  while(tree);
614}
615
616//////////////////////////////////////////////////////////////
617//////////////////////////////////////////////////////////////
618
619// Sample global migration rate parameter from a phylogeo model
620void MCMC_Geo_Tau(t_tree *mixt_tree)
621{
622  t_tree *tree;
623
624  tree = mixt_tree;
625  do
626    {
627      MCMC_Single_Param_Generic(&(tree->geo->tau),
628                                mixt_tree->geo->min_tau,
629                                mixt_tree->geo->max_tau,
630                                mixt_tree->mcmc->num_move_geo_tau,
631                                NULL,&(mixt_tree->geo->c_lnL),
632                                NULL,GEO_Wrap_Lk,
633                                mixt_tree->mcmc->move_type[mixt_tree->mcmc->num_move_geo_tau],
634                                NO,NULL,mixt_tree,NULL);
635
636      tree = tree->next;
637    }
638  while(tree);
639}
640
641//////////////////////////////////////////////////////////////
642//////////////////////////////////////////////////////////////
643
644void MCMC_Geo_Loc(t_tree *tree)
645{
646  int target;
647  int *rec_loc; // recorded locations
648  int i;
649  phydbl cur_lnL, new_lnL;
650  phydbl u, ratio, alpha;
651  phydbl sum;
652  phydbl *probs;
653
654  cur_lnL = tree->geo->c_lnL;
655
656  rec_loc = (int *)mCalloc(2*tree->n_otu-1,sizeof(int));
657
658  For(i,2*tree->n_otu-1) rec_loc[i] = tree->geo->loc[i];
659
660  // Choose an internal node (including the root) at random
661  target = Rand_Int(tree->n_otu,2*tree->n_otu-2); 
662
663  target = 2*tree->n_otu-2;
664 
665  // Root node is special. Select new location uniformly at random
666  if(tree->a_nodes[target] == tree->n_root) 
667    {
668      probs = (phydbl *)mCalloc(tree->geo->ldscape_sz,sizeof(phydbl));
669
670      sum = 0.0;
671      For(i,tree->geo->ldscape_sz) sum += tree->geo->loc_beneath[tree->n_root->num * tree->geo->ldscape_sz + i];
672      For(i,tree->geo->ldscape_sz) probs[i] = tree->geo->loc_beneath[tree->n_root->num * tree->geo->ldscape_sz + i]/sum;
673     
674      tree->geo->loc[tree->n_root->num] = Sample_i_With_Proba_pi(probs,tree->geo->ldscape_sz);     
675
676      Free(probs);
677    }
678
679
680  // Randomize the locations below the selected node
681  GEO_Randomize_Locations(tree->a_nodes[target], 
682                          tree->geo,
683                          tree);
684 
685  new_lnL = GEO_Lk(tree->geo,tree);
686
687  ratio = (new_lnL - cur_lnL);       
688  ratio = EXP(ratio);
689  alpha = MIN(1.,ratio);     
690  u = Uni();
691
692  if(u > alpha) /* Reject */
693    {
694      For(i,2*tree->n_otu-1) tree->geo->loc[i] = rec_loc[i];
695      tree->geo->c_lnL = GEO_Lk(tree->geo,tree); // TO DO: you only need to update the occupation vector here...
696    }
697
698  Free(rec_loc);
699
700}
701
702#endif
703
704//////////////////////////////////////////////////////////////
705//////////////////////////////////////////////////////////////
706
707void MCMC_Sample_Joint_Rates_Prior(t_tree *tree)
708{
709  int i,dim;
710  phydbl T;
711  phydbl *r,*t,*lambda;
712  phydbl *min_r,*max_r;
713  phydbl k;
714
715  dim    = 2*tree->n_otu-2;
716  lambda = tree->rates->_2n_vect1;
717  min_r  = tree->rates->_2n_vect2;
718  max_r  = tree->rates->_2n_vect3;
719  r      = tree->rates->br_r;
720  t      = tree->rates->nd_t;
721
722  For(i,dim) tree->rates->mean_r[i] = 1.0;
723
724  RATES_Fill_Lca_Table(tree);
725  RATES_Covariance_Mu(tree);
726
727  T = .0;
728  For(i,dim) T += (t[tree->a_nodes[i]->num] - t[tree->a_nodes[i]->anc->num]);
729  For(i,dim) lambda[i] = (t[tree->a_nodes[i]->num] - t[tree->a_nodes[i]->anc->num])/T;
730  For(i,dim) r[i] = 1.0;
731  For(i,dim) min_r[i] = tree->rates->min_rate;
732  For(i,dim) max_r[i] = tree->rates->max_rate;
733
734  k = 1.; /* We want \sum_i lambda[i] r[i] = 1 */
735
736  Rnorm_Multid_Trunc_Constraint(tree->rates->mean_r, 
737                                tree->rates->cov_r, 
738                                min_r,max_r, 
739                                lambda,
740                                k, 
741                                r,
742                                dim);
743}
744 
745//////////////////////////////////////////////////////////////
746//////////////////////////////////////////////////////////////
747
748void MCMC_One_Rate(t_node *a, t_node *d, int traversal, t_tree *tree)
749{
750  t_edge *b;
751  int i;
752  phydbl u;
753  phydbl new_lnL_data, cur_lnL_data, new_lnL_rate, cur_lnL_rate;
754  phydbl ratio, alpha;
755  phydbl new_mu, cur_mu;
756  phydbl r_min, r_max;
757  t_edge *b1,*b2,*b3;
758  t_node *v2,*v3;
759  int move_num;
760  phydbl K;
761
762  if(tree->rates->model == STRICTCLOCK) return;
763
764  b = NULL;
765  if(a == tree->n_root) b = tree->e_root;
766  else
767    For(i,3) if(d->v[i] == a) { b = d->b[i]; break; }
768   
769  cur_mu       = tree->rates->br_r[d->num];
770  cur_lnL_data = tree->c_lnL;
771  new_lnL_data = tree->c_lnL;
772  cur_lnL_rate = tree->rates->c_lnL_rates;
773  new_lnL_rate = tree->rates->c_lnL_rates;
774  r_min        = tree->rates->min_rate;
775  r_max        = tree->rates->max_rate;
776  ratio        = 0.0;
777  move_num     = d->num+tree->mcmc->num_move_br_r;
778  K            = tree->mcmc->tune_move[move_num];
779
780  Record_Br_Len(tree);
781 
782  u = Uni();
783 
784  MCMC_Make_Move(&cur_mu,&new_mu,r_min,r_max,&ratio,K,tree->mcmc->move_type[tree->mcmc->num_move_br_r+d->num]);
785
786  /* phydbl dt,sd,mean; */
787  /* int err; */
788  /* dt     = tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num]; */
789  /* sd     = SQRT(tree->rates->nu * dt); */
790  /* mean   = tree->rates->br_r[a->num]; */
791  /* new_mu = Rnorm_Trunc(mean,sd,r_min,r_max,&err); */
792  /* ratio  = Log_Dnorm_Trunc(cur_mu,mean,sd,r_min,r_max,&err) - Log_Dnorm_Trunc(new_mu,mean,sd,r_min,r_max,&err); */
793
794  if(new_mu > r_min && new_mu < r_max)
795    {     
796      tree->rates->br_r[d->num] = new_mu;
797           
798      v2 = v3 = NULL;
799      For(i,3)
800        if((d->v[i] != a) && (d->b[i] != tree->e_root))
801          {
802            if(!v2) { v2 = d->v[i]; }
803            else    { v3 = d->v[i]; }
804          }
805     
806     
807      b1 = NULL;
808      if(a == tree->n_root) b1 = tree->e_root;
809      else For(i,3) if(d->v[i] == a) { b1 = d->b[i]; break; }
810     
811      b2 = b3 = NULL;
812      if(!d->tax)
813        {
814          For(i,3)
815            if((d->v[i] != a) && (d->b[i] != tree->e_root))
816              {
817                if(!b2) { b2 = d->b[i]; }
818                else    { b3 = d->b[i]; }
819              }
820        }
821     
822      tree->rates->br_do_updt[d->num] = YES;
823      if(!d->tax)
824        {
825          tree->rates->br_do_updt[v2->num] = YES;
826          tree->rates->br_do_updt[v3->num] = YES;
827        }
828     
829      RATES_Update_Cur_Bl(tree);
830     
831      /* printf("\n. r0=%f r1=%f cr=%f mean=%f var=%f nu=%f dt=%f", */
832      /*         r0,r1,tree->rates->clock_r,b1->gamma_prior_mean,b1->gamma_prior_var,nu,t1-t0); */
833           
834      if(tree->mcmc->use_data)
835        {
836          if(tree->io->lk_approx == EXACT)
837            {
838              Update_PMat_At_Given_Edge(b1,tree);
839              if(!d->tax)
840                {
841                  Update_PMat_At_Given_Edge(b2,tree);
842                  Update_PMat_At_Given_Edge(b3,tree);
843                }
844              Update_P_Lk(tree,b1,d);
845            }
846          new_lnL_data = Lk(b1,tree);
847         
848          /* tree->both_sides = NO; */
849          /* new_lnL_data = Lk(tree); */
850        }
851     
852      new_lnL_rate = RATES_Lk_Rates(tree);
853     
854      /* Likelihood ratio */
855      if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
856
857      /* Prior ratio */
858      ratio += (new_lnL_rate - cur_lnL_rate);
859           
860      ratio = EXP(ratio);
861      alpha = MIN(1.,ratio);
862     
863      u = Uni();
864     
865      if(u > alpha) /* Reject */
866        {
867          tree->rates->br_r[d->num] = cur_mu;     
868          tree->c_lnL               = cur_lnL_data;
869          tree->rates->c_lnL_rates  = cur_lnL_rate;
870         
871          Restore_Br_Len(tree);
872         
873          if(tree->mcmc->use_data && tree->io->lk_approx == EXACT)
874            {
875              Update_PMat_At_Given_Edge(b1,tree);
876              if(!d->tax)
877                {
878                  Update_PMat_At_Given_Edge(b2,tree);
879                  Update_PMat_At_Given_Edge(b3,tree);
880                }
881              Update_P_Lk(tree,b1,d);
882            }
883         
884          /* tree->both_sides = YES; */
885          /* new_lnL_data = Lk(tree); */
886          /* tree->both_sides = NO; */
887         
888        }
889      else
890        {
891          tree->mcmc->acc_move[tree->mcmc->num_move_br_r+d->num]++;
892        }
893    }
894  tree->mcmc->run_move[tree->mcmc->num_move_br_r+d->num]++;
895 
896  if(traversal == YES)
897    {
898      if(d->tax == YES) return;
899      else
900        {
901          For(i,3)
902            if(d->v[i] != a && d->b[i] != tree->e_root)
903              {
904                if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,d->b[i],d);
905                /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) {tree->both_sides = YES; Lk(tree); } */
906                MCMC_One_Rate(d,d->v[i],YES,tree);
907              }
908        }
909      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,b,d);
910      /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) {tree->both_sides = YES; Lk(tree); } */
911    }
912}
913
914//////////////////////////////////////////////////////////////
915//////////////////////////////////////////////////////////////
916
917void MCMC_One_Node_Rate(t_node *a, t_node *d, int traversal, t_tree *tree)
918{
919  t_edge *b;
920  int i;
921
922  b = NULL;
923  if(a == tree->n_root) b = tree->e_root;
924  else
925    For(i,3) if(d->v[i] == a) { b = d->b[i]; break; }
926
927  /* Only the LOG_RANDWALK move seems to work here. Change with caution then. */
928  tree->rates->br_do_updt[d->num] = YES;
929  MCMC_Single_Param_Generic(&(tree->rates->nd_r[d->num]),
930                            tree->rates->min_rate,
931                            tree->rates->max_rate,
932                            tree->mcmc->num_move_nd_r+d->num,
933                            &(tree->rates->c_lnL_rates),NULL,
934                            Wrap_Lk_Rates,NULL,
935                            tree->mcmc->move_type[tree->mcmc->num_move_nd_r+d->num],
936                            NO,NULL,tree,NULL);
937
938
939  Update_PMat_At_Given_Edge(b,tree);
940
941  if(traversal == YES)
942    {
943      if(d->tax == YES) return;
944      else
945        {
946          For(i,3)
947            if(d->v[i] != a && d->b[i] != tree->e_root)
948              {
949                MCMC_One_Node_Rate(d,d->v[i],YES,tree);
950              }
951        }
952    }
953}
954
955//////////////////////////////////////////////////////////////
956//////////////////////////////////////////////////////////////
957
958void MCMC_One_Time(t_node *a, t_node *d, int traversal, t_tree *tree)
959{
960  phydbl u;
961  phydbl t_min,t_max;
962  phydbl t1_cur, t1_new;
963  phydbl cur_lnL_data, new_lnL_data;
964  phydbl cur_lnL_rate, new_lnL_rate;
965  phydbl cur_lnL_time, new_lnL_time;
966  phydbl ratio,alpha;
967  t_edge *b1,*b2,*b3;
968  int    i;
969  phydbl t0,t2,t3;
970  t_node *v2,*v3;
971  phydbl K;
972  int move_num;
973 
974  if(d->tax) return; /* Won't change time at tip */
975
976  /* if(FABS(tree->rates->t_prior_min[d->num] - tree->rates->t_prior_max[d->num]) < 1.E-10) return; */
977
978  Record_Br_Len(tree);
979  RATES_Record_Rates(tree);
980  RATES_Record_Times(tree);
981 
982  move_num       = d->num-tree->n_otu+tree->mcmc->num_move_nd_t;
983  K              = tree->mcmc->tune_move[move_num];
984  cur_lnL_data   = tree->c_lnL;
985  cur_lnL_rate   = tree->rates->c_lnL_rates;
986  t1_cur         = tree->rates->nd_t[d->num];
987  new_lnL_data   = cur_lnL_data;
988  new_lnL_rate   = cur_lnL_rate;
989  ratio          = 0.0;
990  cur_lnL_time   = tree->rates->c_lnL_times;
991  new_lnL_time   = cur_lnL_time;
992
993  v2 = v3 = NULL;
994  For(i,3)
995    if((d->v[i] != a) && (d->b[i] != tree->e_root))
996      {
997        if(!v2) { v2 = d->v[i]; }
998        else    { v3 = d->v[i]; }
999      }
1000
1001
1002  b1 = NULL;
1003  if(a == tree->n_root) b1 = tree->e_root;
1004  else For(i,3) if(d->v[i] == a) { b1 = d->b[i]; break; }
1005
1006  b2 = b3 = NULL;
1007  For(i,3)
1008    if((d->v[i] != a) && (d->b[i] != tree->e_root))
1009      {
1010        if(!b2) { b2 = d->b[i]; }
1011        else    { b3 = d->b[i]; }
1012      }
1013 
1014  t0 = tree->rates->nd_t[a->num];
1015  t2 = tree->rates->nd_t[v2->num];
1016  t3 = tree->rates->nd_t[v3->num];
1017
1018
1019  /* t_min = MAX(t0,tree->rates->t_prior_min[d->num]);        */
1020  /* t_max = MIN(MIN(t2,t3),tree->rates->t_prior_max[d->num]);*/
1021
1022  t_min = t0;
1023  t_max = MIN(t2,t3);
1024
1025  t_min += tree->rates->min_dt;
1026  t_max -= tree->rates->min_dt;
1027
1028  if(t_min > t_max) 
1029    {
1030      PhyML_Printf("\n. t_min = %f t_max = %f",t_min,t_max);
1031      PhyML_Printf("\n. prior_min = %f prior_max = %f",tree->rates->t_prior_min[d->num],tree->rates->t_prior_max[d->num]);
1032      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1033      /* Exit("\n"); */
1034    }
1035 
1036  MCMC_Make_Move(&t1_cur,&t1_new,t_min,t_max,&ratio,K,tree->mcmc->move_type[move_num]);
1037 
1038
1039  if(t1_new > t_min && t1_new < t_max) 
1040    {
1041      tree->rates->nd_t[d->num] = t1_new;
1042
1043      new_lnL_time = TIMES_Lk_Times(tree);
1044      ratio += (new_lnL_time - cur_lnL_time);
1045
1046      if(isinf(new_lnL_time) == NO) // Proposed value of t is inside its boundary
1047        {
1048          /* Update branch lengths */
1049          tree->rates->br_do_updt[d->num]  = YES;
1050          tree->rates->br_do_updt[v2->num] = YES;
1051          tree->rates->br_do_updt[v3->num] = YES;
1052          RATES_Update_Cur_Bl(tree);
1053         
1054          new_lnL_rate = RATES_Lk_Rates(tree);
1055          ratio += (new_lnL_rate - cur_lnL_rate);
1056
1057          if(tree->mcmc->use_data)
1058            {
1059              if(tree->io->lk_approx == EXACT)
1060                {
1061                  Update_PMat_At_Given_Edge(b1,tree);
1062                  Update_PMat_At_Given_Edge(b2,tree);
1063                  Update_PMat_At_Given_Edge(b3,tree);
1064                  Update_P_Lk(tree,b1,d);
1065                }
1066              new_lnL_data = Lk(b1,tree);
1067             
1068              /* /\* !!!!!!!!!!!!!!!!!!!!1 *\/ */
1069              /* if(FABS(Lk(tree) - new_lnL_data) > 1.E-5) */
1070              /*   { */
1071              /*     PhyML_Printf("\n. b1->l->v=%f b2->l->v=%f b3->l->v=%f",b1->l->v,b2->l->v,b3->l->v); */
1072              /*     PhyML_Printf("\n. a->num=%d d->num=%d root=%d (%d %d)",a->num,d->num,a==tree->n_root,tree->n_root->v[2]->num,tree->n_root->v[1]->num); */
1073              /*     PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max); */
1074              /*     PhyML_Printf("\n. t1_new=%f t1_cur=%f",t1_new,t1_cur); */
1075              /*     PhyML_Printf("\n. %f %f",cur_lnL_data,tree->c_lnL); */
1076              /*     PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
1077              /*     Warn_And_Exit(""); */
1078              /*   } */
1079            }
1080         
1081          if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1082        }
1083     
1084      /* if(d->num == 7) */
1085      /*   { */
1086      /*     printf("\n. nd_t: %f %f new_lnL_rate: %f cur_lnL_rate: %f new_lnL_time: %f cur_lnL_time: %f ratio: %f", */
1087      /*            t1_cur, */
1088      /*            tree->rates->nd_t[d->num], */
1089      /*            new_lnL_rate, */
1090      /*            cur_lnL_rate, */
1091      /*            new_lnL_time, */
1092      /*            cur_lnL_time, */
1093      /*            ratio); */
1094      /*   } */
1095
1096         
1097      ratio = EXP(ratio);
1098      alpha = MIN(1.,ratio);
1099      u = Uni();
1100                   
1101      if(u > alpha) /* Reject */
1102        {
1103          //if(d -> num == 7) PhyML_Printf("\n. t_cur = %f t_rej = %f", t1_cur, t1_new);
1104          tree->rates->nd_t[d->num] = t1_cur;
1105          tree->c_lnL              = cur_lnL_data;
1106          tree->rates->c_lnL_rates = cur_lnL_rate;
1107          tree->rates->c_lnL_times = cur_lnL_time;
1108
1109          if(isinf(new_lnL_time) == NO)
1110            {
1111              Restore_Br_Len(tree);
1112              RATES_Reset_Rates(tree);
1113              RATES_Reset_Times(tree);
1114             
1115              if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) 
1116                {
1117                  Update_PMat_At_Given_Edge(b1,tree);
1118                  Update_PMat_At_Given_Edge(b2,tree);
1119                  Update_PMat_At_Given_Edge(b3,tree);
1120                  Update_P_Lk(tree,b1,d);
1121                }
1122            }
1123        }
1124      else
1125        {
1126          /* printf("\n. A new_lnL_data = %f cur_lnL_data = %f t_new=%f t_cur=%f tmin=%f tmax=%f", */
1127          /*     new_lnL_data,cur_lnL_data,t1_new,t1_cur,t_min,t_max); */
1128          tree->mcmc->acc_move[move_num]++;
1129        }
1130      if(t1_new < t0)
1131        {
1132          t1_new = t0+1.E-4;
1133          PhyML_Printf("\n");
1134          PhyML_Printf("\n== a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
1135          PhyML_Printf("\n== t0 = %f t1_new = %f",t0,t1_new);
1136          PhyML_Printf("\n== t_min=%f t_max=%f",t_min,t_max);
1137          PhyML_Printf("\n== (t1-t0)=%f (t2-t1)=%f",t1_cur-t0,t2-t1_cur);
1138          PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1139          /*       Exit("\n"); */
1140        }
1141      if(t1_new > MIN(t2,t3))
1142        {
1143          PhyML_Printf("\n");
1144          PhyML_Printf("\n== a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
1145          PhyML_Printf("\n== t0 = %f t1_new = %f t1 = %f t2 = %f t3 = %f MIN(t2,t3)=%f",t0,t1_new,t1_cur,t2,t3,MIN(t2,t3));
1146          PhyML_Printf("\n== t_min=%f t_max=%f",t_min,t_max);
1147          PhyML_Printf("\n== (t1-t0)=%f (t2-t1)=%f",t1_cur-t0,t2-t1_cur);
1148          PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1149          /*       Exit("\n"); */
1150        }
1151     
1152      if(isnan(t1_new))
1153        {
1154          PhyML_Printf("\n== run=%d",tree->mcmc->run);
1155          PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1156          /*       Exit("\n"); */
1157        }
1158    }
1159 
1160  tree->mcmc->run_move[move_num]++;
1161
1162  if(traversal == YES)
1163    {
1164      if(d->tax == YES) return;
1165      else
1166        For(i,3)
1167          if(d->v[i] != a && d->b[i] != tree->e_root)
1168            {
1169              if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,d->b[i],d);
1170              /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) {tree->both_sides = YES; Lk(tree); } */
1171              MCMC_One_Time(d,d->v[i],YES,tree);
1172            }
1173      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,b1,d);
1174      /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) {tree->both_sides = YES; Lk(tree); } */
1175    }       
1176}
1177
1178//////////////////////////////////////////////////////////////
1179//////////////////////////////////////////////////////////////
1180
1181void MCMC_Jump_Calibration(t_tree *tree)
1182{
1183
1184  phydbl u;
1185  phydbl *proba_distr, *times_partial_proba; 
1186  phydbl cur_lnL_data, new_lnL_data;
1187  phydbl cur_lnL_rate, new_lnL_rate;
1188  phydbl cur_lnL_time, new_lnL_time;
1189  phydbl cur_lnL_Hastings_ratio, new_lnL_Hastings_ratio;
1190  phydbl ratio,alpha;
1191  //t_edge *b1,*b2,*b3;
1192  int i, j, comb_num, result; //rnd_node
1193  //phydbl t0,t2,t3;
1194  //t_node *v2,*v3;
1195  //phydbl K;
1196  int move_num;
1197  int tot_num;
1198 
1199  tot_num = Number_Of_Comb(tree -> rates -> calib);
1200
1201  if(tot_num > 1)
1202    {
1203      times_partial_proba = tree -> rates -> times_partial_proba;
1204      proba_distr = (phydbl *)mCalloc(tot_num + 1, sizeof(phydbl)); 
1205 
1206      ///////////////////////////////////////////////////////////////////////////////////////////
1207      //for(i = tree -> n_otu; i < 2 * tree -> n_otu - 1; i++) printf("\n. '%f' '%f' \n", tree -> rates -> t_prior_min[i], tree -> rates -> t_prior_max[i]);
1208      //for(i = tree -> n_otu; i < 2 * tree -> n_otu - 1; i++) printf("\n. '%f' \n", tree -> rates -> nd_t[i]);
1209      ///////////////////////////////////////////////////////////////////////////////////////////
1210      Record_Br_Len(tree);
1211      RATES_Record_Rates(tree);
1212      RATES_Record_Times(tree);
1213      TIMES_Record_Prior_Times(tree);
1214     
1215      move_num = tree -> mcmc -> num_move_jump_calibration;
1216      //K = tree -> mcmc -> tune_move[move_num];
1217      cur_lnL_data = tree -> c_lnL;
1218      cur_lnL_rate = tree -> rates -> c_lnL_rates;
1219      cur_lnL_Hastings_ratio = tree -> rates -> c_lnL_Hastings_ratio;
1220      new_lnL_data = cur_lnL_data;
1221      new_lnL_rate = cur_lnL_rate;
1222      new_lnL_Hastings_ratio = cur_lnL_Hastings_ratio;
1223      ratio = 0.0;
1224      cur_lnL_time = tree -> rates -> c_lnL_times;
1225      new_lnL_time = cur_lnL_time;
1226     
1227      //printf("\n. [%d] \n", tot_num);
1228      //For(i, tot_num) printf("\n. [%f] \n", times_partial_proba[i]);
1229
1230      proba_distr[0]  = 0.0;
1231      for(i = 1; i < tot_num + 1; i++)
1232        {
1233          For(j, i)
1234            {
1235              proba_distr[i] = proba_distr[i] + times_partial_proba[j];
1236            }
1237          //printf("\n. [%f] \n", proba_distr[i]);
1238        }
1239      //For(i, tot_num + 1) printf("\n. [%f] \n", proba_distr[i]);
1240      u = Uni();
1241      comb_num = -1;
1242      For(i, tot_num) 
1243        {
1244          if(u > proba_distr[i] && (Are_Equal(proba_distr[i + 1], u, 1.E-10) || u < proba_distr[i + 1]))
1245            {
1246              //printf("\n. u [%f] [%f] [%f]\n", u, proba_distr[i], proba_distr[i + 1]);
1247              comb_num = i + 1;
1248            }
1249        }
1250      //printf("\n. [%d] \n", comb_num);
1251      //Exit("\n");     
1252     
1253      Set_Current_Calibration(comb_num - 1, tree);
1254      TIMES_Set_All_Node_Priors(tree);
1255     
1256      result = TRUE;
1257     
1258      Check_Node_Time(tree -> n_root, tree -> n_root -> v[1], &result, tree);
1259      Check_Node_Time(tree -> n_root, tree -> n_root -> v[2], &result, tree);
1260
1261      ///////////////////////////////////////////////////////////////////////////////////////////
1262      //if((comb_num - 1) == 1) for(i = tree -> n_otu; i < 2 * tree -> n_otu -1; i++) printf("\n. Node number:%d Min:%f Max:%f Cur.time:%f \n", i, tree -> rates -> t_prior_min[i], tree -> rates -> t_prior_max[i], tree -> rates -> nd_t[i]);
1263      //PhyML_Printf("\n. .......................................................................\n");
1264      ///////////////////////////////////////////////////////////////////////////////////////////
1265      //PhyML_Printf("\n. Result:%d \n", result); 
1266
1267      new_lnL_Hastings_ratio = 0.0;
1268      if(result != TRUE)
1269        {
1270          Update_Times_Down_Tree(tree -> n_root, tree -> n_root -> v[1], &new_lnL_Hastings_ratio, tree);
1271          Update_Times_Down_Tree(tree -> n_root, tree -> n_root -> v[2], &new_lnL_Hastings_ratio, tree);
1272          tree -> rates -> c_lnL_Hastings_ratio = new_lnL_Hastings_ratio;
1273        }
1274      else
1275        { 
1276          free(proba_distr);
1277          tree -> mcmc -> acc_move[move_num]++;   
1278          tree -> mcmc -> run_move[move_num]++;
1279          return; 
1280        }
1281      //PhyML_Printf("\n. Hastings Ratio:%f \n", cur_lnL_Hastings_ratio); 
1282      //PhyML_Printf("\n. Hastings Ratio:%f \n", new_lnL_Hastings_ratio); 
1283     
1284      result = TRUE;
1285     
1286      Check_Node_Time(tree -> n_root, tree -> n_root -> v[1], &result, tree);
1287      Check_Node_Time(tree -> n_root, tree -> n_root -> v[2], &result, tree);
1288     
1289      ///////////////////////////////////////////////////////////////////////////////////////////
1290      //for(i = tree -> n_otu; i < 2 * tree -> n_otu -1; i++) printf("\n. Node number:%d Min:%f Max:%f Cur.time:%f \n", i, tree -> rates -> t_prior_min[i], tree -> rates -> t_prior_max[i], tree -> rates -> nd_t[i]);
1291      //PhyML_Printf("\n. .......................................................................\n");
1292      ///////////////////////////////////////////////////////////////////////////////////////////
1293     
1294      if(result != TRUE)
1295        {
1296          PhyML_Printf("\n. ...................... OLD CALIBRATION.....................................\n");
1297          for(i = tree -> n_otu; i < 2 * tree -> n_otu -1; i++) printf("\n. Node number:[%d] Lower bound:[%f] Upper bound:[%f] Node time:[%f]. \n", i, tree -> rates -> t_prior_min_buff[i], tree -> rates -> t_prior_max_buff[i], tree -> rates -> buff_t[i]);
1298          PhyML_Printf("\n. ...........................................................................\n");
1299          PhyML_Printf("\n. ................. NEW PROPOSED CALIBRATION ................................\n");
1300          for(i = tree -> n_otu; i < 2 * tree -> n_otu -1; i++) printf("\n. Node number:[%d] Lower bound:[%f] Upper bound:[%f] Node time:[%f]. \n", i, tree -> rates -> t_prior_min[i], tree -> rates -> t_prior_max[i], tree -> rates -> nd_t[i]);
1301          PhyML_Printf("\n. ...........................................................................\n");
1302          PhyML_Printf("\n==You have a problem with calibration information.\n");
1303          PhyML_Printf("\n==Err in file %s at line %d\n",__FILE__,__LINE__);
1304          Exit("\n");
1305        }
1306         
1307         
1308      ///////////////////////////////////////////////////////////////////////////////////////////
1309      //Exit("\n");
1310      ///////////////////////////////////////////////////////////////////////////////////////////
1311     
1312      For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1313      RATES_Update_Cur_Bl(tree);
1314      if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1315     
1316      new_lnL_rate = RATES_Lk_Rates(tree);
1317      new_lnL_time = TIMES_Lk_Times(tree);
1318     
1319      /* Likelihood ratio */
1320      if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1321     
1322      /* Prior ratio */
1323      ratio += (new_lnL_rate - cur_lnL_rate);
1324      ratio += (new_lnL_time - cur_lnL_time);
1325      ratio += (cur_lnL_Hastings_ratio - new_lnL_Hastings_ratio); //Hastings ratio
1326     
1327      ratio = EXP(ratio);
1328      alpha = MIN(1.,ratio);
1329      u = Uni();
1330     
1331      if(u > alpha)
1332        {
1333          RATES_Reset_Times(tree);
1334          Restore_Br_Len(tree);
1335          TIMES_Reset_Prior_Times(tree);
1336          tree->c_lnL = cur_lnL_data;
1337          tree->rates->c_lnL_rates = cur_lnL_rate;
1338          tree->rates->c_lnL_times = cur_lnL_time;
1339          tree->rates->c_lnL_Hastings_ratio = cur_lnL_Hastings_ratio;
1340        }
1341      else
1342        {
1343          tree->mcmc->acc_move[move_num]++;
1344        }     
1345      tree->mcmc->run_move[move_num]++;
1346      free(proba_distr);
1347    }
1348  else 
1349    return;
1350}
1351
1352//////////////////////////////////////////////////////////////
1353//////////////////////////////////////////////////////////////
1354
1355void MCMC_Root_Time(t_tree *tree)
1356{
1357  phydbl u;
1358  phydbl t_min,t_max;
1359  phydbl t1_cur, t1_new;
1360  phydbl cur_lnL_data, new_lnL_data;
1361  phydbl cur_lnL_rate, new_lnL_rate;
1362  phydbl cur_lnL_time, new_lnL_time;
1363  phydbl ratio,alpha;
1364  t_edge *b1;
1365  phydbl t0,t2,t3;
1366  t_node *v2,*v3;
1367  phydbl K;
1368  int move_num;
1369  t_node *root;
1370
1371
1372  root = tree->n_root;
1373
1374  if(FABS(tree->rates->t_prior_min[root->num] - tree->rates->t_prior_max[root->num]) < 1.E-10) return;
1375
1376  Record_Br_Len(tree);
1377  RATES_Record_Rates(tree);
1378  RATES_Record_Times(tree);
1379 
1380  move_num       = root->num-tree->n_otu+tree->mcmc->num_move_nd_t;
1381  K              = tree->mcmc->tune_move[move_num];
1382  cur_lnL_data   = tree->c_lnL;
1383  cur_lnL_rate   = tree->rates->c_lnL_rates;
1384  t1_cur         = tree->rates->nd_t[root->num];
1385  new_lnL_data   = cur_lnL_data;
1386  new_lnL_rate   = cur_lnL_rate;
1387  ratio          = 0.0;
1388  cur_lnL_time   = tree->rates->c_lnL_times;
1389  new_lnL_time   = cur_lnL_time;
1390
1391 
1392  v2 = root->v[2];
1393  v3 = root->v[1];
1394
1395  b1 = tree->e_root;
1396 
1397  t0 = tree->rates->t_prior_min[root->num];
1398  t2 = tree->rates->nd_t[v2->num];
1399  t3 = tree->rates->nd_t[v3->num];
1400
1401  t_min = t0;
1402  /* t_max = MIN(MIN(t2,t3),tree->rates->t_prior_max[root->num]); */
1403  t_max = MIN(t2,t3);
1404
1405  t_min += tree->rates->min_dt;
1406  t_max -= tree->rates->min_dt;
1407
1408  if(t_min > t_max) 
1409    {
1410      PhyML_Printf("\n== t_min = %f t_max = %f",t_min,t_max);
1411      PhyML_Printf("\n== prior_min = %f prior_max = %f",tree->rates->t_prior_min[root->num],tree->rates->t_prior_max[root->num]);
1412      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1413      /* Exit("\n"); */
1414    }
1415
1416
1417  MCMC_Make_Move(&t1_cur,&t1_new,t_min,t_max,&ratio,K,tree->mcmc->move_type[move_num]);
1418 
1419  if(t1_new > t_min && t1_new < t_max) 
1420    {
1421      tree->rates->nd_t[root->num] = t1_new;
1422
1423      /* Update branch lengths */
1424      RATES_Update_Cur_Bl(tree);
1425
1426      if(tree->mcmc->use_data) new_lnL_data = Lk(b1,tree);
1427
1428      new_lnL_rate = RATES_Lk_Rates(tree);
1429      new_lnL_time = TIMES_Lk_Times(tree); 
1430
1431      if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1432      ratio += (new_lnL_rate - cur_lnL_rate);
1433      ratio += (new_lnL_time - cur_lnL_time);
1434
1435      ratio = EXP(ratio);
1436      alpha = MIN(1.,ratio);
1437      u = Uni();
1438         
1439      if(u > alpha) /* Reject */
1440        {
1441          tree->rates->nd_t[root->num] = t1_cur;
1442          tree->c_lnL              = cur_lnL_data;
1443          tree->rates->c_lnL_rates = cur_lnL_rate;
1444          tree->rates->c_lnL_times = cur_lnL_time;
1445          Restore_Br_Len(tree);
1446          RATES_Reset_Rates(tree);
1447          RATES_Reset_Times(tree);
1448        }
1449      else
1450        {
1451          tree->mcmc->acc_move[move_num]++;
1452        }
1453     
1454      if(t1_new < t0)
1455        {
1456          t1_new = t0+1.E-4;
1457          PhyML_Printf("\n");
1458          PhyML_Printf("\n. t0 = %f t1_new = %f",t0,t1_new);
1459          PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
1460          PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1_cur-t0,t2-t1_cur);
1461          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1462          /*       Exit("\n"); */
1463        }
1464      if(t1_new > MIN(t2,t3))
1465        {
1466          PhyML_Printf("\n");
1467          PhyML_Printf("\n. t0 = %f t1_new = %f t1 = %f t2 = %f t3 = %f MIN(t2,t3)=%f",t0,t1_new,t1_cur,t2,t3,MIN(t2,t3));
1468          PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
1469          PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1_cur-t0,t2-t1_cur);
1470          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1471          /*       Exit("\n"); */
1472        }
1473     
1474      if(isnan(t1_new))
1475        {
1476          PhyML_Printf("\n. run=%d",tree->mcmc->run);
1477          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1478          /*       Exit("\n"); */
1479        }
1480    }
1481 
1482  tree->mcmc->run_move[move_num]++;
1483
1484}
1485
1486//////////////////////////////////////////////////////////////
1487//////////////////////////////////////////////////////////////
1488
1489void MCMC_Tree_Height(t_tree *tree)
1490{
1491  int i;
1492  phydbl K,mult,u,alpha,ratio;
1493  phydbl cur_lnL_data,new_lnL_data;
1494  phydbl cur_lnL_rate,new_lnL_rate;
1495  phydbl cur_lnL_time,new_lnL_time;
1496  phydbl floor;
1497  int n_nodes;
1498 
1499  if(FABS(tree->rates->t_prior_max[tree->n_root->num] - tree->rates->t_prior_min[tree->n_root->num]) < 1.E-10) return;
1500
1501  RATES_Record_Times(tree);
1502  Record_Br_Len(tree);
1503
1504  K            = tree->mcmc->tune_move[tree->mcmc->num_move_tree_height];
1505  cur_lnL_data = tree->c_lnL;
1506  new_lnL_data = tree->c_lnL;
1507  ratio        = 0.0;
1508  cur_lnL_rate = tree->rates->c_lnL_rates;
1509  new_lnL_rate = tree->rates->c_lnL_rates;
1510  cur_lnL_time = tree->rates->c_lnL_times;
1511 
1512  u = Uni();
1513  mult = EXP(K*(u-0.5));
1514 
1515 /* WARNING: It must not be floor = tree->rates->t_prior_max[tree->n_root->num];
1516     floor is the maximum value a node height can take when one ignores the
1517     calibration nodes, i.e., floor is set by the height of the tips
1518  */
1519  /* floor = 0.0; */
1520  floor = tree->rates->t_floor[tree->n_root->num];
1521  /* floor = tree->rates->t_prior_max[tree->n_root->num]; */
1522
1523  Scale_Subtree_Height(tree->n_root,mult,floor,&n_nodes,tree);
1524
1525  /* For(i,2*tree->n_otu-1) */
1526  /*   { */
1527  /*     if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] || */
1528  /*     tree->rates->nd_t[i] < tree->rates->t_prior_min[i]) */
1529  /*    { */
1530  /*      RATES_Reset_Times(tree); */
1531  /*      Restore_Br_Len(tree); */
1532  /*      tree->mcmc->run_move[tree->mcmc->num_move_tree_height]++; */
1533  /*         return; */
1534  /*    } */
1535  /*   } */
1536 
1537  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1538  RATES_Update_Cur_Bl(tree);
1539  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1540
1541  new_lnL_rate = RATES_Lk_Rates(tree);
1542  new_lnL_time = TIMES_Lk_Times(tree);
1543
1544  /* The Hastings ratio is actually mult^(n) when changing the absolute
1545     node heights. When considering the relative heights, this ratio combined
1546     to the Jacobian for the change of variable ends up to being equal to mult.
1547  */
1548  /* ratio += LOG(mult); */
1549  ratio += (phydbl)(n_nodes)*LOG(mult);
1550
1551  /* Likelihood ratio */
1552  if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1553
1554  /* Prior ratio */
1555  ratio += (new_lnL_rate - cur_lnL_rate);
1556  ratio += (new_lnL_time - cur_lnL_time);
1557
1558  /* phydbl diff = (new_lnL_time - cur_lnL_time)+(phydbl)(n_nodes-1.)*LOG(mult); */
1559  /* printf("\n. diff_lk = %12f -(n-1)*log(mult) = %12f n_nodes=%d [%12f] log(mult)=%f", */
1560  /*     (new_lnL_time - cur_lnL_time), */
1561  /*     -(phydbl)(n_nodes-1.)*LOG(mult), */
1562  /*     n_nodes, */
1563  /*     diff, */
1564  /*     LOG(mult)); */
1565
1566  /* !!!!!!!!!!!! */
1567  /* ratio += LOG(Dexp(FABS(new_height-floor),1./10.) / Dexp(FABS(cur_height-floor),1./10.)); */
1568 
1569  ratio = EXP(ratio);
1570  alpha = MIN(1.,ratio);
1571  u = Uni();
1572
1573 
1574 
1575  if(u > alpha)
1576    {
1577      RATES_Reset_Times(tree);
1578      Restore_Br_Len(tree);
1579      tree->c_lnL = cur_lnL_data;
1580      tree->rates->c_lnL_rates = cur_lnL_rate;
1581      tree->rates->c_lnL_times = cur_lnL_time;
1582    }
1583  else
1584    {
1585      tree->mcmc->acc_move[tree->mcmc->num_move_tree_height]++;
1586      tree->mcmc->acc_move[tree->mcmc->num_move_nd_t+tree->n_root->num-tree->n_otu]++;
1587    }
1588
1589  tree->mcmc->run_move[tree->mcmc->num_move_tree_height]++;
1590  tree->mcmc->run_move[tree->mcmc->num_move_nd_t+tree->n_root->num-tree->n_otu]++;
1591}
1592
1593//////////////////////////////////////////////////////////////
1594//////////////////////////////////////////////////////////////
1595
1596void MCMC_Updown_T_Cr(t_tree *tree)
1597{
1598  /*! TO DO: make sure to change the values of clock_r across
1599    the different data partitions
1600  */
1601 
1602  int i;
1603  phydbl K,mult,u,alpha,ratio;
1604  phydbl cur_lnL_data,new_lnL_data;
1605  phydbl cur_lnL_rate,new_lnL_rate;
1606  phydbl cur_lnL_time,new_lnL_time;
1607  phydbl floor;
1608  int n_nodes;
1609
1610  /*! Check that sequences are isochronous. */
1611  For(i,tree->n_otu-1) if(!Are_Equal(tree->rates->nd_t[i+1],tree->rates->nd_t[i],1.E-10)) return; 
1612
1613  if(FABS(tree->rates->t_prior_max[tree->n_root->num] - tree->rates->t_prior_min[tree->n_root->num]) < 1.E-10) return;
1614
1615  RATES_Record_Times(tree);
1616  Record_Br_Len(tree);
1617
1618  K            = tree->mcmc->tune_move[tree->mcmc->num_move_updown_t_cr];
1619  cur_lnL_data = tree->c_lnL;
1620  new_lnL_data = tree->c_lnL;
1621  ratio        = 0.0;
1622  cur_lnL_rate = tree->rates->c_lnL_rates;
1623  new_lnL_rate = tree->rates->c_lnL_rates;
1624  cur_lnL_time = tree->rates->c_lnL_times;
1625
1626  u = Uni();
1627  mult = EXP(K*(u-0.5));
1628
1629  floor = 0.0;
1630
1631  Scale_Subtree_Height(tree->n_root,mult,floor,&n_nodes,tree);
1632
1633  For(i,2*tree->n_otu-1)
1634    {
1635      if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] ||
1636         tree->rates->nd_t[i] < tree->rates->t_prior_min[i])
1637        {
1638          RATES_Reset_Times(tree);
1639          Restore_Br_Len(tree);
1640          tree->mcmc->run_move[tree->mcmc->num_move_updown_t_cr]++;
1641          return;
1642        }
1643    }
1644
1645  if(RATES_Check_Node_Times(tree))
1646    {
1647      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1648      Warn_And_Exit("");
1649    }
1650
1651  tree->rates->clock_r /= mult;
1652  if(tree->rates->clock_r < tree->rates->min_clock || tree->rates->clock_r > tree->rates->max_clock)
1653    {
1654      tree->rates->clock_r *= mult;
1655      RATES_Reset_Times(tree);
1656      Restore_Br_Len(tree);
1657      tree->mcmc->run_move[tree->mcmc->num_move_updown_t_cr]++;
1658      return;
1659    }
1660 
1661  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1662  RATES_Update_Cur_Bl(tree);
1663  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1664
1665  new_lnL_rate = RATES_Lk_Rates(tree);
1666  new_lnL_time = TIMES_Lk_Times(tree);
1667
1668  /* The Hastings ratio is actually mult^(n) when changing the absolute
1669     node heights. When considering the relative heights, this ratio combined
1670     to the Jacobian for the change of variable ends up to being equal to mult.
1671  */
1672  ratio += (n_nodes - 1)*LOG(mult);
1673  /* ratio += -LOG(mult) + LOG(Dgamma(1./mult,1./K,K)/Dgamma(mult,1./K,K)); */
1674
1675  /* Likelihood ratio */
1676  if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1677
1678  /* Prior ratio */
1679  ratio += (new_lnL_rate - cur_lnL_rate);
1680  ratio += (new_lnL_time - cur_lnL_time);
1681
1682  /* !!!!!!!!!!!!1 */
1683  /* ratio += LOG(Dexp(FABS(new_height-floor),1./10.) / Dexp(FABS(cur_height-floor),1./10.)); */
1684 
1685  ratio = EXP(ratio);
1686  alpha = MIN(1.,ratio);
1687  u = Uni();
1688 
1689  /* printf("\n. t_old = %f t_new = %f cr_old = %f cr_new = %f", */
1690  /*        tree->rates->nd_t[tree->n_root->num]/mult, */
1691  /*        tree->rates->nd_t[tree->n_root->num], */
1692  /*        tree->rates->clock_r*mult, */
1693  /*        tree->rates->clock_r); */
1694
1695  if(u > alpha)
1696    {
1697      RATES_Reset_Times(tree);
1698      tree->rates->clock_r *= mult;
1699      Restore_Br_Len(tree);
1700      tree->c_lnL = cur_lnL_data;
1701      tree->rates->c_lnL_rates = cur_lnL_rate;
1702      tree->rates->c_lnL_times = cur_lnL_time;
1703    }
1704  else
1705    {
1706      tree->mcmc->acc_move[tree->mcmc->num_move_updown_t_cr]++;
1707    }
1708
1709  tree->mcmc->run_move[tree->mcmc->num_move_updown_t_cr]++;
1710}
1711
1712//////////////////////////////////////////////////////////////
1713//////////////////////////////////////////////////////////////
1714
1715void MCMC_Updown_T_Br(t_tree *tree)
1716{
1717 
1718  int i;
1719  phydbl K,mult,u,alpha,ratio;
1720  phydbl cur_lnL_data,new_lnL_data;
1721  phydbl cur_lnL_rate,new_lnL_rate;
1722  phydbl cur_lnL_time,new_lnL_time;
1723  phydbl floor;
1724  int n_nodes;
1725
1726  /*! Check that sequences are isochronous. */
1727  For(i,tree->n_otu-1) if(!Are_Equal(tree->rates->nd_t[i+1],tree->rates->nd_t[i],1.E-10)) return; 
1728
1729  if(FABS(tree->rates->t_prior_max[tree->n_root->num] - tree->rates->t_prior_min[tree->n_root->num]) < 1.E-10) return;
1730
1731  RATES_Record_Times(tree);
1732  Record_Br_Len(tree);
1733
1734  K            = tree->mcmc->tune_move[tree->mcmc->num_move_updown_t_br];
1735  cur_lnL_data = tree->c_lnL;
1736  new_lnL_data = tree->c_lnL;
1737  ratio        = 0.0;
1738  cur_lnL_rate = tree->rates->c_lnL_rates;
1739  new_lnL_rate = tree->rates->c_lnL_rates;
1740  cur_lnL_time = tree->rates->c_lnL_times;
1741
1742  u = Uni();
1743  mult = EXP(K*(u-0.5));
1744
1745
1746  floor = 0.0;
1747
1748  Scale_Subtree_Height(tree->n_root,mult,floor,&n_nodes,tree);
1749
1750  For(i,2*tree->n_otu-1)
1751    {
1752      if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] ||
1753         tree->rates->nd_t[i] < tree->rates->t_prior_min[i])
1754        {
1755          RATES_Reset_Times(tree);
1756          Restore_Br_Len(tree);
1757          tree->mcmc->run_move[tree->mcmc->num_move_updown_t_br]++;
1758          return;
1759        }
1760    }
1761
1762  if(RATES_Check_Node_Times(tree))
1763    {
1764      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
1765      Warn_And_Exit("");
1766    }
1767
1768  tree->rates->birth_rate /= mult;
1769
1770  if(tree->rates->birth_rate < tree->rates->birth_rate_min || tree->rates->birth_rate > tree->rates->birth_rate_max)
1771    {
1772      tree->rates->birth_rate *= mult;
1773      RATES_Reset_Times(tree);
1774      Restore_Br_Len(tree);
1775      tree->mcmc->run_move[tree->mcmc->num_move_updown_t_br]++;
1776      return;
1777    }
1778 
1779  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1780  RATES_Update_Cur_Bl(tree);
1781  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1782
1783  new_lnL_rate = RATES_Lk_Rates(tree);
1784  new_lnL_time = TIMES_Lk_Times(tree);
1785
1786  /* The Hastings ratio is actually mult^(n) when changing the absolute
1787     node heights. When considering the relative heights, this ratio combined
1788     to the Jacobian for the change of variable ends up to being equal to mult.
1789  */
1790  ratio += (n_nodes - 1)*LOG(mult);
1791  /* ratio += -LOG(mult) + LOG(Dgamma(1./mult,1./K,K)/Dgamma(mult,1./K,K)); */
1792
1793  /* Likelihood ratio */
1794  if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1795
1796  /* Prior ratio */
1797  ratio += (new_lnL_rate - cur_lnL_rate);
1798  ratio += (new_lnL_time - cur_lnL_time);
1799
1800  /* !!!!!!!!!!!!1 */
1801  /* ratio += LOG(Dexp(FABS(new_height-floor),1./10.) / Dexp(FABS(cur_height-floor),1./10.)); */
1802 
1803  ratio = EXP(ratio);
1804  alpha = MIN(1.,ratio);
1805  u = Uni();
1806 
1807  /* printf("\n. t_old = %f t_new = %f br_old = %f br_new = %f mult = %f K=%f", */
1808  /*        tree->rates->nd_t[tree->n_root->num]/mult, */
1809  /*        tree->rates->nd_t[tree->n_root->num], */
1810  /*        tree->rates->birth_rate*mult, */
1811  /*        tree->rates->birth_rate,mult,K); */
1812
1813  if(u > alpha)
1814    {
1815      RATES_Reset_Times(tree);
1816      tree->rates->birth_rate *= mult;
1817      Restore_Br_Len(tree);
1818      tree->c_lnL = cur_lnL_data;
1819      tree->rates->c_lnL_rates = cur_lnL_rate;
1820      tree->rates->c_lnL_times = cur_lnL_time;
1821    }
1822  else
1823    {
1824      tree->mcmc->acc_move[tree->mcmc->num_move_updown_t_br]++;
1825    }
1826
1827  tree->mcmc->run_move[tree->mcmc->num_move_updown_t_br]++;
1828}
1829
1830//////////////////////////////////////////////////////////////
1831//////////////////////////////////////////////////////////////
1832
1833void MCMC_Subtree_Height(t_tree *tree)
1834{
1835  int i;
1836  phydbl K,mult,u,alpha,ratio;
1837  phydbl cur_lnL_data,new_lnL_data;
1838  phydbl cur_lnL_rate,new_lnL_rate;
1839  phydbl cur_lnL_time,new_lnL_time;
1840  phydbl floor;
1841  int target;
1842  int n_nodes;
1843
1844  RATES_Record_Times(tree);
1845  Record_Br_Len(tree);
1846
1847  K = tree->mcmc->tune_move[tree->mcmc->num_move_subtree_height];
1848  cur_lnL_data = tree->c_lnL;
1849  new_lnL_data = tree->c_lnL;
1850  cur_lnL_rate = tree->rates->c_lnL_rates;
1851  new_lnL_rate = tree->rates->c_lnL_rates;
1852  ratio        = 0.0;
1853  cur_lnL_time = tree->rates->c_lnL_times;
1854
1855  u = Uni();
1856  mult = EXP(K*(u-0.5));
1857  /* mult = Rgamma(1./K,K); */
1858
1859  target = Rand_Int(tree->n_otu,2*tree->n_otu-3);
1860
1861  floor = tree->rates->t_floor[target];
1862
1863  if(tree->a_nodes[target] == tree->n_root)
1864    {
1865      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1866      Warn_And_Exit("");
1867    }
1868
1869  if(!Scale_Subtree_Height(tree->a_nodes[target],mult,floor,&n_nodes,tree))
1870    {
1871      RATES_Reset_Times(tree);
1872      Restore_Br_Len(tree);
1873      tree->mcmc->run_move[tree->mcmc->num_move_subtree_height]++;
1874      return;
1875    }
1876
1877 
1878  For(i,2*tree->n_otu-1)
1879    {
1880      if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] ||
1881         tree->rates->nd_t[i] < tree->rates->t_prior_min[i])
1882        {
1883          RATES_Reset_Times(tree);
1884          Restore_Br_Len(tree);
1885          tree->mcmc->run_move[tree->mcmc->num_move_subtree_height]++;
1886          return;
1887        }
1888    }
1889
1890     
1891  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1892  RATES_Update_Cur_Bl(tree);
1893  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1894
1895  new_lnL_rate = RATES_Lk_Rates(tree);
1896  new_lnL_time = TIMES_Lk_Times(tree);
1897
1898  /* The Hastings ratio here is mult^(n_nodes) and the ratio of the prior joint densities
1899     of the modified node heigths given the unchanged one is 1. This is different from the
1900     case where all the nodes, including the root node, are scaled.
1901  */
1902  ratio += (phydbl)(n_nodes)*LOG(mult);
1903
1904  /* Likelihood ratio */
1905  if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
1906
1907  /* Prior ratio */
1908  ratio += (new_lnL_rate - cur_lnL_rate);
1909  ratio += (new_lnL_time - cur_lnL_time);
1910
1911
1912  ratio = EXP(ratio);
1913  alpha = MIN(1.,ratio);
1914  u = Uni();
1915 
1916  if(u > alpha)
1917    {
1918      RATES_Reset_Times(tree);
1919      Restore_Br_Len(tree);
1920      tree->c_lnL = cur_lnL_data;
1921      tree->rates->c_lnL_rates = cur_lnL_rate;
1922      tree->rates->c_lnL_times = cur_lnL_time;
1923    }
1924  else
1925    {
1926      tree->mcmc->acc_move[tree->mcmc->num_move_subtree_height]++;
1927    }
1928
1929  tree->mcmc->run_move[tree->mcmc->num_move_subtree_height]++;
1930
1931}
1932
1933//////////////////////////////////////////////////////////////
1934//////////////////////////////////////////////////////////////
1935
1936
1937void MCMC_Tree_Rates(t_tree *tree)
1938{
1939  phydbl K,mult,u,alpha,ratio;
1940  phydbl cur_lnL_rate,new_lnL_rate;
1941  phydbl cur_lnL_data,new_lnL_data;
1942  int n_nodes;
1943  phydbl init_clock;
1944 
1945  if(tree->rates->model == STRICTCLOCK) return;
1946
1947  RATES_Record_Rates(tree);
1948  Record_Br_Len(tree);
1949
1950  K             = tree->mcmc->tune_move[tree->mcmc->num_move_tree_rates];
1951  cur_lnL_data  = tree->c_lnL;
1952  new_lnL_data  = tree->c_lnL;
1953  cur_lnL_rate  = tree->rates->c_lnL_rates;
1954  new_lnL_rate  = tree->rates->c_lnL_rates;
1955  init_clock    = tree->rates->clock_r;
1956  ratio         = 0.0;
1957
1958  u = Uni();
1959  mult = EXP(K*(u-0.5));
1960  /* mult = Rgamma(1./K,K); */
1961   
1962  /* Multiply branch rates (or add to log of rates) */
1963  if(!Scale_Subtree_Rates(tree->n_root,mult,&n_nodes,tree))
1964    {
1965      RATES_Reset_Rates(tree);
1966      Restore_Br_Len(tree);
1967      tree->mcmc->run_move[tree->mcmc->num_move_tree_rates]++;
1968      return;
1969    }
1970
1971  if(n_nodes != 2*tree->n_otu-2)
1972    {
1973      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
1974      Warn_And_Exit("");
1975    }
1976
1977  /* Divide clock_r */
1978  tree->rates->clock_r /= mult;
1979  if(tree->rates->clock_r < tree->rates->min_clock || tree->rates->clock_r > tree->rates->max_clock)
1980    {
1981      tree->rates->clock_r = init_clock;
1982      RATES_Reset_Rates(tree);
1983      Restore_Br_Len(tree);
1984      tree->mcmc->run_move[tree->mcmc->num_move_tree_rates]++;
1985      return;
1986    }
1987
1988/*   if(tree->rates->model == GUINDON) */
1989/*     { */
1990      int i;
1991      For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
1992      RATES_Update_Cur_Bl(tree);
1993      if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
1994/*     } */
1995
1996  new_lnL_rate = RATES_Lk_Rates(tree);
1997
1998
1999  /* Proposal ratio: 2n-2=> number of multiplications, 1=>number of divisions */
2000  ratio += (+(2*tree->n_otu-2)-1)*LOG(mult);
2001  /* ratio += (+(2*tree->n_otu-2)-1-2)*LOG(mult) + LOG(Dgamma(1./mult,1./K,K)/Dgamma(mult,1./K,K)); */
2002
2003  /* If modelling log of rates instead of rates */
2004  if(tree->rates->model_log_rates == YES) ratio -= (2*tree->n_otu-2)*LOG(mult);
2005
2006  /* Prior density ratio */
2007  ratio += (new_lnL_rate - cur_lnL_rate);
2008
2009  /* Likelihood density ratio */
2010  ratio += (new_lnL_data - cur_lnL_data);
2011
2012  ratio = EXP(ratio);
2013  alpha = MIN(1.,ratio);
2014  u = Uni();
2015 
2016/*   printf("\n. cur_lnL=%f new_lnL=%f ratio=%f mult=%f %f [%f %f]", */
2017/*       cur_lnL_rate,new_lnL_rate,ratio,mult,(+(2*tree->n_otu-2)-1-2)*LOG(mult) + LOG(Dgamma(1./mult,1./K,K)/Dgamma(mult,1./K,K)),new_lnL_data,cur_lnL_data); */
2018
2019  if(u > alpha)
2020    {
2021/*       PhyML_Printf("\n. Reject mult=%f",mult); */
2022      tree->rates->clock_r = init_clock;
2023      RATES_Reset_Rates(tree);
2024      Restore_Br_Len(tree);
2025      tree->rates->c_lnL_rates = cur_lnL_rate;
2026      tree->c_lnL        = cur_lnL_data;
2027    }
2028  else
2029    {
2030/*       PhyML_Printf("\n. Accept mult=%f",mult); */
2031      tree->mcmc->acc_move[tree->mcmc->num_move_tree_rates]++;
2032    }
2033
2034  tree->mcmc->run_move[tree->mcmc->num_move_tree_rates]++;
2035}
2036
2037//////////////////////////////////////////////////////////////
2038//////////////////////////////////////////////////////////////
2039
2040
2041void MCMC_Subtree_Rates(t_tree *tree)
2042{  phydbl K,mult,u,alpha,ratio;
2043  phydbl cur_lnL_rate,new_lnL_rate;
2044  phydbl cur_lnL_data,new_lnL_data;
2045  int target;
2046  int n_nodes;
2047 
2048  if(tree->rates->model == STRICTCLOCK) return;
2049
2050  RATES_Record_Rates(tree);
2051  Record_Br_Len(tree);
2052
2053  K             = tree->mcmc->tune_move[tree->mcmc->num_move_subtree_rates];
2054  cur_lnL_rate  = tree->rates->c_lnL_rates;
2055  new_lnL_rate  = cur_lnL_rate;
2056  cur_lnL_data  = tree->c_lnL;
2057  new_lnL_data  = cur_lnL_data;
2058  ratio         = 0.0;
2059
2060  u = Uni();
2061  mult = EXP(K*(u-0.5));
2062  /* mult = Rgamma(1./K,K); */
2063
2064  target = Rand_Int(tree->n_otu,2*tree->n_otu-3);
2065
2066  /* Multiply branch rates */
2067  if(!Scale_Subtree_Rates(tree->a_nodes[target],mult,&n_nodes,tree))
2068    {
2069      RATES_Reset_Rates(tree);
2070      Restore_Br_Len(tree);
2071      tree->mcmc->run_move[tree->mcmc->num_move_subtree_rates]++;
2072      return;
2073    }
2074
2075  new_lnL_rate = RATES_Lk_Rates(tree);
2076  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
2077
2078  /* Proposal ratio: 2n-2=> number of multiplications, 1=>number of divisions */
2079  ratio += (+n_nodes)*LOG(mult);
2080  /* ratio += (n_nodes-2)*LOG(mult) + LOG(Dgamma(1./mult,1./K,K)/Dgamma(mult,1./K,K)); */
2081 
2082  /* If modelling log of rates instead of rates */
2083  if(tree->rates->model_log_rates == YES) ratio -= (n_nodes)*LOG(mult);
2084
2085  /* Prior density ratio */
2086  ratio += (new_lnL_rate - cur_lnL_rate);
2087
2088  /* Likelihood density ratio */
2089  ratio += (new_lnL_data - cur_lnL_data);
2090
2091
2092  ratio = EXP(ratio);
2093  alpha = MIN(1.,ratio);
2094  u = Uni();
2095 
2096  if(u > alpha)
2097    {
2098      RATES_Reset_Rates(tree);
2099      Restore_Br_Len(tree);
2100      tree->rates->c_lnL_rates = cur_lnL_rate;
2101      tree->c_lnL        = cur_lnL_data;
2102    }
2103  else
2104    {
2105      tree->mcmc->acc_move[tree->mcmc->num_move_subtree_rates]++;
2106    }
2107
2108  tree->mcmc->run_move[tree->mcmc->num_move_subtree_rates]++;
2109}
2110
2111//////////////////////////////////////////////////////////////
2112//////////////////////////////////////////////////////////////
2113
2114
2115void MCMC_Swing(t_tree *tree)
2116{
2117  int i;
2118  phydbl K,mult,u,alpha,ratio;
2119  phydbl cur_lnL_data,new_lnL_data;
2120  phydbl cur_lnL_rate,new_lnL_rate;
2121
2122  if(tree->rates->model == STRICTCLOCK) return;
2123
2124  RATES_Record_Times(tree);
2125  RATES_Record_Rates(tree);
2126  Record_Br_Len(tree);
2127
2128  K             = 3.;
2129  cur_lnL_data  = tree->c_lnL;
2130  new_lnL_data  = cur_lnL_data;
2131  cur_lnL_rate  = tree->rates->c_lnL_rates;
2132  new_lnL_rate  = cur_lnL_rate;
2133  ratio         = 0.0;
2134
2135  u = Uni();
2136  /* mult = EXP(K*(u-0.5)); */
2137  mult = u*(K - 1./K) + 1./K;
2138
2139
2140  For(i,2*tree->n_otu-1)
2141    {
2142      if(tree->a_nodes[i]->tax == NO) 
2143        {
2144          tree->rates->nd_t[i] *= mult;
2145        }
2146
2147      if(tree->rates->nd_t[i] > tree->rates->t_prior_max[i] ||
2148         tree->rates->nd_t[i] < tree->rates->t_prior_min[i])
2149        {
2150          RATES_Reset_Times(tree);
2151          Restore_Br_Len(tree);
2152          return;
2153        }
2154    }
2155
2156  For(i,2*tree->n_otu-2)
2157    {
2158      tree->rates->br_r[i] /= mult;
2159
2160      if(tree->rates->br_r[i] > tree->rates->max_rate ||
2161         tree->rates->br_r[i] < tree->rates->min_rate)
2162        {
2163          RATES_Reset_Times(tree);
2164          RATES_Reset_Rates(tree);
2165          Restore_Br_Len(tree);
2166          return;
2167        }
2168    }
2169
2170  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
2171  RATES_Update_Cur_Bl(tree);
2172  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
2173
2174  new_lnL_rate = RATES_Lk_Rates(tree);
2175
2176  ratio += (-(tree->n_otu-1.)-2.)*LOG(mult);
2177  ratio += (new_lnL_rate - cur_lnL_rate);
2178  if(tree->mcmc->use_data) ratio += (new_lnL_data - cur_lnL_data);
2179
2180  ratio = EXP(ratio);
2181  alpha = MIN(1.,ratio);
2182  u = Uni();
2183
2184  if(u > alpha)
2185    {
2186      RATES_Reset_Times(tree);
2187      RATES_Reset_Rates(tree);
2188      Restore_Br_Len(tree);
2189      tree->c_lnL = cur_lnL_data;
2190      tree->rates->c_lnL_rates = cur_lnL_rate;
2191      /* printf("\n. Reject %8f",mult); */
2192    }
2193  else
2194    {
2195      /* printf("\n. Accept %8f",mult); */
2196    }
2197
2198
2199}
2200
2201//////////////////////////////////////////////////////////////
2202//////////////////////////////////////////////////////////////
2203
2204
2205void MCMC_Updown_Nu_Cr(t_tree *tree)
2206{
2207  phydbl K,mult,u,alpha,ratio;
2208  phydbl cur_lnL_rate,new_lnL_rate;
2209  phydbl cur_lnL_data,new_lnL_data;
2210  int i;
2211 
2212  RATES_Record_Rates(tree);
2213  Record_Br_Len(tree);
2214
2215  K             = tree->mcmc->tune_move[tree->mcmc->num_move_updown_nu_cr];
2216  cur_lnL_data  = tree->c_lnL;
2217  new_lnL_data  = tree->c_lnL;
2218  cur_lnL_rate  = tree->rates->c_lnL_rates;
2219  new_lnL_rate  = tree->rates->c_lnL_rates;
2220 
2221  u = Uni();
2222  mult = EXP(K*(u-0.5));
2223
2224  /* Multiply branch rates */
2225  /* if(!Scale_Subtree_Rates(tree->n_root,mult,&n_nodes,tree)) */
2226  /*   { */
2227  /*     RATES_Reset_Rates(tree); */
2228  /*     Restore_Br_Len(tree); */
2229  /*     tree->mcmc->run_move[tree->mcmc->num_move_updown_nu_cr]++; */
2230  /*     return; */
2231  /*   } */
2232
2233  tree->rates->clock_r /= mult;
2234  if(tree->rates->clock_r < tree->rates->min_clock || tree->rates->clock_r > tree->rates->max_clock)
2235    {
2236      tree->rates->clock_r *= mult;
2237      RATES_Reset_Rates(tree);
2238      Restore_Br_Len(tree);
2239      tree->mcmc->run_move[tree->mcmc->num_move_updown_nu_cr]++;
2240      return;
2241    }
2242
2243  tree->rates->nu *= mult;
2244  if(tree->rates->nu < tree->rates->min_nu || tree->rates->nu > tree->rates->max_nu)
2245    {
2246      tree->rates->nu      /= mult;
2247      tree->rates->clock_r *= mult;
2248      RATES_Reset_Rates(tree);
2249      Restore_Br_Len(tree);
2250      tree->mcmc->run_move[tree->mcmc->num_move_updown_nu_cr]++;
2251      return;
2252    }
2253 
2254  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
2255  RATES_Update_Cur_Bl(tree);
2256  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
2257
2258  new_lnL_rate = RATES_Lk_Rates(tree);
2259
2260  ratio = 0.0;
2261  /* Proposal ratio: 2n-2=> number of multiplications, 1=>number of divisions */
2262  /* ratio += n_nodes*LOG(mult); /\* (1-1)*LOG(mult); *\/ */
2263  ratio += 0.0*LOG(mult); /* (1-1)*LOG(mult); */
2264  /* Prior density ratio */
2265  ratio += (new_lnL_rate - cur_lnL_rate);
2266  /* Likelihood density ratio */
2267  ratio += (new_lnL_data - cur_lnL_data);
2268
2269  ratio = EXP(ratio);
2270  alpha = MIN(1.,ratio);
2271  u = Uni();
2272 
2273  if(u > alpha)
2274    {
2275      tree->rates->clock_r *= mult;
2276      tree->rates->nu /= mult;
2277      RATES_Reset_Rates(tree);
2278      Restore_Br_Len(tree);
2279      tree->rates->c_lnL_rates = cur_lnL_rate;
2280      tree->c_lnL        = cur_lnL_data;
2281    }
2282  else
2283    {
2284      tree->mcmc->acc_move[tree->mcmc->num_move_updown_nu_cr]++;
2285    }
2286
2287  tree->mcmc->run_move[tree->mcmc->num_move_updown_nu_cr]++;
2288}
2289
2290//////////////////////////////////////////////////////////////
2291//////////////////////////////////////////////////////////////
2292
2293
2294void MCMC_Print_Param_Stdin(t_mcmc *mcmc, t_tree *tree)
2295{
2296  time_t cur_time;
2297  phydbl min;
2298  int i;
2299  time(&cur_time);
2300 
2301  min = MDBL_MAX;
2302  For(i,tree->n_otu-1)
2303    {
2304      /*       printf("\n. %d %f %f %f %f",i, */
2305      /*             tree->mcmc->new_param_val[tree->mcmc->num_move_nd_t+i], */
2306      /*             tree->mcmc->old_param_val[tree->mcmc->num_move_nd_t+i], */
2307      /*             tree->mcmc->ess[tree->mcmc->num_move_nd_t+i], */
2308      /*             tree->mcmc->sum_val[tree->mcmc->num_move_nd_t+i]); */
2309
2310      if(tree->mcmc->ess[tree->mcmc->num_move_nd_t+i] < min)
2311        min = tree->mcmc->ess[tree->mcmc->num_move_nd_t+i];
2312    }
2313
2314
2315  if(mcmc->run == 1)
2316    {
2317      PhyML_Printf("\n\n");
2318      PhyML_Printf("%9s","Run");
2319      PhyML_Printf("  %5s","Time");
2320      PhyML_Printf("  %10s","Likelihood");
2321      PhyML_Printf("  %10s","Prior");
2322      PhyML_Printf("  %19s","SubstRate[ ESS ]");
2323      PhyML_Printf("  %17s","TreeHeight[ ESS ]");   
2324      if(tree->rates->model == THORNE || tree->rates->model == GUINDON) PhyML_Printf("  %16s","AutoCor[ ESS ]");   
2325      else PhyML_Printf("  %16s","RateVar[ ESS ]");
2326      PhyML_Printf("  %15s","BirthR[ ESS ]");
2327      PhyML_Printf("  %8s","MinESS");   
2328    }
2329
2330  if((cur_time - mcmc->t_last_print) >  mcmc->print_every)
2331    {
2332      mcmc->t_last_print = cur_time;
2333      PhyML_Printf("\n");
2334      PhyML_Printf("%9d",tree->mcmc->run);
2335      PhyML_Printf("  %5d",(int)(cur_time-mcmc->t_beg));
2336      PhyML_Printf("  %10.2f",tree->c_lnL);
2337      PhyML_Printf("  %10.2f",(tree->rates ? tree->rates->c_lnL_rates+tree->rates->c_lnL_times : +1));
2338      PhyML_Printf("  %12.6f[%5.0f]",RATES_Average_Substitution_Rate(tree),tree->mcmc->ess[tree->mcmc->num_move_clock_r]);
2339      /* PhyML_Printf("\t%12.6f[%5.0f]",tree->rates->clock_r,tree->mcmc->ess[tree->mcmc->num_move_clock_r]); */
2340      PhyML_Printf("  %10.1f[%5.0f]",
2341                   (tree->rates ? tree->rates->nd_t[tree->n_root->num] : -1.),
2342                   tree->mcmc->ess[tree->mcmc->num_move_nd_t+tree->n_root->num-tree->n_otu]);
2343      PhyML_Printf("  %9f[%5.0f]",
2344                   (tree->rates ? tree->rates->nu : -1.),
2345                   tree->mcmc->ess[tree->mcmc->num_move_nu]);
2346      PhyML_Printf("  %8f[%5.0f]",
2347                   (tree->rates ? tree->rates->birth_rate : -1.),
2348                   tree->mcmc->ess[tree->mcmc->num_move_birth_rate]);
2349      PhyML_Printf("  %8.0f",min);
2350    }
2351}
2352
2353//////////////////////////////////////////////////////////////
2354//////////////////////////////////////////////////////////////
2355
2356
2357void MCMC_Print_Param(t_mcmc *mcmc, t_tree *tree)
2358{
2359  int i;
2360  FILE *fp;
2361  char *s;
2362  int orig_approx;
2363  phydbl orig_lnL;
2364  char *s_tree;
2365 
2366  if(tree->mcmc->run > mcmc->chain_len) return;
2367
2368  s = (char *)mCalloc(100,sizeof(char));
2369
2370  fp = mcmc->out_fp_stats;
2371
2372/*   if(tree->mcmc->run == 0) */
2373/*     { */
2374/*       PhyML_Fprintf(stdout," ["); */
2375/*       fflush(NULL); */
2376/*     } */
2377 
2378/*   if(!(mcmc->run%(mcmc->chain_len/10))) */
2379/*     { */
2380/*       PhyML_Fprintf(stdout,"."); */
2381/*       fflush(NULL); */
2382/*     } */
2383
2384/*   if(tree->mcmc->run == mcmc->chain_len) */
2385/*     { */
2386/*       PhyML_Fprintf(stdout,"]"); */
2387/*       fflush(NULL); */
2388/*     } */
2389
2390
2391/*   MCMC_Print_Means(mcmc,tree); */
2392/*   MCMC_Print_Last(mcmc,tree); */
2393
2394
2395
2396  if(!(mcmc->run%mcmc->sample_interval)) 
2397    {
2398
2399      MCMC_Copy_To_New_Param_Val(tree->mcmc,tree);
2400
2401      For(i,tree->mcmc->n_moves)
2402        {
2403          if((tree->mcmc->acc_rate[i] > .1) && (tree->mcmc->start_ess[i] == NO)) tree->mcmc->start_ess[i] = YES;
2404          if(tree->mcmc->start_ess[i] == YES) MCMC_Update_Effective_Sample_Size(i,tree->mcmc,tree);
2405          if(tree->mcmc->run > (int) tree->mcmc->chain_len * 0.1) tree->mcmc->adjust_tuning[i] = NO;
2406        }
2407
2408      if(tree->mcmc->run == 0)
2409        {
2410          time(&(mcmc->t_beg));
2411          time(&(mcmc->t_last_print));
2412
2413          PhyML_Fprintf(fp,"# Random seed: %d",tree->io->r_seed);
2414          PhyML_Fprintf(fp,"\n");
2415          PhyML_Fprintf(fp,"Run\t");
2416/*        PhyML_Fprintf(fp,"Time\t"); */
2417          /* PhyML_Fprintf(fp,"MeanRate\t"); */
2418/*        PhyML_Fprintf(fp,"NormFact\t"); */
2419
2420          For(i,mcmc->n_moves)
2421            {
2422              strcpy(s,"Acc.");
2423              PhyML_Fprintf(fp,"%s%d\t",strcat(s,mcmc->move_name[i]),i);
2424            }
2425
2426          For(i,mcmc->n_moves)
2427            {
2428              strcpy(s,"Tune.");
2429              PhyML_Fprintf(fp,"%s%d\t",strcat(s,mcmc->move_name[i]),i);
2430            }
2431
2432          /* For(i,mcmc->n_moves) */
2433          /*   { */
2434          /*     strcpy(s,"Run."); */
2435          /*     PhyML_Fprintf(fp,"%s\t",strcat(s,mcmc->move_name[i])); */
2436          /*   } */
2437
2438          PhyML_Fprintf(fp,"LnLike[Exact]\t");
2439          PhyML_Fprintf(fp,"LnLike[Approx]\t");
2440          PhyML_Fprintf(fp,"LnRates\t");
2441          PhyML_Fprintf(fp,"LnTimes\t");
2442          PhyML_Fprintf(fp,"LnPosterior\t");
2443          PhyML_Fprintf(fp,"ClockRate\t");
2444          PhyML_Fprintf(fp,"EvolRate\t");
2445          PhyML_Fprintf(fp,"Nu\t");
2446          PhyML_Fprintf(fp,"BirthRate\t");
2447          PhyML_Fprintf(fp,"TsTv\t");
2448
2449
2450          if(tree->mod->ras->n_catg > 1)
2451            {
2452              if(tree->mod->ras->free_mixt_rates == NO) PhyML_Fprintf(fp,"Alpha\t");
2453              else
2454                {
2455                  For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"p%d\t",i);
2456                  For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"r%d\t",i);
2457                }
2458            }
2459
2460
2461          if(tree->mod->m4mod->n_h > 1 && tree->mod->use_m4mod == YES)
2462            {
2463              For(i,tree->mod->m4mod->n_h) PhyML_Fprintf(fp,"cov_p%d\t",i);
2464              For(i,tree->mod->m4mod->n_h) PhyML_Fprintf(fp,"cov_r%d\t",i);
2465              PhyML_Fprintf(fp,"cov_switch\t");
2466            }
2467
2468
2469          if(fp != stdout)
2470            {
2471              for(i=tree->n_otu;i<2*tree->n_otu-1;i++)
2472                {
2473                  if(tree->a_nodes[i] == tree->n_root->v[2])
2474                    PhyML_Fprintf(fp,"T%d%s\t",i,"[LeftRoot]");
2475                  else if(tree->a_nodes[i] == tree->n_root->v[1])
2476                    PhyML_Fprintf(fp,"T%d%s\t",i,"[RightRoot]");
2477                  else if(tree->a_nodes[i] == tree->n_root)
2478                    PhyML_Fprintf(fp,"T%d%s\t",i,"[Root]");
2479                  else
2480                    PhyML_Fprintf(fp,"T%d[%d]\t",i,tree->a_nodes[i]->anc->num);
2481                }
2482            }
2483
2484
2485/*        if(fp != stdout) */
2486/*          { */
2487/*            For(i,2*tree->n_otu-1) */
2488/*              { */
2489/*                if(tree->a_nodes[i] == tree->n_root->v[2]) */
2490/*                  PhyML_Fprintf(fp,"R%d[LeftRoot]\t",i); */
2491/*                else if(tree->a_nodes[i] == tree->n_root->v[1]) */
2492/*                  PhyML_Fprintf(fp,"R%d[RightRoot]\t",i); */
2493/*                else if(tree->a_nodes[i] != tree->n_root) */
2494/*                  PhyML_Fprintf(fp," R%d[%d]\t",i,tree->a_nodes[i]->anc->num); */
2495/*                else */
2496/*                  PhyML_Fprintf(fp," R%d[Root]\t",i); */
2497
2498/* /\*            PhyML_Fprintf(fp," R%d[%f]\t",i,tree->rates->mean_l[i]); *\/ */
2499/*              } */
2500/*          } */
2501
2502
2503          if(fp != stdout)
2504            {
2505              For(i,2*tree->n_otu-2)
2506                {
2507                  if(tree->a_nodes[i] == tree->n_root->v[2])
2508                    PhyML_Fprintf(fp,"B%d[LeftRoot]\t",i);
2509                  else if(tree->a_nodes[i] == tree->n_root->v[1])
2510                    PhyML_Fprintf(fp,"B%d[RightRoot]\t",i);
2511                  else
2512                    PhyML_Fprintf(fp," B%d[%d]\t",i,tree->a_nodes[i]->anc->num);
2513
2514/*                PhyML_Fprintf(fp," R%d[%f]\t",i,tree->rates->mean_l[i]); */
2515                }
2516            }
2517         
2518
2519      /*          if(fp != stdout) */
2520      /*            { */
2521      /*              For(i,2*tree->n_otu-2) */
2522      /*                { */
2523      /*                  if(tree->a_nodes[i] == tree->n_root->v[2]) */
2524      /*                    PhyML_Fprintf(fp,"G%d[LeftRoot]\t",i); */
2525      /*                  else if(tree->a_nodes[i] == tree->n_root->v[1]) */
2526      /*                    PhyML_Fprintf(fp,"G%d[RightRoot]\t",i); */
2527      /*                  else */
2528      /*                    PhyML_Fprintf(fp," G%d[%f]\t",i,tree->rates->ml_l[i]); */
2529
2530
2531      /* /\*              PhyML_Fprintf(fp," R%d[%f]\t",i,tree->rates->mean_l[i]); *\/ */
2532      /*                } */
2533      /*            } */
2534
2535
2536          /* if(fp != stdout) */
2537          /*   { */
2538          /*     For(i,2*tree->n_otu-3) */
2539          /*       { */
2540          /*         if(tree->a_edges[i] == tree->e_root) */
2541          /*           PhyML_Fprintf(fp,"*L[%f]%d\t",i,tree->rates->u_ml_l[i]); */
2542          /*         else */
2543          /*           PhyML_Fprintf(fp," L[%f]%d\t",i,tree->rates->u_ml_l[i]); */
2544          /*       } */
2545          /*   } */
2546
2547
2548          PhyML_Fprintf(mcmc->out_fp_trees,"#NEXUS\n");
2549          PhyML_Fprintf(mcmc->out_fp_trees,"BEGIN TREES;\n");
2550          PhyML_Fprintf(mcmc->out_fp_trees,"\tTRANSLATE\n");
2551          For(i,tree->n_otu-1) PhyML_Fprintf(mcmc->out_fp_trees,"\t%3d\t%s,\n",tree->a_nodes[i]->num+1,tree->a_nodes[i]->name);
2552          PhyML_Fprintf(mcmc->out_fp_trees,"\t%3d\t%s;\n",tree->a_nodes[i]->num+1,tree->a_nodes[i]->name);
2553          tree->write_tax_names = NO;
2554        }
2555
2556      PhyML_Fprintf(fp,"\n");
2557
2558      PhyML_Fprintf(fp,"%6d\t",tree->mcmc->run);
2559
2560/*       time(&mcmc->t_cur); */
2561/*       PhyML_Fprintf(fp,"%6d\t",(int)(mcmc->t_cur-mcmc->t_beg)); */
2562     
2563/*       RATES_Update_Cur_Bl(tree); */
2564/*       PhyML_Fprintf(fp,"%f\t",RATES_Check_Mean_Rates(tree)); */
2565
2566/*       PhyML_Fprintf(fp,"%f\t",tree->rates->norm_fact); */
2567      For(i,tree->mcmc->n_moves) PhyML_Fprintf(fp,"%f\t",tree->mcmc->acc_rate[i]);
2568      For(i,tree->mcmc->n_moves) PhyML_Fprintf(fp,"%f\t",(phydbl)(tree->mcmc->tune_move[i]));
2569/*       For(i,tree->mcmc->n_moves) PhyML_Fprintf(fp,"%d\t",(int)(tree->mcmc->run_move[i])); */
2570
2571      orig_approx = tree->io->lk_approx;
2572      orig_lnL = tree->c_lnL;
2573      tree->io->lk_approx = EXACT;
2574      if(tree->mcmc->use_data)  Lk(NULL,tree);  else tree->c_lnL = 0.0;
2575      PhyML_Fprintf(fp,"%.1f\t",tree->c_lnL);
2576      tree->io->lk_approx = NORMAL;
2577      tree->c_lnL = 0.0;
2578      if(tree->mcmc->use_data)  Lk(NULL,tree);  else tree->c_lnL = 0.0;
2579      PhyML_Fprintf(fp,"%.1f\t",tree->c_lnL);
2580      tree->io->lk_approx = orig_approx;
2581      tree->c_lnL = orig_lnL;
2582
2583/*       PhyML_Fprintf(fp,"0\t0\t"); */
2584
2585      PhyML_Fprintf(fp,"%G\t",tree->rates->c_lnL_rates);
2586      PhyML_Fprintf(fp,"%G\t",tree->rates->c_lnL_times);
2587      PhyML_Fprintf(fp,"%G\t",tree->c_lnL+tree->rates->c_lnL_rates+tree->rates->c_lnL_times);
2588      PhyML_Fprintf(fp,"%G\t",tree->rates->clock_r);
2589      PhyML_Fprintf(fp,"%G\t",RATES_Average_Substitution_Rate(tree));
2590      PhyML_Fprintf(fp,"%G\t",tree->rates->nu);
2591      PhyML_Fprintf(fp,"%G\t",tree->rates->birth_rate);
2592      PhyML_Fprintf(fp,"%G\t",tree->mod->kappa->v);
2593     
2594      if(tree->mod->ras->n_catg > 1)
2595        {
2596          if(tree->mod->ras->free_mixt_rates == NO)
2597            PhyML_Fprintf(fp,"%G\t",tree->mod->ras->alpha->v);
2598          else
2599            {
2600              For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"%G\t",tree->mod->ras->gamma_r_proba->v[i]);
2601              For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"%G\t",tree->mod->ras->gamma_rr->v[i]);
2602              /* For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"%G\t",tree->mod->ras->gamma_r_proba_unscaled[i]); */
2603              /* For(i,tree->mod->ras->n_catg) PhyML_Fprintf(fp,"%G\t",tree->mod->ras->gamma_rr_unscaled[i]); */
2604            }
2605        }
2606
2607      if(tree->mod->m4mod->n_h > 1 && tree->mod->use_m4mod == YES)
2608        {
2609          For(i,tree->mod->m4mod->n_h) PhyML_Fprintf(fp,"%G\t",tree->mod->m4mod->h_fq[i]);
2610          For(i,tree->mod->m4mod->n_h) PhyML_Fprintf(fp,"%G\t",tree->mod->m4mod->multipl[i]);
2611          PhyML_Fprintf(fp,"%G\t",tree->mod->m4mod->delta);
2612        }
2613
2614
2615      char *format = (char *)mCalloc(100,sizeof(char));
2616      sprintf(format,"%%.%df\t",tree->mcmc->nd_t_digits);
2617      for(i=tree->n_otu;i<2*tree->n_otu-1;i++) PhyML_Fprintf(fp,format,tree->rates->nd_t[i]);
2618      Free(format);
2619
2620      /* for(i=0;i<2*tree->n_otu-1;i++) PhyML_Fprintf(fp,"%.4f\t",LOG(tree->rates->nd_r[i])); */
2621     
2622      // Average rate along edges: length divided by elapsed time
2623      For(i,2*tree->n_otu-2)
2624        PhyML_Fprintf(fp,"%.4f\t",
2625                      tree->rates->cur_l[i]/(tree->rates->nd_t[tree->a_nodes[i]->num] - tree->rates->nd_t[tree->a_nodes[i]->anc->num]));
2626
2627      /* fp_pred = fopen("predict.txt","a");       */
2628      /* for(i=0;i<2*tree->n_otu-2;i++)  */
2629      /*        PhyML_Fprintf(fp_pred,"B%d\t%12f\t%12f\t%4d\n",i,EXP(tree->rates->br_r[i]),tree->rates->nd_t[i],tree->rates->has_survived[i]); */
2630      /* fclose(fp_pred); */
2631
2632
2633      /* phydbl p,sd,mean; */
2634     
2635      /* for(i=0;i<2*tree->n_otu-2;i++) */
2636      /*        { */
2637      /*          sd = tree->rates->nu * (tree->rates->nd_t[i] - tree->rates->nd_t[tree->a_nodes[i]->anc->num]); */
2638      /*          mean = tree->rates->br_r[tree->a_nodes[i]->anc->num] - .5*sd*sd; */
2639      /*          p = Pnorm(tree->rates->br_r[i],mean,sd); */
2640      /*          PhyML_Fprintf(fp,"%f\t",p); */
2641      /*          tree->rates->mean_r[i] += p; */
2642      /*        } */
2643
2644      /* for(i=0;i<2*tree->n_otu-2;i++) PhyML_Fprintf(fp,"%.4f\t",tree->rates->cur_gamma_prior_mean[i]); */
2645      /* if(fp != stdout) for(i=tree->n_otu;i<2*tree->n_otu-1;i++) PhyML_Fprintf(fp,"%G\t",tree->rates->t_prior[i]); */
2646/*       For(i,2*tree->n_otu-3) PhyML_Fprintf(fp,"%f\t",EXP(tree->a_edges[i]->l->v)); */
2647
2648
2649      /* RATES_Update_Cur_Bl(tree); */
2650      /* For(i,2*tree->n_otu-3) PhyML_Fprintf(fp,"%f\t",tree->a_edges[i]->l->v); */
2651
2652
2653      For(i,2*tree->n_otu-2) tree->rates->mean_r[i] = EXP(tree->rates->br_r[i]);
2654      For(i,2*tree->n_otu-1) tree->rates->mean_t[i] = tree->rates->nd_t[i];
2655     
2656      /* Time_To_Branch(tree); */
2657      /* char *s; */
2658      /* s = (char *)mCalloc(T_MAX_NAME,sizeof(char)); */
2659      /* strcpy(s,mcmc->io->in_align_file); */
2660      /* strcat(s,"_"); */
2661      /* strcat(s,mcmc->out_filename); */
2662      /* strcat(s,".ps"); */
2663      /* DR_Draw_Tree(s,tree); */
2664      /* Free(s);        */
2665
2666      /* FILE *fp; */
2667      /* int j; */
2668      /* t_node *d, *v1, *v2; */
2669      /* int n1, n2; */
2670      /* phydbl r1, r2; */
2671
2672      /* s = (char *)mCalloc(T_MAX_NAME,sizeof(char)); */
2673      /* strcpy(s,mcmc->io->in_align_file); */
2674      /* strcat(s,"_"); */
2675      /* strcat(s,mcmc->out_filename); */
2676      /* strcat(s,".corr"); */
2677      /* fp = fopen(s,"w"); */
2678
2679      /* n1 = n2 = 0; */
2680      /* r1 = r2 = 0.f; */
2681      /* For(i,2*tree->n_otu-3) */
2682      /*        { */
2683      /*          if(tree->a_nodes[i]->tax == NO) */
2684      /*            { */
2685      /*              d = tree->a_nodes[i]; */
2686      /*              v1 = v2 = NULL; */
2687      /*              For(j,3) */
2688      /*                { */
2689      /*                  if(d->v[j] != d->anc && d->b[j] != tree->e_root) */
2690      /*                    { */
2691      /*                      if(!v1) v1 = d->v[j]; */
2692      /*                      else    v2 = d->v[j]; */
2693      /*                    } */
2694      /*                } */
2695
2696      /*              For(j,3) */
2697      /*                { */
2698      /*                  if(v1->v[j] && v1->v[j] == d) */
2699      /*                    { */
2700      /*                      n1 = v1->bip_size[j]; */
2701      /*                      /\* r1 = RATES_Get_Mean_Rate_In_Subtree(v1,tree); *\/ */
2702      /*                      r1 = EXP(tree->rates->br_r[v1->num]); */
2703      /*                      break; */
2704      /*                    } */
2705      /*                } */
2706
2707
2708      /*              For(j,3) */
2709      /*                { */
2710      /*                  if(v2->v[j] && v2->v[j] == d) */
2711      /*                    { */
2712      /*                      n2 = v2->bip_size[j]; */
2713      /*                      /\* r2 = RATES_Get_Mean_Rate_In_Subtree(v2,tree); *\/ */
2714      /*                      r2 = EXP(tree->rates->br_r[v2->num]); */
2715      /*                      break; */
2716      /*                    } */
2717      /*                } */
2718
2719      /*              fprintf(fp,"\n%4d %4d %15f %15f",n1,n2,r1,r2); */
2720      /*            } */
2721      /*        } */
2722
2723      /* fclose(fp); */
2724     
2725      // TREES
2726      Time_To_Branch(tree);
2727      tree->bl_ndigits = 3;
2728      s_tree = Write_Tree(tree,NO);
2729      tree->bl_ndigits = 7;
2730      PhyML_Fprintf(mcmc->out_fp_trees,"TREE %8d [%f] = [&R] %s\n",mcmc->run,tree->c_lnL,s_tree);
2731      Free(s_tree);
2732    }
2733
2734  if(tree->mcmc->run == mcmc->chain_len) PhyML_Fprintf(mcmc->out_fp_trees,"END;\n");
2735
2736  fflush(NULL);   
2737
2738  Free(s);
2739 
2740}
2741
2742//////////////////////////////////////////////////////////////
2743//////////////////////////////////////////////////////////////
2744
2745
2746void MCMC_Print_Means(t_mcmc *mcmc, t_tree *tree)
2747{
2748
2749  if(!(mcmc->run%mcmc->sample_interval)) 
2750    {
2751      int i;
2752      char *s;
2753
2754      s = (char *)mCalloc(T_MAX_FILE,sizeof(char));
2755
2756      strcpy(s,tree->mcmc->out_filename);
2757      strcat(s,".means");
2758     
2759      fclose(mcmc->out_fp_means);
2760
2761      mcmc->out_fp_means = fopen(s,"w");
2762     
2763      PhyML_Fprintf(mcmc->out_fp_means,"#");
2764      for(i=tree->n_otu;i<2*tree->n_otu-1;i++) PhyML_Fprintf(mcmc->out_fp_means,"T%d\t",i);       
2765
2766      PhyML_Fprintf(mcmc->out_fp_means,"\n");     
2767
2768      for(i=tree->n_otu;i<2*tree->n_otu-1;i++) tree->rates->t_mean[i] *= (phydbl)(mcmc->run / mcmc->sample_interval);
2769
2770      for(i=tree->n_otu;i<2*tree->n_otu-1;i++)
2771        {
2772          tree->rates->t_mean[i] += tree->rates->nd_t[i];
2773          tree->rates->t_mean[i] /= (phydbl)(mcmc->run / mcmc->sample_interval + 1);
2774
2775/*        PhyML_Fprintf(tree->mcmc->out_fp_means,"%d\t",tree->mcmc->run / tree->mcmc->sample_interval);    */
2776          PhyML_Fprintf(tree->mcmc->out_fp_means,"%.1f\t",tree->rates->t_mean[i]);
2777        }
2778
2779      PhyML_Fprintf(tree->mcmc->out_fp_means,"\n");
2780      fflush(NULL);
2781
2782      Free(s);
2783    }
2784}
2785
2786//////////////////////////////////////////////////////////////
2787//////////////////////////////////////////////////////////////
2788
2789
2790void MCMC_Print_Last(t_mcmc *mcmc, t_tree *tree)
2791{
2792
2793  if(!(mcmc->run%mcmc->sample_interval)) 
2794    {
2795      int i;
2796      char *s;
2797
2798      s = (char *)mCalloc(T_MAX_FILE,sizeof(char));
2799
2800      strcpy(s,tree->mcmc->out_filename);
2801      strcat(s,".lasts");
2802     
2803      fclose(mcmc->out_fp_last);
2804
2805      mcmc->out_fp_last = fopen(s,"w");
2806
2807/*       rewind(mcmc->out_fp_last); */
2808
2809      PhyML_Fprintf(mcmc->out_fp_last,"#");
2810      PhyML_Fprintf(tree->mcmc->out_fp_last,"Time\t");
2811
2812      for(i=tree->n_otu;i<2*tree->n_otu-1;i++)
2813        PhyML_Fprintf(tree->mcmc->out_fp_last,"T%d\t",i);
2814
2815      PhyML_Fprintf(tree->mcmc->out_fp_last,"\n");
2816
2817      if(mcmc->run)
2818        {
2819          time(&(mcmc->t_cur));
2820          PhyML_Fprintf(tree->mcmc->out_fp_last,"%d\t",(int)(mcmc->t_cur-mcmc->t_beg));
2821/*        PhyML_Fprintf(tree->mcmc->out_fp_last,"%d\t",(int)(mcmc->t_beg)); */
2822        }
2823
2824      for(i=tree->n_otu;i<2*tree->n_otu-1;i++) PhyML_Fprintf(tree->mcmc->out_fp_last,"%.1f\t",tree->rates->nd_t[i]);
2825
2826      PhyML_Fprintf(tree->mcmc->out_fp_last,"\n");
2827      fflush(NULL);
2828
2829      Free(s);
2830  }
2831}
2832
2833
2834//////////////////////////////////////////////////////////////
2835//////////////////////////////////////////////////////////////
2836
2837void MCMC_Pause(t_mcmc *mcmc)
2838{
2839  char choice;
2840  char *s;
2841  int len;
2842
2843  s = (char *)mCalloc(100,sizeof(char));
2844
2845 
2846  if(!(mcmc->run%mcmc->chain_len) && (mcmc->is_burnin == NO))
2847    {
2848      PhyML_Printf("\n. Do you wish to stop the analysis [N/y] ");
2849      if(!scanf("%c",&choice)) Exit("\n");
2850      if(choice == '\n') choice = 'N';
2851      else getchar(); /* \n */
2852     
2853      Uppercase(&choice);
2854         
2855      switch(choice)
2856        {
2857        case 'N': 
2858          {         
2859            len = 1E+4;
2860            PhyML_Printf("\n. How many extra generations is required [default: 1E+4] ");
2861            Getstring_Stdin(s);
2862            if(s[0] == '\0') len = 1E+4; 
2863            else len = (int)atof(s); 
2864
2865            if(len < 0)
2866              {
2867                PhyML_Printf("\n. The value entered must be an integer greater than 0.\n");
2868                Exit("\n");
2869              }     
2870            mcmc->chain_len += len;
2871            break;
2872          }
2873             
2874        case 'Y': 
2875          {
2876            PhyML_Printf("\n. Ok. Done.\n");
2877            Exit("\n");
2878            break;
2879          }
2880         
2881        default: 
2882          {
2883            PhyML_Printf("\n. Please enter 'Y' or 'N'.\n");
2884            Exit("\n");
2885          }
2886        }
2887    }
2888
2889  Free(s);
2890
2891}
2892
2893
2894//////////////////////////////////////////////////////////////
2895//////////////////////////////////////////////////////////////
2896
2897
2898void MCMC_Terminate()
2899{
2900  char choice;
2901  PhyML_Printf("\n\n. Do you really want to terminate [Y/n]: ");
2902  if(!scanf("%c",&choice)) Exit("\n");
2903  if(choice == '\n') choice = 'Y';
2904  else getchar(); /* \n */     
2905  Uppercase(&choice);     
2906  if(choice == 'Y') raise(SIGTERM);
2907}
2908
2909
2910//////////////////////////////////////////////////////////////
2911//////////////////////////////////////////////////////////////
2912
2913void MCMC_Init_MCMC_Struct(char *filename, option *io, t_mcmc *mcmc)
2914{
2915  int pid;
2916
2917  mcmc->io               = io;
2918  mcmc->is               = NO;
2919  mcmc->use_data         = YES;
2920  mcmc->run              = 0;
2921  mcmc->sample_interval  = 1E+3;
2922  mcmc->chain_len        = 1E+6;
2923  mcmc->chain_len_burnin = 1E+5;
2924  mcmc->randomize        = YES;
2925  mcmc->norm_freq        = 1E+3;
2926  mcmc->max_tune         = 1.E+20;
2927  mcmc->min_tune         = 1.E-10;
2928  mcmc->print_every      = 2;
2929  mcmc->is_burnin        = NO;
2930  mcmc->nd_t_digits      = 4;
2931
2932  if(filename)
2933    {
2934      char *s;
2935
2936      s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
2937
2938      strcpy(mcmc->out_filename,filename);
2939      pid = getpid();
2940      sprintf(mcmc->out_filename+strlen(mcmc->out_filename),".%d",pid);
2941
2942      strcpy(s,mcmc->io->in_align_file);
2943      strcat(s,"_");
2944      strcat(s,mcmc->out_filename);
2945      strcat(s,".stats");
2946      mcmc->out_fp_stats = fopen(s,"w");
2947
2948      strcpy(s,mcmc->io->in_align_file);
2949      strcat(s,"_");
2950      strcat(s,mcmc->out_filename);
2951      strcat(s,".trees");
2952      mcmc->out_fp_trees = fopen(s,"w");
2953
2954      strcpy(s,mcmc->io->in_align_file);
2955      strcat(s,"_");
2956      strcat(s,mcmc->out_filename);
2957      strcat(s,".constree");
2958      mcmc->out_fp_constree = fopen(s,"w");
2959
2960/*       strcpy(s,tree->mcmc->out_filename); */
2961/*       strcat(s,".means"); */
2962/*       tree->mcmc->out_fp_means = fopen(s,"w"); */
2963
2964/*       strcpy(s,tree->mcmc->out_filename); */
2965/*       strcat(s,".lasts"); */
2966/*       tree->mcmc->out_fp_last  = fopen(s,"w"); */
2967
2968      Free(s);
2969    }
2970  else 
2971    {
2972      mcmc->out_fp_stats = stderr;
2973      mcmc->out_fp_trees = stderr;
2974      /* tree->mcmc->out_fp_means = stderr; */
2975      /* tree->mcmc->out_fp_last  = stderr; */
2976    }
2977}
2978
2979//////////////////////////////////////////////////////////////
2980//////////////////////////////////////////////////////////////
2981
2982
2983void MCMC_Copy_MCMC_Struct(t_mcmc *ori, t_mcmc *cpy, char *filename)
2984{
2985  int pid;
2986  int i;
2987 
2988  cpy->use_data           = ori->use_data        ;
2989  cpy->sample_interval    = ori->sample_interval ;
2990  cpy->chain_len          = ori->chain_len       ;
2991  cpy->randomize          = ori->randomize       ;
2992  cpy->norm_freq          = ori->norm_freq       ;
2993  cpy->n_moves            = ori->n_moves         ;
2994  cpy->max_tune           = ori->max_tune        ;
2995  cpy->min_tune           = ori->min_tune        ;
2996  cpy->print_every        = ori->print_every     ;
2997  cpy->is_burnin          = ori->is_burnin       ;
2998  cpy->is                 = ori->is              ;
2999  cpy->io                 = ori->io              ;
3000  cpy->in_fp_par          = ori->in_fp_par       ;
3001  cpy->nd_t_digits        = ori->nd_t_digits     ;
3002
3003  For(i,cpy->n_moves) 
3004    {
3005      cpy->old_param_val[i]      = ori->old_param_val[i];
3006      cpy->new_param_val[i]      = ori->new_param_val[i];
3007      cpy->start_ess[i]          = ori->start_ess[i];
3008      cpy->ess_run[i]            = ori->ess_run[i];
3009      cpy->ess[i]                = ori->ess[i];
3010      cpy->sum_val[i]            = ori->sum_val[i];
3011      cpy->sum_valsq[i]          = ori->sum_valsq[i];
3012      cpy->first_val[i]          = ori->first_val[i];
3013      cpy->sum_curval_nextval[i] = ori->sum_curval_nextval[i];
3014      cpy->move_weight[i]        = ori->move_weight[i];
3015      cpy->run_move[i]           = ori->run_move[i];
3016      cpy->acc_move[i]           = ori->acc_move[i];
3017      cpy->prev_run_move[i]      = ori->prev_run_move[i];
3018      cpy->prev_acc_move[i]      = ori->prev_acc_move[i];
3019      cpy->acc_rate[i]           = ori->acc_rate[i];
3020      cpy->tune_move[i]          = ori->tune_move[i];
3021      strcpy(cpy->move_name[i],ori->move_name[i]);
3022      cpy->adjust_tuning[i]      = ori->adjust_tuning[i];
3023    }
3024
3025
3026 
3027  if(filename) 
3028    {
3029      char *s;
3030
3031      s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
3032
3033      strcpy(cpy->out_filename,filename);
3034      pid = getpid();
3035      sprintf(cpy->out_filename+strlen(cpy->out_filename),".%d",pid);
3036
3037      strcpy(s,cpy->io->in_align_file);
3038      strcat(s,"_");
3039      strcat(s,cpy->out_filename);
3040      strcat(s,".stats");
3041      cpy->out_fp_stats = fopen(s,"w");
3042
3043      strcpy(s,cpy->io->in_align_file);
3044      strcat(s,"_");
3045      strcat(s,cpy->out_filename);
3046      strcat(s,".trees");
3047      cpy->out_fp_trees = fopen(s,"w");
3048
3049      strcpy(s,cpy->io->in_align_file);
3050      strcat(s,"_");
3051      strcat(s,cpy->out_filename);
3052      strcat(s,".constree");
3053      cpy->out_fp_constree = fopen(s,"w");
3054 
3055      Free(s);
3056    }
3057  else 
3058    {
3059      cpy->out_fp_stats = stderr;
3060      cpy->out_fp_trees = stderr;
3061    }
3062}
3063
3064
3065//////////////////////////////////////////////////////////////
3066//////////////////////////////////////////////////////////////
3067
3068
3069void MCMC_Close_MCMC(t_mcmc *mcmc)
3070{
3071  fclose(mcmc->out_fp_trees);
3072  fclose(mcmc->out_fp_stats);
3073  fclose(mcmc->out_fp_constree);
3074  /* fclose(mcmc->out_fp_means); */
3075  /* fclose(mcmc->out_fp_last); */
3076}
3077
3078//////////////////////////////////////////////////////////////
3079//////////////////////////////////////////////////////////////
3080
3081
3082void MCMC_Randomize_Kappa(t_tree *tree)
3083{
3084  tree->mod->kappa->v = Uni()*5.;
3085}
3086
3087//////////////////////////////////////////////////////////////
3088//////////////////////////////////////////////////////////////
3089
3090
3091void MCMC_Randomize_Rate_Across_Sites(t_tree *tree)
3092{
3093  if(tree->mod->ras->n_catg == 1) return;
3094 
3095  if(tree->mod->ras->free_mixt_rates == YES)
3096    {
3097      int i;
3098      For(i,tree->mod->ras->n_catg-1) tree->mod->ras->gamma_r_proba_unscaled->v[i] = Uni()*100.;
3099      tree->mod->ras->gamma_r_proba_unscaled->v[tree->mod->ras->n_catg-1] = 100.;
3100      For(i,tree->mod->ras->n_catg) tree->mod->ras->gamma_rr_unscaled->v[i] = (phydbl)i+1.; /* Do not randomize those as their ordering matter */
3101    }
3102  else
3103    {
3104      tree->mod->ras->alpha->v = Uni()*5.;
3105    }
3106}
3107
3108//////////////////////////////////////////////////////////////
3109//////////////////////////////////////////////////////////////
3110
3111
3112void MCMC_Randomize_Covarion_Rates(t_tree *tree)
3113{
3114  if(tree->mod->use_m4mod == NO) return;
3115 
3116  if(tree->mod->m4mod->n_h == 1) return;
3117 
3118  int i;
3119
3120  For(i,tree->mod->m4mod->n_h)
3121    {
3122      tree->mod->m4mod->multipl_unscaled[i] = (phydbl)i+1.;
3123      tree->mod->m4mod->h_fq_unscaled[i] = Uni()*(100.-0.01) + 0.01;
3124    }
3125}
3126
3127//////////////////////////////////////////////////////////////
3128//////////////////////////////////////////////////////////////
3129
3130
3131void MCMC_Randomize_Covarion_Switch(t_tree *tree)
3132{
3133  if(tree->mod->use_m4mod == NO) return;
3134
3135  tree->mod->m4mod->delta = Uni()*(10.-0.01)+0.01;
3136}
3137
3138//////////////////////////////////////////////////////////////
3139//////////////////////////////////////////////////////////////
3140
3141
3142void MCMC_Randomize_Branch_Lengths(t_tree *tree)
3143{
3144  int i;
3145
3146  if(tree->mod->log_l == NO)
3147    For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = Rexp(10.);
3148  else
3149    For(i,2*tree->n_otu-3) tree->a_edges[i]->l->v = -4* Uni();
3150}
3151
3152//////////////////////////////////////////////////////////////
3153//////////////////////////////////////////////////////////////
3154
3155
3156void MCMC_Randomize_Node_Rates(t_tree *tree)
3157{
3158
3159  int i,err;
3160  phydbl mean_r, var_r;
3161  phydbl min_r, max_r;
3162
3163  mean_r = 1.0;
3164  var_r  = 0.5;
3165  min_r  = tree->rates->min_rate;
3166  max_r  = tree->rates->max_rate;
3167
3168  For(i,2*tree->n_otu-2)
3169    if(tree->a_nodes[i] != tree->n_root)
3170      tree->rates->nd_r[i] = Rnorm_Trunc(mean_r,SQRT(var_r),min_r,max_r,&err);
3171}
3172
3173//////////////////////////////////////////////////////////////
3174//////////////////////////////////////////////////////////////
3175
3176
3177void MCMC_Randomize_Rates(t_tree *tree)
3178{
3179
3180  /* Should be called once t_node times have been determined */
3181
3182  int i;
3183
3184  For(i,2*tree->n_otu-2) tree->rates->br_r[i] = 1.0;
3185
3186/*   For(i,2*tree->n_otu-2) */
3187/*     { */
3188/*       u = Uni(); */
3189/*       u = u * (r_max-r_min) + r_min; */
3190/*       tree->rates->br_r[i] = u; */
3191
3192/*       if(tree->rates->br_r[i] < tree->rates->min_rate) tree->rates->br_r[i] = tree->rates->min_rate;  */
3193/*       if(tree->rates->br_r[i] > tree->rates->max_rate) tree->rates->br_r[i] = tree->rates->max_rate;  */
3194/*     } */
3195
3196  MCMC_Randomize_Rates_Pre(tree->n_root,tree->n_root->v[2],tree);
3197  MCMC_Randomize_Rates_Pre(tree->n_root,tree->n_root->v[1],tree);
3198}
3199
3200//////////////////////////////////////////////////////////////
3201//////////////////////////////////////////////////////////////
3202
3203
3204void MCMC_Randomize_Rates_Pre(t_node *a, t_node *d, t_tree *tree)
3205{
3206  int i;
3207  phydbl mean_r, var_r;
3208  phydbl min_r, max_r;
3209  int err;
3210
3211/*   mean_r = tree->rates->br_r[a->num]; */
3212/*   var_r  = tree->rates->nu * (tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num]); */
3213  mean_r = 1.0;
3214  var_r  = 0.5;
3215  min_r  = tree->rates->min_rate;
3216  max_r  = tree->rates->max_rate;
3217 
3218  tree->rates->br_r[d->num] = Rnorm_Trunc(mean_r,SQRT(var_r),min_r,max_r,&err);
3219
3220  if(d->tax) return;
3221  else
3222    {
3223      For(i,3)
3224        if(d->v[i] != a && d->b[i] != tree->e_root)
3225          MCMC_Randomize_Rates_Pre(d,d->v[i],tree);
3226    }
3227}
3228
3229//////////////////////////////////////////////////////////////
3230//////////////////////////////////////////////////////////////
3231
3232void MCMC_Randomize_Birth(t_tree *tree)
3233{
3234  phydbl min_b,max_b;
3235  phydbl u;
3236
3237  min_b = tree->rates->birth_rate_min;
3238  max_b = MIN(0.5,tree->rates->birth_rate_max);
3239 
3240  u = Uni();
3241  tree->rates->birth_rate = (max_b - min_b) * u + min_b;
3242}
3243
3244//////////////////////////////////////////////////////////////
3245//////////////////////////////////////////////////////////////
3246
3247void MCMC_Randomize_Nu(t_tree *tree)
3248{
3249  phydbl min_nu,max_nu;
3250  phydbl u;
3251
3252  /* It is preferable to start with small values of nu
3253     as if is difficult for the MCMC sampler to sample
3254     equal rates on edge (i.e., molecular clock) since
3255     such combination of rate lies on the boundary of
3256     the space of all edge rate combination. We give
3257     here a bit of help to the sampler by considering
3258     starting points close to the molecular clock
3259     constraint.
3260  */
3261  min_nu = tree->rates->min_nu;
3262  max_nu = tree->rates->max_nu/10.;
3263 
3264  u = Uni();
3265  tree->rates->nu = (max_nu - min_nu) * u + min_nu;
3266}
3267
3268//////////////////////////////////////////////////////////////
3269//////////////////////////////////////////////////////////////
3270
3271
3272void MCMC_Randomize_Clock_Rate(t_tree *tree)
3273{
3274  phydbl u;
3275  u = Uni();
3276  tree->rates->clock_r = u * (tree->rates->max_clock - tree->rates->min_clock) + tree->rates->min_clock;
3277}
3278
3279//////////////////////////////////////////////////////////////
3280//////////////////////////////////////////////////////////////
3281
3282
3283void MCMC_Randomize_Alpha(t_tree *tree)
3284{
3285  phydbl u;
3286
3287  u = Uni();
3288  tree->rates->alpha = u*6.0+1.0;
3289}
3290
3291//////////////////////////////////////////////////////////////
3292//////////////////////////////////////////////////////////////
3293
3294
3295void MCMC_Randomize_Node_Times(t_tree *tree)
3296{
3297  phydbl t_sup, t_inf;
3298  phydbl u;
3299  int iter;
3300  int i;
3301  phydbl dt,min_dt;
3302  int min_node;
3303
3304  t_inf = tree->rates->t_prior_min[tree->n_root->num];
3305  t_sup = tree->rates->t_prior_max[tree->n_root->num];
3306
3307  u = Uni();
3308  u *= (t_sup - t_inf);
3309  u += t_inf;
3310 
3311  tree->rates->nd_t[tree->n_root->num] = u;
3312
3313  MCMC_Randomize_Node_Times_Top_Down(tree->n_root,tree->n_root->v[2],tree);
3314  MCMC_Randomize_Node_Times_Top_Down(tree->n_root,tree->n_root->v[1],tree);
3315
3316  min_node = -1;
3317  iter = 0;
3318  do
3319    {
3320      min_dt = MDBL_MAX;
3321      For(i,2*tree->n_otu-2) 
3322        {
3323          dt = tree->rates->nd_t[i] - tree->rates->nd_t[tree->a_nodes[i]->anc->num];
3324          if(dt < min_dt)
3325            {
3326              min_dt = dt;
3327              min_node = i;
3328            }
3329        }
3330
3331      if(min_dt > 0.01 * FABS(tree->rates->nd_t[tree->n_root->num])/(phydbl)(tree->n_otu-1)) break;
3332
3333      RATES_Record_Times(tree);
3334      For(i,2*tree->n_otu-1) 
3335        {
3336          if(tree->a_nodes[i]->tax == NO) 
3337            tree->rates->nd_t[i] -= 0.1*FABS(tree->rates->nd_t[tree->n_root->num])/(phydbl)(tree->n_otu-1);
3338
3339          if(tree->rates->nd_t[i] < tree->rates->t_prior_min[i] || 
3340             tree->rates->nd_t[i] > tree->rates->t_prior_max[i])
3341            {
3342              RATES_Reset_Times(tree);
3343              break;
3344            }
3345        }
3346
3347      MCMC_Randomize_Node_Times_Bottom_Up(tree->n_root,tree->n_root->v[2],tree);
3348      MCMC_Randomize_Node_Times_Bottom_Up(tree->n_root,tree->n_root->v[1],tree);
3349     
3350      iter++;
3351    }
3352  while(iter < 1000);
3353
3354  if(iter == 1000)
3355    {     
3356      PhyML_Printf("\n== min_dt = %f",min_dt);
3357      PhyML_Printf("\n== min->t=%f min->anc->t=%f",tree->rates->nd_t[min_node],tree->rates->nd_t[tree->a_nodes[min_node]->anc->num]);
3358      PhyML_Printf("\n== d up=%f down=%f",tree->rates->t_prior_min[min_node],tree->rates->t_prior_max[min_node]);
3359      PhyML_Printf("\n== a up=%f down=%f",tree->rates->t_prior_min[tree->a_nodes[min_node]->anc->num],tree->rates->t_prior_max[tree->a_nodes[min_node]->anc->num]);
3360      PhyML_Printf("\n== up=%f down=%f",tree->rates->t_prior_min[min_node],tree->rates->t_floor[tree->a_nodes[min_node]->anc->num]);
3361      PhyML_Printf("\n== min_node = %d",min_node);
3362      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3363      Warn_And_Exit("");
3364    }
3365
3366
3367/*   PhyML_Printf("\n. Needed %d iterations to randomize node heights.",iter); */
3368/*   TIMES_Print_Node_Times(tree->n_root,tree->n_root->v[2],tree); */
3369/*   TIMES_Print_Node_Times(tree->n_root,tree->n_root->v[1],tree); */
3370
3371  if(RATES_Check_Node_Times(tree))
3372    {
3373      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
3374      Warn_And_Exit("");
3375    }
3376}
3377
3378//////////////////////////////////////////////////////////////
3379//////////////////////////////////////////////////////////////
3380
3381
3382void MCMC_Randomize_Node_Times_Bottom_Up(t_node *a, t_node *d, t_tree *tree)
3383{
3384  if(d->tax) return;
3385  else
3386    {
3387      int i;     
3388      phydbl u;
3389      phydbl t_inf, t_sup;
3390      t_node *v1, *v2;
3391
3392
3393      For(i,3)
3394        {
3395          if((d->v[i] != a) && (d->b[i] != tree->e_root))
3396            {
3397              MCMC_Randomize_Node_Times_Bottom_Up(d,d->v[i],tree);
3398            }
3399        }
3400
3401      v1 = v2 = NULL;
3402      For(i,3)
3403        {
3404          if(d->v[i] != a && d->b[i] != tree->e_root)
3405            {
3406              if(!v1) v1 = d->v[i];
3407              else    v2 = d->v[i];
3408            }
3409        }
3410
3411      t_sup = MIN(tree->rates->nd_t[v1->num],tree->rates->nd_t[v2->num]);
3412      t_inf = tree->rates->nd_t[a->num];
3413
3414      u = Uni();
3415      u *= (t_sup - t_inf);
3416      u += t_inf;
3417     
3418 
3419      if(u > tree->rates->t_prior_min[d->num] && u < tree->rates->t_prior_max[d->num])
3420        tree->rates->nd_t[d->num] = u;
3421    } 
3422}
3423
3424//////////////////////////////////////////////////////////////
3425//////////////////////////////////////////////////////////////
3426
3427
3428void MCMC_Randomize_Node_Times_Top_Down(t_node *a, t_node *d, t_tree *tree)
3429{
3430  if(d->tax) return;
3431  else
3432    {
3433      int i;     
3434      phydbl u;
3435      phydbl t_inf, t_sup;
3436
3437      t_inf = MAX(tree->rates->nd_t[a->num],tree->rates->t_prior_min[d->num]);
3438      t_sup = tree->rates->t_prior_max[d->num];
3439
3440      u = Uni();
3441      u *= (t_sup - t_inf);
3442      u += t_inf;
3443     
3444      tree->rates->nd_t[d->num] = u;
3445
3446      For(i,3)
3447        {
3448          if((d->v[i] != a) && (d->b[i] != tree->e_root))
3449            {
3450              MCMC_Randomize_Node_Times_Top_Down(d,d->v[i],tree);
3451            }
3452        }
3453    }
3454}
3455
3456//////////////////////////////////////////////////////////////
3457//////////////////////////////////////////////////////////////
3458
3459
3460void MCMC_Get_Acc_Rates(t_mcmc *mcmc)
3461{
3462  int i;
3463  phydbl eps;
3464  int lag;
3465
3466
3467  if(mcmc->run < (int)(0.01*mcmc->chain_len)) lag = 100;
3468  else lag = 100;
3469
3470  eps = 1.E-6;
3471
3472  For(i,mcmc->n_moves)
3473    {
3474      if(mcmc->run_move[i] - mcmc->prev_run_move[i] > lag)
3475        {
3476          mcmc->acc_rate[i] = 
3477            (phydbl)(mcmc->acc_move[i] - mcmc->prev_acc_move[i] + eps) / 
3478            (phydbl)(mcmc->run_move[i] - mcmc->prev_run_move[i] + eps) ;
3479
3480
3481          mcmc->prev_run_move[i] = mcmc->run_move[i];
3482          mcmc->prev_acc_move[i] = mcmc->acc_move[i];
3483         
3484          MCMC_Adjust_Tuning_Parameter(i,mcmc);
3485        }
3486    }
3487}
3488
3489//////////////////////////////////////////////////////////////
3490//////////////////////////////////////////////////////////////
3491
3492
3493void MCMC_Adjust_Tuning_Parameter(int move, t_mcmc *mcmc)
3494{
3495  if(mcmc->adjust_tuning[move])
3496    {
3497      phydbl scale;
3498      phydbl rate;
3499      phydbl rate_inf,rate_sup;
3500     
3501      if(mcmc->run < (int)(0.01*mcmc->chain_len)) scale = 1.5;
3502      else scale = 1.2;
3503
3504      if(!strcmp(mcmc->move_name[move],"tree_height"))
3505        {
3506          rate_inf = 0.1;
3507          rate_sup = 0.1;
3508          /* rate_inf = 0.3; */
3509          /* rate_sup = 0.3; */
3510          /* rate_inf = 0.7; */
3511          /* rate_sup = 0.7; */
3512        }
3513      else if(!strcmp(mcmc->move_name[move],"subtree_height"))
3514        {
3515          rate_inf = 0.2;
3516          rate_sup = 0.2;
3517        }
3518      else if(!strcmp(mcmc->move_name[move],"updown_t_cr"))
3519        {
3520          rate_inf = 0.1;
3521          rate_sup = 0.1;
3522        }
3523      else if(!strcmp(mcmc->move_name[move],"clock"))
3524        {
3525          rate_inf = 0.1;
3526          rate_sup = 0.1;
3527        }
3528      /* if(!strcmp(mcmc->move_name[move],"tree_rates")) */
3529      /*        { */
3530      /*          rate_inf = 0.05; */
3531      /*          rate_sup = 0.05; */
3532      /*        } */
3533      else
3534        {
3535          rate_inf = 0.234; // Gareth Robert's magic number !
3536          rate_sup = 0.234;
3537        }
3538
3539      /* PhyML_Printf("\n. %s acc=%d run=%d tune=%f", */
3540      /*                   mcmc->move_name[move], */
3541      /*                   mcmc->acc_move[move], */
3542      /*                   mcmc->run_move[move], */
3543      /*                   mcmc->tune_move[move]); */
3544
3545      rate = mcmc->acc_rate[move];
3546     
3547      if(rate < rate_inf)     
3548        {
3549          mcmc->tune_move[move] /= scale;
3550        }
3551
3552      else if(rate > rate_sup) 
3553        {
3554          mcmc->tune_move[move] *= scale;
3555        }
3556
3557      if(mcmc->tune_move[move] > mcmc->max_tune) mcmc->tune_move[move] = mcmc->max_tune;
3558      if(mcmc->tune_move[move] < mcmc->min_tune) mcmc->tune_move[move] = mcmc->min_tune;
3559         
3560    }
3561}
3562
3563//////////////////////////////////////////////////////////////
3564//////////////////////////////////////////////////////////////
3565
3566
3567void MCMC_One_Length(t_edge *b, t_tree *tree)
3568{
3569  phydbl u;
3570  phydbl new_lnL_data, cur_lnL_data;
3571  phydbl ratio, alpha;
3572  phydbl new_l, cur_l;
3573  phydbl K,mult;
3574
3575
3576  cur_l        = b->l->v;
3577  cur_lnL_data = tree->c_lnL;
3578  K            = 0.1;
3579 
3580  u = Uni();
3581  mult = EXP(K*(u-0.5));
3582  /* mult = u*(K-1./K)+1./K; */
3583  new_l = cur_l * mult;
3584
3585  if(new_l < tree->mod->l_min || new_l > tree->mod->l_max) return;
3586
3587  b->l->v = new_l;
3588  new_lnL_data = Lk(b,tree);
3589/*   tree->both_sides = NO; */
3590/*   new_lnL_data = Lk(tree); */
3591
3592
3593  ratio =
3594    (new_lnL_data - cur_lnL_data) +
3595    (LOG(mult));
3596
3597
3598  ratio = EXP(ratio);
3599  alpha = MIN(1.,ratio);
3600 
3601  u = Uni();
3602 
3603  if(u > alpha) /* Reject */
3604    {
3605      b->l->v = cur_l;
3606      Update_PMat_At_Given_Edge(b,tree);
3607      tree->c_lnL = cur_lnL_data;
3608    }
3609
3610}
3611
3612//////////////////////////////////////////////////////////////
3613//////////////////////////////////////////////////////////////
3614
3615
3616void MCMC_Scale_Br_Lens(t_tree *tree)
3617{
3618  phydbl u;
3619  phydbl new_lnL_data, cur_lnL_data;
3620  phydbl ratio, alpha;
3621  phydbl K,mult;
3622  int i;
3623
3624  Record_Br_Len(tree);
3625
3626  cur_lnL_data = tree->c_lnL;
3627  K            = 1.2;
3628 
3629  u = Uni();
3630  mult = u*(K-1./K)+1./K;
3631
3632  For(i,2*tree->n_otu-3) 
3633    {
3634      tree->a_edges[i]->l->v *= mult;
3635      if(tree->a_edges[i]->l->v < tree->mod->l_min || 
3636         tree->a_edges[i]->l->v > tree->mod->l_max) return;
3637    }
3638
3639  Set_Both_Sides(NO,tree);
3640  new_lnL_data = Lk(NULL,tree);
3641
3642  ratio =
3643    (new_lnL_data - cur_lnL_data) +
3644    (2*tree->n_otu-5) * (LOG(mult));
3645
3646  ratio = EXP(ratio);
3647  alpha = MIN(1.,ratio);
3648 
3649  u = Uni();
3650 
3651  if(u > alpha) /* Reject */
3652    {
3653      Restore_Br_Len(tree);
3654      tree->c_lnL = cur_lnL_data;
3655    }
3656}
3657
3658//////////////////////////////////////////////////////////////
3659//////////////////////////////////////////////////////////////
3660
3661
3662void MCMC_Br_Lens(t_tree *tree)
3663{
3664  MCMC_Br_Lens_Pre(tree->a_nodes[0],
3665                   tree->a_nodes[0]->v[0],
3666                   tree->a_nodes[0]->b[0],tree);
3667
3668  /* int i; */
3669  /* For(i,2*tree->n_otu-3) */
3670  /*   { */
3671  /*     MCMC_One_Length(tree->a_edges[Rand_Int(0,2*tree->n_otu-4)],acc,run,tree); */
3672  /*   } */
3673}
3674
3675//////////////////////////////////////////////////////////////
3676//////////////////////////////////////////////////////////////
3677
3678
3679void MCMC_Br_Lens_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
3680{
3681  int i;
3682
3683
3684  if(a == tree->n_root || d == tree->n_root)
3685    {
3686      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
3687      Exit("\n");
3688    }
3689
3690  MCMC_One_Length(b,tree);
3691  if(d->tax) return;
3692  else 
3693    {
3694      For(i,3) 
3695        if(d->v[i] != a)
3696          {
3697            Update_P_Lk(tree,d->b[i],d);
3698            MCMC_Br_Lens_Pre(d,d->v[i],d->b[i],tree);
3699          }
3700      Update_P_Lk(tree,b,d);
3701    } 
3702}
3703
3704//////////////////////////////////////////////////////////////
3705//////////////////////////////////////////////////////////////
3706
3707
3708void MCMC_Kappa(t_tree *tree)
3709{
3710  int change,i;
3711
3712  change = NO;
3713
3714  Switch_Eigen(YES,tree->mod);
3715               
3716  if(tree->io->lk_approx == NORMAL)
3717    {
3718      tree->io->lk_approx = EXACT;
3719      if(tree->mcmc->use_data == YES) Lk(NULL,tree);
3720      change = YES;
3721    }
3722 
3723  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = NO;
3724  MCMC_Single_Param_Generic(&(tree->mod->kappa->v),0.,15.,tree->mcmc->num_move_kappa,
3725                            NULL,&(tree->c_lnL),
3726                            NULL,Wrap_Lk,tree->mcmc->move_type[tree->mcmc->num_move_kappa],NO,NULL,tree,NULL);
3727         
3728  if(change == YES)
3729    {
3730      tree->io->lk_approx = NORMAL;
3731      if(tree->mcmc->use_data == YES) Lk(NULL,tree);
3732    }
3733 
3734  Switch_Eigen(NO,tree->mod);
3735}
3736
3737//////////////////////////////////////////////////////////////
3738//////////////////////////////////////////////////////////////
3739
3740
3741void MCMC_Rate_Across_Sites(t_tree *tree)
3742{
3743  if(tree->mod->ras->n_catg == 1) return;
3744
3745  if(tree->mod->ras->free_mixt_rates == YES)
3746    {
3747      MCMC_Free_Mixt_Rate(tree);
3748    }
3749  else
3750    {
3751      MCMC_Alpha(tree);
3752    }
3753}
3754
3755//////////////////////////////////////////////////////////////
3756//////////////////////////////////////////////////////////////
3757
3758
3759void MCMC_Alpha(t_tree *tree)
3760{
3761  int i;
3762 
3763  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = NO;
3764  MCMC_Single_Param_Generic(&(tree->mod->ras->alpha->v),0.,100.,tree->mcmc->num_move_ras,
3765                            NULL,&(tree->c_lnL),
3766                            NULL,Wrap_Lk,tree->mcmc->move_type[tree->mcmc->num_move_ras],NO,NULL,tree,NULL);
3767}
3768
3769//////////////////////////////////////////////////////////////
3770//////////////////////////////////////////////////////////////
3771
3772
3773void MCMC_Free_Mixt_Rate(t_tree *tree)
3774{
3775  phydbl num,denom;
3776  phydbl Jnow,Jthen;
3777  phydbl *z,*y;
3778  phydbl y_cur,z_cur;
3779  phydbl low_bound,up_bound;
3780  int c,i;
3781  int c2updt; 
3782  phydbl u;
3783  phydbl hr;
3784  int n_moves;
3785  phydbl cur_lnL_data, new_lnL_data;
3786  phydbl ratio,alpha;
3787 
3788  cur_lnL_data = tree->c_lnL;
3789  new_lnL_data = tree->c_lnL;
3790
3791  c = tree->mod->ras->n_catg;
3792
3793  z = tree->mod->ras->gamma_rr_unscaled->v;
3794  y = tree->mod->ras->gamma_r_proba_unscaled->v;
3795
3796  num = z[c-1]*(y[c-1]-y[c-2]);
3797  denom = z[0]*y[0];
3798  for(i=1;i<c;i++) denom += z[i]*(y[i]-y[i-1]);
3799  denom = POW(denom,c);
3800  Jthen = num/denom;
3801
3802  n_moves = 0;
3803  do
3804    {
3805      n_moves++;
3806
3807      // Update frequencies
3808
3809      // Choose the class freq to update at random.
3810      c2updt = Rand_Int(0,c-2); 
3811
3812      // Proposal is uniform. Determine upper and lower bounds.
3813      u = Uni();
3814      low_bound = (c2updt==0)?(.0):(y[c2updt-1]);
3815      up_bound = (c2updt==c-1)?(100.0):(y[c2updt+1]);
3816      y_cur = y[c2updt];
3817      y[c2updt] = low_bound + u*(up_bound - low_bound);
3818     
3819      // Calculate the Jacobian for the change of variable from unscaled
3820      // frequencies to the frequencies themselves.
3821      num = z[c-1]*(y[c-1]-y[c-2]);
3822      denom = z[0]*y[0];
3823      for(i=1;i<c;i++) denom += z[i]*(y[i]-y[i-1]);
3824      denom = POW(denom,c);
3825      Jnow = num/denom;
3826 
3827      hr = Jnow/Jthen;
3828
3829      new_lnL_data = Lk(NULL,tree);
3830
3831      // Metropolis-Hastings step
3832      ratio = 0.;
3833      if(tree->mcmc->use_data == YES) ratio += (new_lnL_data - cur_lnL_data);
3834      ratio += LOG(hr);
3835      ratio = EXP(ratio);
3836      alpha = MIN(1.,ratio);
3837     
3838      /* printf("\n. class=%d new_val=%f cur_val=%f ratio=%f hr=%f y=%f denom=%f",c2updt,y[c2updt],y_cur,ratio,Jthen/Jnow,y[c-1],denom); */
3839
3840      u = Uni();
3841      if(u > alpha) // Reject
3842        {
3843          y[c2updt] = y_cur;
3844          tree->c_lnL = cur_lnL_data;
3845        }
3846      else // Accept
3847        {
3848          cur_lnL_data = new_lnL_data;
3849          // Update the Jacobian
3850          num = z[c-1]*(y[c-1]-y[c-2]);
3851          denom = z[0]*y[0];
3852          for(i=1;i<c;i++) denom += z[i]*(y[i]-y[i-1]);
3853          denom = POW(denom,c);
3854          Jthen = num/denom;
3855        }
3856     
3857
3858
3859
3860      // Update rates
3861
3862      // Choose the class freq to update at random.
3863      c2updt = Rand_Int(0,tree->mod->ras->n_catg-1);
3864
3865      // Proposal move.
3866      u = Uni();
3867
3868      /* K = tree->mcmc->tune_move[tree->mcmc->num_move_ras+c2updt+c]; */
3869      /* z_cur = z[c2updt]; */
3870      /* mult = EXP(K*(u-0.5)); */
3871      /* z[c2updt] *= mult; */
3872
3873      u = Uni();
3874      low_bound = (c2updt==0)?(.0):(z[c2updt-1]);
3875      up_bound = (c2updt==c-1)?(100.0):(z[c2updt+1]);
3876      z_cur = z[c2updt];
3877      z[c2updt] = low_bound + u*(up_bound - low_bound);
3878
3879     
3880      // Calculate the Jacobian for the change of variable from unscaled
3881      // frequencies to the frequencies themselves.
3882      num = z[c-1]*(y[c-1]-y[c-2]);
3883      denom = z[0]*y[0];
3884      for(i=1;i<c;i++) denom += z[i]*(y[i]-y[i-1]);
3885      denom = POW(denom,c);
3886      Jnow = num/denom;
3887
3888      hr = Jnow/Jthen;
3889     
3890      new_lnL_data = Lk(NULL,tree);
3891
3892      // Metropolis-Hastings step
3893      ratio = 0.;
3894      if(tree->mcmc->use_data == YES) ratio += (new_lnL_data - cur_lnL_data);
3895      ratio += LOG(hr);
3896      /* ratio += LOG(mult); */
3897
3898      ratio = EXP(ratio);
3899      alpha = MIN(1.,ratio);
3900
3901      u = Uni();
3902      if(u > alpha)
3903        {
3904          z[c2updt] = z_cur;
3905          tree->c_lnL = cur_lnL_data;
3906        }
3907      else
3908        {
3909          cur_lnL_data = new_lnL_data;
3910
3911          // Update the Jacobian
3912          num = z[c-1]*(y[c-1]-y[c-2]);
3913          denom = z[0]*y[0];
3914          for(i=1;i<c;i++) denom += z[i]*(y[i]-y[i-1]);
3915          denom = POW(denom,c);
3916          Jthen = num/denom;
3917        }
3918
3919    }while(n_moves != c);
3920}
3921
3922//////////////////////////////////////////////////////////////
3923//////////////////////////////////////////////////////////////
3924
3925
3926void MCMC_Covarion_Rates(t_tree *tree)
3927{
3928  int i, class;
3929  phydbl u;
3930  phydbl min,max;
3931
3932  if(tree->mod->use_m4mod == NO) return;
3933
3934  Switch_Eigen(YES,tree->mod);
3935
3936  For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
3937
3938  class = Rand_Int(0,tree->mod->m4mod->n_h-1);
3939
3940
3941  min = 0.01;
3942  max = 100.;
3943  u = Uni();
3944  if(u < .5)
3945    {
3946      if(!class)
3947        {
3948          min = 0.01;
3949          max = tree->mod->m4mod->multipl_unscaled[1];
3950        }
3951      else if(class == tree->mod->m4mod->n_h-1)
3952        {
3953          min = tree->mod->m4mod->multipl_unscaled[tree->mod->m4mod->n_h-2];
3954          max = +100.;
3955        }
3956      else
3957        {
3958          min = MIN(tree->mod->m4mod->multipl_unscaled[class-1],tree->mod->m4mod->multipl_unscaled[class+1]);
3959          max = MAX(tree->mod->m4mod->multipl_unscaled[class-1],tree->mod->m4mod->multipl_unscaled[class+1]);
3960        }
3961
3962      MCMC_Single_Param_Generic(&(tree->mod->m4mod->multipl_unscaled[class]),min,max,tree->mcmc->num_move_cov_rates+class+tree->mod->m4mod->n_h,
3963                                NULL,&(tree->c_lnL),
3964                                NULL,Wrap_Lk,tree->mcmc->move_type[tree->mcmc->num_move_cov_rates+class+tree->mod->m4mod->n_h],NO,NULL,tree,NULL);
3965    }
3966  else
3967    {
3968      MCMC_Single_Param_Generic(&(tree->mod->m4mod->h_fq_unscaled[class]),0.01,+100.,tree->mcmc->num_move_cov_rates+class,
3969                                NULL,&(tree->c_lnL),
3970                                NULL,Wrap_Lk,tree->mcmc->move_type[tree->mcmc->num_move_cov_rates+class],NO,NULL,tree,NULL);
3971    }
3972
3973  Switch_Eigen(NO,tree->mod);
3974
3975}
3976
3977//////////////////////////////////////////////////////////////
3978//////////////////////////////////////////////////////////////
3979
3980
3981void MCMC_Covarion_Switch(t_tree *tree)
3982{
3983  if(tree->mod->use_m4mod == NO) return;
3984
3985  Switch_Eigen(YES,tree->mod);
3986  MCMC_Single_Param_Generic(&(tree->mod->m4mod->delta),0.01,+100.,tree->mcmc->num_move_cov_switch,
3987                            NULL,&(tree->c_lnL),
3988                            NULL,Wrap_Lk,tree->mcmc->move_type[tree->mcmc->num_move_cov_switch],NO,NULL,tree,NULL); 
3989  Switch_Eigen(NO,tree->mod);
3990}
3991
3992//////////////////////////////////////////////////////////////
3993//////////////////////////////////////////////////////////////
3994
3995void MCMC_Birth_Rate(t_tree *tree)
3996{
3997  MCMC_Single_Param_Generic(&(tree->rates->birth_rate),
3998                            tree->rates->birth_rate_min,
3999                            tree->rates->birth_rate_max,
4000                            tree->mcmc->num_move_birth_rate,
4001                            &(tree->rates->c_lnL_times),NULL,
4002                            Wrap_Lk_Times,NULL,tree->mcmc->move_type[tree->mcmc->num_move_birth_rate],NO,NULL,tree,NULL); 
4003}
4004
4005//////////////////////////////////////////////////////////////
4006//////////////////////////////////////////////////////////////
4007
4008void MCMC_Nu(t_tree *tree)
4009{
4010  phydbl cur_nu,new_nu,cur_lnL_rate,new_lnL_rate;
4011  phydbl u,alpha,ratio;
4012  phydbl min_nu,max_nu;
4013  phydbl K;
4014  phydbl cur_lnL_data, new_lnL_data;
4015  int i;
4016
4017  Record_Br_Len(tree);
4018
4019  cur_nu        = -1.0;
4020  new_nu        = -1.0;
4021  ratio         =  0.0;
4022
4023  K = tree->mcmc->tune_move[tree->mcmc->num_move_nu];
4024
4025  cur_lnL_rate = tree->rates->c_lnL_rates;
4026  new_lnL_rate = tree->rates->c_lnL_rates;
4027
4028  cur_lnL_data = tree->c_lnL;
4029  new_lnL_data = tree->c_lnL;
4030 
4031  cur_nu       = tree->rates->nu;
4032
4033  min_nu = tree->rates->min_nu;
4034  max_nu = tree->rates->max_nu;
4035
4036  MCMC_Make_Move(&cur_nu,&new_nu,min_nu,max_nu,&ratio,K,tree->mcmc->move_type[tree->mcmc->num_move_nu]);
4037 
4038  if(new_nu < max_nu && new_nu > min_nu) 
4039    {
4040      tree->rates->nu = new_nu;
4041     
4042      new_lnL_rate = RATES_Lk_Rates(tree);     
4043
4044      if((tree->rates->model == GUINDON) && (tree->mcmc->use_data))
4045        {
4046          For(i,2*tree->n_otu-2) tree->rates->br_do_updt[i] = YES;
4047          new_lnL_data = Lk(NULL,tree);
4048        }
4049     
4050      ratio +=
4051        (new_lnL_rate - cur_lnL_rate);
4052
4053      ratio +=
4054        (new_lnL_data - cur_lnL_data);
4055
4056
4057      /* !!!!!!!!!!!!!!!! */
4058      /* Modelling exp(nu) and making move on nu */
4059      /* ratio += (new_nu - cur_nu); */
4060         
4061      /* Exponential prior on nu */
4062      /* ratio += LOG(Dexp(new_nu,10.) / Dexp(cur_nu,10.)); */
4063
4064 
4065      ratio = EXP(ratio);
4066      alpha = MIN(1.,ratio);
4067     
4068      u = Uni();
4069      if(u > alpha) /* Reject */
4070        {
4071          tree->rates->nu    = cur_nu;
4072          tree->rates->c_lnL_rates = cur_lnL_rate;
4073          tree->c_lnL        = cur_lnL_data;
4074          Restore_Br_Len(tree);
4075        }
4076      else
4077        {
4078          tree->mcmc->acc_move[tree->mcmc->num_move_nu]++;
4079        }
4080    }
4081  tree->mcmc->run_move[tree->mcmc->num_move_nu]++;
4082}
4083
4084//////////////////////////////////////////////////////////////
4085//////////////////////////////////////////////////////////////
4086
4087
4088void MCMC_All_Rates(t_tree *tree)
4089{
4090  phydbl cur_lnL_data, new_lnL_data, cur_lnL_rate;
4091  phydbl u, ratio, alpha;
4092
4093  new_lnL_data = tree->c_lnL;
4094  cur_lnL_data = tree->c_lnL;
4095  cur_lnL_rate = tree->rates->c_lnL_rates;
4096  ratio        = 0.0;
4097 
4098  Record_Br_Len(tree);
4099  RATES_Record_Rates(tree);
4100 
4101  MCMC_Sim_Rate(tree->n_root,tree->n_root->v[2],tree);
4102  MCMC_Sim_Rate(tree->n_root,tree->n_root->v[1],tree);
4103
4104  new_lnL_data = Lk(NULL,tree);
4105 
4106  ratio += (new_lnL_data - cur_lnL_data);
4107  ratio = EXP(ratio);
4108
4109  alpha = MIN(1.,ratio);
4110  u = Uni();
4111
4112  if(u > alpha) /* Reject */
4113    {
4114      Restore_Br_Len(tree);
4115      RATES_Reset_Rates(tree);
4116      tree->rates->c_lnL_rates = cur_lnL_rate;
4117    }
4118  else
4119    {
4120      tree->rates->c_lnL_rates = RATES_Lk_Rates(tree);;
4121    }
4122}
4123
4124//////////////////////////////////////////////////////////////
4125//////////////////////////////////////////////////////////////
4126
4127
4128/* Only works when simulating from prior */
4129void MCMC_Sim_Rate(t_node *a, t_node *d, t_tree *tree)
4130{
4131  int err;
4132  phydbl mean,sd,br_r_a,dt_d;
4133
4134  br_r_a = tree->rates->br_r[a->num];
4135  dt_d   = tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num];
4136  sd     = SQRT(dt_d*tree->rates->nu);
4137 
4138  if(tree->rates->model_log_rates == YES)
4139    {
4140      mean = br_r_a - .5*sd*sd;
4141    }
4142  else
4143    {
4144      mean = br_r_a;
4145    }
4146     
4147  if(tree->rates->model == STRICTCLOCK) tree->rates->br_r[d->num] = 1.0;
4148  else
4149    tree->rates->br_r[d->num] = Rnorm_Trunc(mean,sd,tree->rates->min_rate,tree->rates->max_rate,&err);
4150
4151  if(d->tax) return;
4152  else
4153    {
4154      int i;
4155
4156      For(i,3)
4157        if(d->v[i] != a && d->b[i] != tree->e_root)
4158          MCMC_Sim_Rate(d,d->v[i],tree);
4159    }
4160}
4161
4162//////////////////////////////////////////////////////////////
4163//////////////////////////////////////////////////////////////
4164
4165
4166void MCMC_Complete_MCMC(t_mcmc *mcmc, t_tree *tree)
4167{
4168  int i;
4169  phydbl sum;
4170
4171  mcmc->io      = tree->io;
4172  mcmc->n_moves = 0;
4173
4174  mcmc->num_move_nd_r = mcmc->n_moves;
4175  mcmc->n_moves += 2*tree->n_otu-1;
4176
4177  mcmc->num_move_br_r = mcmc->n_moves;
4178  mcmc->n_moves += 2*tree->n_otu-2;
4179
4180  mcmc->num_move_nd_t = mcmc->n_moves;
4181  mcmc->n_moves += tree->n_otu-1;
4182
4183  mcmc->num_move_nu                    = mcmc->n_moves; mcmc->n_moves += 1;
4184  mcmc->num_move_clock_r               = mcmc->n_moves; mcmc->n_moves += 1;
4185  mcmc->num_move_tree_height           = mcmc->n_moves; mcmc->n_moves += 1;
4186  mcmc->num_move_subtree_height        = mcmc->n_moves; mcmc->n_moves += 1;
4187  mcmc->num_move_kappa                 = mcmc->n_moves; mcmc->n_moves += 1;
4188  mcmc->num_move_tree_rates            = mcmc->n_moves; mcmc->n_moves += 1;
4189  mcmc->num_move_subtree_rates         = mcmc->n_moves; mcmc->n_moves += 1;
4190  mcmc->num_move_updown_nu_cr          = mcmc->n_moves; mcmc->n_moves += 1;
4191  mcmc->num_move_ras                   = mcmc->n_moves; mcmc->n_moves += (tree->mod ? 2*tree->mod->ras->n_catg : 1);
4192  mcmc->num_move_updown_t_cr           = mcmc->n_moves; mcmc->n_moves += 1;
4193  mcmc->num_move_cov_rates             = mcmc->n_moves; mcmc->n_moves += (tree->mod ? 2*tree->mod->m4mod->n_h : 1);
4194  mcmc->num_move_cov_switch            = mcmc->n_moves; mcmc->n_moves += 1;
4195  mcmc->num_move_birth_rate            = mcmc->n_moves; mcmc->n_moves += 1;
4196  mcmc->num_move_updown_t_br           = mcmc->n_moves; mcmc->n_moves += 1;
4197  mcmc->num_move_jump_calibration      = mcmc->n_moves; mcmc->n_moves += 1;
4198  mcmc->num_move_geo_lambda            = mcmc->n_moves; mcmc->n_moves += 1;
4199  mcmc->num_move_geo_sigma             = mcmc->n_moves; mcmc->n_moves += 1;
4200  mcmc->num_move_geo_tau               = mcmc->n_moves; mcmc->n_moves += 1;
4201  mcmc->num_move_geo_updown_tau_lbda   = mcmc->n_moves; mcmc->n_moves += 1;
4202  mcmc->num_move_geo_updown_lbda_sigma = mcmc->n_moves; mcmc->n_moves += 1;
4203
4204  mcmc->run_move           = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4205  mcmc->acc_move           = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4206  mcmc->prev_run_move      = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4207  mcmc->prev_acc_move      = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4208  mcmc->acc_rate           = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4209  mcmc->move_weight        = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4210  mcmc->move_type          = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4211 
4212  /* TO DO: instead of n_moves here we should have something like n_param */
4213  mcmc->sum_val            = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4214  mcmc->first_val          = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4215  mcmc->sum_valsq          = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4216  mcmc->sum_curval_nextval = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4217  mcmc->ess                = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4218  mcmc->ess_run            = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4219  mcmc->start_ess          = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4220  mcmc->new_param_val      = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4221  mcmc->old_param_val      = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4222  mcmc->adjust_tuning      = (int *)mCalloc(mcmc->n_moves,sizeof(int));
4223  mcmc->tune_move          = (phydbl *)mCalloc(mcmc->n_moves,sizeof(phydbl));
4224
4225  mcmc->move_name = (char **)mCalloc(mcmc->n_moves,sizeof(char *));
4226  For(i,mcmc->n_moves) mcmc->move_name[i] = (char *)mCalloc(50,sizeof(char));
4227
4228  For(i,mcmc->n_moves) mcmc->adjust_tuning[i] = YES;
4229
4230  for(i=mcmc->num_move_br_r;i<mcmc->num_move_br_r+2*tree->n_otu-2;i++) strcpy(mcmc->move_name[i],"br_rate"); 
4231  for(i=mcmc->num_move_nd_r;i<mcmc->num_move_nd_r+2*tree->n_otu-1;i++) strcpy(mcmc->move_name[i],"nd_rate");
4232  for(i=mcmc->num_move_nd_t;i<mcmc->num_move_nd_t+tree->n_otu-1;i++)   strcpy(mcmc->move_name[i],"time");
4233  strcpy(mcmc->move_name[mcmc->num_move_nu],"nu");
4234  strcpy(mcmc->move_name[mcmc->num_move_clock_r],"clock");
4235  strcpy(mcmc->move_name[mcmc->num_move_tree_height],"tree_height");
4236  strcpy(mcmc->move_name[mcmc->num_move_subtree_height],"subtree_height");
4237  strcpy(mcmc->move_name[mcmc->num_move_kappa],"kappa");
4238  strcpy(mcmc->move_name[mcmc->num_move_tree_rates],"tree_rates");
4239  strcpy(mcmc->move_name[mcmc->num_move_subtree_rates],"subtree_rates");
4240  strcpy(mcmc->move_name[mcmc->num_move_updown_nu_cr],"updown_nu_cr");
4241  for(i=mcmc->num_move_ras;i<mcmc->num_move_ras+(tree->mod ? 2*tree->mod->ras->n_catg : 1);i++) 
4242    strcpy(mcmc->move_name[i],"ras"); 
4243  strcpy(mcmc->move_name[mcmc->num_move_updown_t_cr],"updown_t_cr");
4244  for(i=mcmc->num_move_cov_rates;i<mcmc->num_move_cov_rates+(tree->mod ? 2*tree->mod->m4mod->n_h : 1);i++) 
4245    strcpy(mcmc->move_name[i],"cov_rates"); 
4246  strcpy(mcmc->move_name[mcmc->num_move_cov_switch],"cov_switch");
4247  strcpy(mcmc->move_name[mcmc->num_move_birth_rate],"birth_rate");
4248  strcpy(mcmc->move_name[mcmc->num_move_updown_t_br],"updown_t_br");
4249  strcpy(mcmc->move_name[mcmc->num_move_jump_calibration],"jump_calibration");
4250  strcpy(mcmc->move_name[mcmc->num_move_geo_lambda],"geo_lambda");
4251  strcpy(mcmc->move_name[mcmc->num_move_geo_sigma],"geo_sigma");
4252  strcpy(mcmc->move_name[mcmc->num_move_geo_tau],"geo_tau");
4253  strcpy(mcmc->move_name[mcmc->num_move_geo_updown_tau_lbda],"geo_updown_tau_lbda");
4254  strcpy(mcmc->move_name[mcmc->num_move_geo_updown_lbda_sigma],"geo_updown_lbda_sigma");
4255
4256 
4257  if(tree->rates->model_log_rates == YES)
4258    for(i=mcmc->num_move_br_r;i<mcmc->num_move_br_r+2*tree->n_otu-2;i++) mcmc->move_type[i] = MCMC_MOVE_RANDWALK_NORMAL; 
4259  else
4260    for(i=mcmc->num_move_br_r;i<mcmc->num_move_br_r+2*tree->n_otu-2;i++) mcmc->move_type[i] = MCMC_MOVE_SCALE_THORNE; 
4261
4262  for(i=mcmc->num_move_nd_r;i<mcmc->num_move_nd_r+2*tree->n_otu-1;i++) mcmc->move_type[i] = MCMC_MOVE_SCALE_THORNE;
4263  for(i=mcmc->num_move_nd_t;i<mcmc->num_move_nd_t+tree->n_otu-1;i++)   mcmc->move_type[i] = MCMC_MOVE_RANDWALK_UNIFORM;
4264  /* for(i=mcmc->num_move_nd_t;i<mcmc->num_move_nd_t+tree->n_otu-1;i++)   mcmc->move_type[i] = MCMC_MOVE_SCALE_THORNE; */
4265  mcmc->move_type[mcmc->num_move_nu] = MCMC_MOVE_SCALE_THORNE;
4266  mcmc->move_type[mcmc->num_move_clock_r] = MCMC_MOVE_SCALE_THORNE;
4267  mcmc->move_type[mcmc->num_move_tree_height] = MCMC_MOVE_SCALE_THORNE;
4268  mcmc->move_type[mcmc->num_move_subtree_height] = MCMC_MOVE_SCALE_THORNE;
4269  mcmc->move_type[mcmc->num_move_kappa] = MCMC_MOVE_SCALE_THORNE;
4270  mcmc->move_type[mcmc->num_move_tree_rates] = MCMC_MOVE_SCALE_THORNE;
4271  mcmc->move_type[mcmc->num_move_subtree_rates] = MCMC_MOVE_SCALE_THORNE;
4272  mcmc->move_type[mcmc->num_move_updown_nu_cr] = MCMC_MOVE_RANDWALK_NORMAL;
4273  for(i=mcmc->num_move_ras;i<mcmc->num_move_ras+(tree->mod ? 2*tree->mod->ras->n_catg : 1);i++) mcmc->move_type[i] = MCMC_MOVE_RANDWALK_NORMAL; 
4274  mcmc->move_type[mcmc->num_move_updown_t_cr] = MCMC_MOVE_SCALE_THORNE;
4275  for(i=mcmc->num_move_cov_rates;i<mcmc->num_move_cov_rates+(tree-> mod ? 2*tree->mod->m4mod->n_h : 1);i++) mcmc->move_type[i] = MCMC_MOVE_SCALE_THORNE;
4276  mcmc->move_type[mcmc->num_move_cov_switch] = MCMC_MOVE_SCALE_THORNE;
4277  mcmc->move_type[mcmc->num_move_birth_rate] = MCMC_MOVE_SCALE_THORNE;
4278  mcmc->move_type[mcmc->num_move_updown_t_br] = MCMC_MOVE_SCALE_THORNE;
4279  mcmc->move_type[mcmc->num_move_jump_calibration] = MCMC_MOVE_SCALE_THORNE;
4280  mcmc->move_type[mcmc->num_move_geo_lambda] = MCMC_MOVE_SCALE_THORNE;
4281  mcmc->move_type[mcmc->num_move_geo_sigma] = MCMC_MOVE_SCALE_THORNE;
4282  mcmc->move_type[mcmc->num_move_geo_tau] = MCMC_MOVE_SCALE_THORNE;
4283  mcmc->move_type[mcmc->num_move_geo_updown_tau_lbda] = MCMC_MOVE_SCALE_THORNE;
4284  mcmc->move_type[mcmc->num_move_geo_updown_lbda_sigma] = MCMC_MOVE_SCALE_THORNE;
4285
4286  /* We start with small tuning parameter values in order to have inflated ESS
4287     for clock_r */
4288  For(i,mcmc->n_moves)
4289    {
4290      switch(mcmc->move_type[i])
4291        {
4292        case MCMC_MOVE_RANDWALK_NORMAL:
4293          {
4294            /* mcmc->tune_move[i] = 1.E-1; */
4295            mcmc->tune_move[i] = 100.;
4296            break;
4297          }
4298        case MCMC_MOVE_RANDWALK_UNIFORM:
4299          {
4300            /* mcmc->tune_move[i] = 2.0; */
4301            mcmc->tune_move[i] = 20.0;
4302            break;
4303          }
4304        case MCMC_MOVE_SCALE_GAMMA:
4305          {
4306            /* mcmc->tune_move[i] = 1.0; */
4307            mcmc->tune_move[i] = 10.0;
4308            break;
4309          }
4310        case MCMC_MOVE_SCALE_THORNE:
4311          {
4312            /* mcmc->tune_move[i] = 1.0; */
4313            mcmc->tune_move[i] = 10.0;
4314            break;
4315          }
4316        }
4317    }
4318 
4319  for(i=mcmc->num_move_br_r;i<mcmc->num_move_br_r+2*tree->n_otu-2;i++) mcmc->move_weight[i] = (phydbl)(1./(2.*tree->n_otu-2)); /* Rates */
4320  for(i=mcmc->num_move_nd_r;i<mcmc->num_move_nd_r+2*tree->n_otu-1;i++) mcmc->move_weight[i] = 0.0; /* Node rates */
4321  for(i=mcmc->num_move_nd_t;i<mcmc->num_move_nd_t+tree->n_otu-1;i++)   mcmc->move_weight[i] = (phydbl)(1./(tree->n_otu-1));  /* Times */
4322  mcmc->move_weight[mcmc->num_move_clock_r]          = 1.0;
4323  mcmc->move_weight[mcmc->num_move_tree_height]      = 2.0;
4324  mcmc->move_weight[mcmc->num_move_subtree_height]   = 0.0;
4325  mcmc->move_weight[mcmc->num_move_nu]               = 2.0;
4326  mcmc->move_weight[mcmc->num_move_kappa]            = 0.5;
4327  mcmc->move_weight[mcmc->num_move_tree_rates]       = 1.0;
4328  mcmc->move_weight[mcmc->num_move_subtree_rates]    = 0.5;
4329  mcmc->move_weight[mcmc->num_move_updown_nu_cr]     = 0.0;
4330  for(i=mcmc->num_move_ras;i<mcmc->num_move_ras+(tree->mod ? 2*tree->mod->ras->n_catg : 1);i++) mcmc->move_weight[i] = 0.5*(1./(tree->mod ? (phydbl)tree->mod->ras->n_catg : 1));
4331  mcmc->move_weight[mcmc->num_move_updown_t_cr]      = 0.0; /* Does not seem to work well (does not give uniform prior on root height
4332                                                              when sampling from prior) */
4333  for(i=mcmc->num_move_cov_rates;i<mcmc->num_move_cov_rates+(tree->mod ? 2*tree->mod->m4mod->n_h : 1);i++) mcmc->move_weight[i] = 0.5*(1./(tree->mod ? (phydbl)tree->mod->m4mod->n_h : 1));
4334  mcmc->move_weight[mcmc->num_move_cov_switch]            = 1.0;
4335  mcmc->move_weight[mcmc->num_move_birth_rate]            = 2.0;
4336  mcmc->move_weight[mcmc->num_move_updown_t_br]           = 1.0;
4337#if defined (SERGEII)
4338  mcmc->move_weight[mcmc->num_move_jump_calibration]      = 0.0;
4339#else
4340  mcmc->move_weight[mcmc->num_move_jump_calibration]      = 0.0;
4341#endif
4342  mcmc->move_weight[mcmc->num_move_geo_lambda]            = 1.0;
4343  mcmc->move_weight[mcmc->num_move_geo_sigma]             = 1.0;
4344  mcmc->move_weight[mcmc->num_move_geo_tau]               = 1.0;
4345  mcmc->move_weight[mcmc->num_move_geo_updown_tau_lbda]   = 1.0;
4346  mcmc->move_weight[mcmc->num_move_geo_updown_lbda_sigma] = 1.0;
4347
4348
4349
4350
4351/*   for(i=mcmc->num_move_br_r;i<mcmc->num_move_br_r+2*tree->n_otu-2;i++) mcmc->move_weight[i] = 0.0; /\* Rates *\/ */
4352/*   for(i=mcmc->num_move_nd_r;i<mcmc->num_move_nd_r+2*tree->n_otu-1;i++) mcmc->move_weight[i] = 0.0; /\* Node rates *\/ */
4353/*   for(i=mcmc->num_move_nd_t;i<mcmc->num_move_nd_t+tree->n_otu-1;i++)   mcmc->move_weight[i] = (phydbl)(1./(tree->n_otu-1));  /\* Times *\/ */
4354/*   mcmc->move_weight[mcmc->num_move_clock_r]          = 0.0; */
4355/*   mcmc->move_weight[mcmc->num_move_tree_height]      = 0.0; */
4356/*   mcmc->move_weight[mcmc->num_move_subtree_height]   = 0.0; */
4357/*   mcmc->move_weight[mcmc->num_move_tree_height]      = 0.0; */
4358/*   mcmc->move_weight[mcmc->num_move_subtree_height]   = 0.0; */
4359/*   mcmc->move_weight[mcmc->num_move_nu]               = 0.0; */
4360/*   mcmc->move_weight[mcmc->num_move_kappa]            = 0.0; */
4361/*   mcmc->move_weight[mcmc->num_move_tree_rates]       = 0.0; */
4362/*   mcmc->move_weight[mcmc->num_move_subtree_rates]    = 0.0; */
4363/*   mcmc->move_weight[mcmc->num_move_updown_nu_cr]     = 0.0; */
4364/*   for(i=mcmc->num_move_ras;i<mcmc->num_move_ras+(tree->mod ? 2*tree->mod->ras->n_catg : 1);i++) mcmc->move_weight[i] = 0.0; */
4365/*   mcmc->move_weight[mcmc->num_move_updown_t_cr]      = 0.0; /\* Does not seem to work well (does not give uniform prior on root height */
4366/*                                                            when sampling from prior) *\/ */
4367/*   for(i=mcmc->num_move_cov_rates;i<mcmc->num_move_cov_rates+(tree->mod ? 2*tree->mod->m4mod->n_h : 1);i++) mcmc->move_weight[i] = 0.5*(1./(tree->mod ? (phydbl)tree->mod->m4mod->n_h : 1)); */
4368/*   mcmc->move_weight[mcmc->num_move_cov_switch]       = 0.0; */
4369/*   mcmc->move_weight[mcmc->num_move_birth_rate]       = 0.0; */
4370/*   mcmc->move_weight[mcmc->num_move_updown_t_br]      = 0.0; */
4371/* #if defined (SERGEII) */
4372/*   mcmc->move_weight[mcmc->num_move_jump_calibration] = 0.0; */
4373/* #else */
4374/*   mcmc->move_weight[mcmc->num_move_jump_calibration] = 0.0; */
4375/* #endif */
4376/*   mcmc->move_weight[mcmc->num_move_geo_lambda]       = 0.0; */
4377/*   mcmc->move_weight[mcmc->num_move_geo_sigma]        = 0.0; */
4378/*   mcmc->move_weight[mcmc->num_move_geo_tau]          = 0.0; */
4379
4380
4381  sum = 0.0;
4382  For(i,mcmc->n_moves) sum += mcmc->move_weight[i];
4383  For(i,mcmc->n_moves) mcmc->move_weight[i] /= sum;
4384  for(i=1;i<mcmc->n_moves;i++) mcmc->move_weight[i] += mcmc->move_weight[i-1];
4385}
4386
4387//////////////////////////////////////////////////////////////
4388//////////////////////////////////////////////////////////////
4389
4390
4391void MCMC_Initialize_Param_Val(t_mcmc *mcmc, t_tree *tree)
4392{
4393  int i;
4394
4395  mcmc->old_param_val[mcmc->num_move_nu]          = tree->rates->nu;
4396  mcmc->old_param_val[mcmc->num_move_clock_r]     = tree->rates->clock_r;
4397  mcmc->old_param_val[mcmc->num_move_tree_height] = tree->rates->nd_t[tree->n_root->num];
4398  mcmc->old_param_val[mcmc->num_move_kappa]       = tree->mod->kappa->v;
4399
4400  For(i,2*tree->n_otu-2)
4401    mcmc->old_param_val[mcmc->num_move_br_r+i] = tree->rates->br_r[i];
4402 
4403  For(i,tree->n_otu-1)
4404    mcmc->old_param_val[mcmc->num_move_nd_t+i] = tree->rates->nd_t[tree->n_otu+i];
4405
4406  For(i,2*tree->n_otu-1)
4407    mcmc->old_param_val[mcmc->num_move_nd_r+i] = tree->rates->nd_r[i];
4408
4409  For(i,mcmc->n_moves) tree->mcmc->new_param_val[i] = tree->mcmc->old_param_val[i];
4410}
4411
4412//////////////////////////////////////////////////////////////
4413//////////////////////////////////////////////////////////////
4414
4415
4416void MCMC_Copy_To_New_Param_Val(t_mcmc *mcmc, t_tree *tree)
4417{
4418  int i;
4419
4420  mcmc->new_param_val[mcmc->num_move_nu]          = tree->rates->nu;
4421  mcmc->new_param_val[mcmc->num_move_clock_r]     = tree->rates->clock_r;
4422  mcmc->new_param_val[mcmc->num_move_tree_height] = tree->rates->nd_t[tree->n_root->num];
4423  mcmc->new_param_val[mcmc->num_move_kappa]       = tree->mod ? tree->mod->kappa->v : -1.;
4424  mcmc->new_param_val[mcmc->num_move_birth_rate]  = tree->rates->birth_rate;
4425
4426
4427  mcmc->new_param_val[mcmc->num_move_geo_tau]    = tree->geo ? tree->geo->tau   : -1.;
4428  mcmc->new_param_val[mcmc->num_move_geo_lambda] = tree->geo ? tree->geo->lbda  : -1.;
4429  mcmc->new_param_val[mcmc->num_move_geo_sigma]  = tree->geo ? tree->geo->sigma : -1.;
4430
4431
4432  For(i,2*tree->n_otu-2)
4433    mcmc->new_param_val[mcmc->num_move_br_r+i] = tree->rates->br_r[i];
4434 
4435  For(i,tree->n_otu-1)
4436    mcmc->new_param_val[mcmc->num_move_nd_t+i] = tree->rates->nd_t[tree->n_otu+i];
4437
4438  For(i,2*tree->n_otu-1)
4439    mcmc->new_param_val[mcmc->num_move_nd_r+i] = tree->rates->nd_r[i];
4440 
4441}
4442
4443//////////////////////////////////////////////////////////////
4444//////////////////////////////////////////////////////////////
4445
4446
4447void MCMC_Slice_One_Rate(t_node *a, t_node *d, int traversal, t_tree *tree)
4448{
4449  phydbl L,R; /* Left and Right limits of the slice */
4450  phydbl w; /* window width */
4451  phydbl u;
4452  phydbl x0,x1;
4453  phydbl logy;
4454  t_edge *b;
4455  int i;
4456
4457
4458
4459  b = NULL;
4460  if(a == tree->n_root) b = tree->e_root;
4461  else For(i,3) if(d->v[i] == a) { b = d->b[i]; break; }
4462 
4463  w = 0.05;
4464  /* w = 10.; */
4465
4466  x0 = tree->rates->br_r[d->num];
4467  logy = tree->c_lnL+tree->rates->c_lnL_rates - Rexp(1.);
4468 
4469  u = Uni();
4470
4471  L = x0 - w*u; 
4472  R = L + w;
4473 
4474  do
4475    {
4476      tree->rates->br_r[d->num] = L;
4477      tree->rates->br_do_updt[d->num] = YES;
4478      RATES_Update_Cur_Bl(tree);
4479      Lk(b,tree);
4480      RATES_Lk_Rates(tree);
4481      if(L < tree->rates->min_rate) { L = tree->rates->min_rate - w; break;}
4482      L = L - w;
4483    }
4484  while(tree->c_lnL + tree->rates->c_lnL_rates > logy);
4485  L = L + w;
4486
4487
4488  do
4489    {
4490      tree->rates->br_r[d->num] = R;
4491      tree->rates->br_do_updt[d->num] = YES;
4492      RATES_Update_Cur_Bl(tree);
4493      Lk(b,tree);
4494      RATES_Lk_Rates(tree);
4495      if(R > tree->rates->max_rate) { R = tree->rates->max_rate + w; break;}
4496      R = R + w;
4497    }
4498  while(tree->c_lnL + tree->rates->c_lnL_rates > logy);
4499  R = R - w;
4500
4501
4502  for(;;)
4503    {
4504      u = Uni();
4505      x1 = L + u*(R-L);
4506
4507      tree->rates->br_r[d->num] = x1;
4508      tree->rates->br_do_updt[d->num] = YES;
4509      RATES_Update_Cur_Bl(tree);
4510      Lk(b,tree);
4511      RATES_Lk_Rates(tree);
4512     
4513      if(tree->c_lnL + tree->rates->c_lnL_rates > logy) break;
4514     
4515      if(x1 < x0) L = x1;
4516      else        R = x1;
4517    }
4518
4519
4520  if(traversal == YES)
4521    {
4522      if(d->tax == YES) return;
4523      else
4524        {
4525          For(i,3)
4526            if(d->v[i] != a && d->b[i] != tree->e_root)
4527              {
4528                if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,d->b[i],d);
4529                /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) { tree->both_sides = YES; Lk(tree); } */
4530                MCMC_Slice_One_Rate(d,d->v[i],YES,tree);
4531              }
4532        }
4533      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,b,d);
4534      /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) { tree->both_sides = YES; Lk(tree); } */
4535    }
4536}
4537
4538//////////////////////////////////////////////////////////////
4539//////////////////////////////////////////////////////////////
4540
4541void MCMC_Make_Move(phydbl *cur, phydbl *new, phydbl inf, phydbl sup, phydbl *loghr, phydbl tune, int move_type)
4542{
4543  phydbl u;
4544
4545
4546     
4547  switch(move_type)
4548    {
4549
4550    case MCMC_MOVE_RANDWALK_NORMAL:
4551      {
4552        *new   = *cur + Rnorm(0.0,tune);
4553        /* Do not use reflection here */
4554        *loghr = 0.0;
4555
4556        break;
4557      }
4558
4559    case MCMC_MOVE_RANDWALK_UNIFORM:
4560      {
4561        u = Uni(); 
4562        /* *new   = u * (2.*tune) + (*cur) - tune; */
4563        /* *new   = Reflect(*new,inf,sup); */
4564        *new   = u*(sup-inf)+inf; 
4565        *loghr = 0.0;
4566        break;
4567      }
4568
4569    case MCMC_MOVE_SCALE_THORNE:
4570      {
4571        u = Uni();
4572        *new   = (*cur) * EXP(tune*(u-.5));
4573        *loghr = LOG((*new)/(*cur));
4574        break;
4575      }
4576
4577    case MCMC_MOVE_SCALE_GAMMA:
4578      {
4579        phydbl r;
4580        *new = (*cur) * Rgamma(1./tune,tune);
4581        r = (*new)/(*cur);
4582        *loghr = -LOG(r) + LOG(Dgamma(1./r,1./tune,tune)/Dgamma(r,1./tune,tune));
4583        break;
4584      }
4585
4586    default : 
4587      {
4588        PhyML_Printf("\n. Move not implemented");
4589        Exit("");
4590        break;
4591      }
4592    }
4593}
4594
4595//////////////////////////////////////////////////////////////
4596//////////////////////////////////////////////////////////////
4597
4598void MCMC_Read_Param_Vals(t_tree *tree)
4599{
4600  char *token;
4601  int sizemax;
4602  FILE *in_fp;
4603  phydbl val;
4604  int i,v;
4605 
4606  in_fp = tree->mcmc->in_fp_par;
4607 
4608
4609  sizemax = T_MAX_LINE;
4610
4611  token = (char *)mCalloc(sizemax,sizeof(char));
4612 
4613  if(fgets(token,sizemax,in_fp)); // Skip first line
4614  if(fgets(token,sizemax,in_fp)); // Skip second
4615 
4616  v=fscanf(in_fp,"%lf\t",&val); // Run
4617  /* PhyML_Printf("\n. Run = %d",(int)val); */
4618  v=fscanf(in_fp,"%lf\t",&val); // LnLike[Exact]
4619  /* PhyML_Printf("\n. LnLike = %f",val); */
4620  v=fscanf(in_fp,"%lf\t",&val); // LnLike[Approx]
4621  /* PhyML_Printf("\n. LnLike = %f",val); */
4622  v=fscanf(in_fp,"%lf\t",&val); // LnPriorRate
4623  /* PhyML_Printf("\n. LnPrior = %f",val); */
4624  v=fscanf(in_fp,"%lf\t",&val); // LnPriorTime
4625  /* PhyML_Printf("\n. LnPrior = %f",val); */
4626  v=fscanf(in_fp,"%lf\t",&val); // LnPosterior
4627  /* PhyML_Printf("\n. LnPost = %f",val); */
4628
4629  v=fscanf(in_fp,"%lf\t",&val); // ClockRate
4630  tree->rates->clock_r = val;
4631  /* PhyML_Printf("\n. Clock = %f",val); */
4632
4633  v=fscanf(in_fp,"%lf\t",&val); // EvolRate
4634
4635  v=fscanf(in_fp,"%lf\t",&val); // Nu
4636  tree->rates->nu = val;
4637  /* PhyML_Printf("\n. Nu = %f",val); */
4638
4639  v=fscanf(in_fp,"%lf\t",&val); // Birth rate
4640  tree->rates->birth_rate = val;
4641
4642  v=fscanf(in_fp,"%lf\t",&val); // TsTv
4643  tree->mod->kappa->v = val;
4644  /* PhyML_Printf("\n. TsTv = %f",val); */
4645
4646  For(i,tree->n_otu-1)
4647    {
4648      v=fscanf(in_fp,"%lf\t",&val); // Node heights
4649      tree->rates->nd_t[i+tree->n_otu] = val;
4650    }
4651
4652  For(i,2*tree->n_otu-2)
4653    {
4654      v=fscanf(in_fp,"%lf\t",&val); // Edge average rates
4655      tree->rates->br_r[i] = LOG(val);
4656    }
4657 
4658  v++;
4659
4660  Free(token);
4661
4662}
4663
4664//////////////////////////////////////////////////////////////
4665//////////////////////////////////////////////////////////////
4666
4667//////////////////////////////////////////////////////////////
4668//////////////////////////////////////////////////////////////
4669
4670//////////////////////////////////////////////////////////////
4671//////////////////////////////////////////////////////////////
4672
4673//////////////////////////////////////////////////////////////
4674//////////////////////////////////////////////////////////////
4675
4676//////////////////////////////////////////////////////////////
4677//////////////////////////////////////////////////////////////
4678
4679//////////////////////////////////////////////////////////////
4680//////////////////////////////////////////////////////////////
4681
Note: See TracBrowser for help on using the repository browser.