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

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

added most recent version of phyml
File size: 43.5 KB
Line 
1/*
2
3PhyML:  a program that  computes maximum likelihood phylogenies from
4DNA or AA homologous sequences.
5
6Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8All parts of the source except where indicated are distributed under
9the GNU public licence. See http://www.opensource.org for details.
10
11*/
12
13
14/* Routines for molecular clock trees and molecular dating */
15
16
17#include "times.h"
18
19//////////////////////////////////////////////////////////////
20//////////////////////////////////////////////////////////////
21
22#ifdef PHYTIME
23
24int TIMES_main(int argc, char **argv)
25{
26  align **data;
27  calign *cdata;
28  option *io;
29  t_tree *tree;
30  int num_data_set;
31  int num_tree,num_rand_tree;
32  t_mod *mod;
33  time_t t_beg,t_end;
34  int r_seed;
35  char *most_likely_tree;
36  int i;
37  int user_lk_approx;
38 
39#ifdef MPI
40  int rc;
41  rc = MPI_Init(&argc,&argv);
42  if (rc != MPI_SUCCESS) {
43    PhyML_Printf("\n. Error starting MPI program. Terminating.\n");
44    MPI_Abort(MPI_COMM_WORLD, rc);
45  }
46  MPI_Comm_size(MPI_COMM_WORLD,&Global_numTask);
47  MPI_Comm_rank(MPI_COMM_WORLD,&Global_myRank);
48#endif
49
50#ifdef QUIET
51  setvbuf(stdout,NULL,_IOFBF,2048);
52#endif
53
54  tree             = NULL;
55  mod              = NULL;
56  data             = NULL;
57  most_likely_tree = NULL;
58
59  io = (option *)Get_Input(argc,argv);
60  r_seed = (io->r_seed < 0)?(time(NULL)):(io->r_seed);
61  /* !!!!!!!!!!!!!!!!!!!!!!!! */
62  /* r_seed = 1289246338; */
63/*   r_seed = 1289266727; */
64/*   r_seed = 1289422815; */
65/*   r_seed = 1289443891; */
66/*   r_seed = 1290652518; */
67/*   r_seed = 1292195490; */
68  /* r_seed =  1298284669; */
69  /* r_seed  = 1298403366; */
70  /* r_seed = 1298509108; */
71  /* sys = system("sleep 5s"); */
72  /* r_seed = 1299649586; */
73  /* r_seed = 1302160422; */
74  /* r_seed = 1302576741; */
75  /* r_seed = 1302588678; */
76  /* r_seed = 1303247709; */
77  /* r_seed =  1303970631; */
78  /* r_seed = 1304059976; */
79  /* r_seed = 1306315195; */
80  /* r_seed = 1308263660; */
81  /* r_seed = 1313356025; */
82
83
84  /* phydbl mean; */
85  /* phydbl T,A,B,u; */
86  /* phydbl K; */
87
88  /* u = 1.E-4; */
89
90  /* K = 1.; */
91  /* T = 9 * K; */
92  /* A = LOG(1/K); */
93  /* B = LOG(2/K); */
94
95  /* Integrated_Geometric_Brownian_Bridge_Mean(T,A,B,u,&mean); */
96  /* printf("\n. mean = %f",mean*T); */
97
98  /* K = 1.E+2; */
99  /* T = 9 * K; */
100  /* A = LOG(1/K); */
101  /* B = LOG(2/K); */
102
103  /* Integrated_Geometric_Brownian_Bridge_Mean(T,A,B,u,&mean); */
104  /* printf("\n. mean = %f",mean*T); */
105
106
107  Exit("\n");
108
109
110
111  io->r_seed = r_seed;
112
113  srand(r_seed); rand();
114  PhyML_Printf("\n. Seed: %d\n",r_seed);
115  PhyML_Printf("\n. Pid: %d\n",getpid());
116  Make_Model_Complete(io->mod);
117  mod = io->mod;
118  if(io->in_tree == 2) Test_Multiple_Data_Set_Format(io);
119  else io->n_trees = 1;
120
121  if((io->n_data_sets > 1) && (io->n_trees > 1))
122    {
123      io->n_data_sets = MIN(io->n_trees,io->n_data_sets);
124      io->n_trees     = MIN(io->n_trees,io->n_data_sets);
125    }
126
127  For(num_data_set,io->n_data_sets)
128    {
129      data = Get_Seq(io);
130
131      if(data)
132        {
133          if(io->n_data_sets > 1) PhyML_Printf("\n. Data set [#%d]\n",num_data_set+1);
134          PhyML_Printf("\n. Compressing sequences...\n");
135          cdata = Compact_Data(data,io);
136          Free_Seq(data,cdata->n_otu);
137          Check_Ambiguities(cdata,io->mod->io->datatype,io->state_len);
138
139          for(num_tree=(io->n_trees == 1)?(0):(num_data_set);num_tree < io->n_trees;num_tree++)
140            {
141              if(!io->mod->s_opt->random_input_tree) io->mod->s_opt->n_rand_starts = 1;
142
143              For(num_rand_tree,io->mod->s_opt->n_rand_starts)
144                {
145                  if((io->mod->s_opt->random_input_tree) && (io->mod->s_opt->topo_search != NNI_MOVE))
146                    PhyML_Printf("\n. [Random start %3d/%3d]\n",num_rand_tree+1,io->mod->s_opt->n_rand_starts);
147
148
149                  Init_Model(cdata,mod,io);
150
151                  if(io->mod->use_m4mod) M4_Init_Model(mod->m4mod,cdata,mod);
152
153                  /* A BioNJ tree is built here */
154                  if(!io->in_tree) tree = Dist_And_BioNJ(cdata,mod,io);
155                  /* A user-given tree is used here instead of BioNJ */
156                  else             tree = Read_User_Tree(cdata,mod,io);
157
158
159                  if(io->fp_in_constraint_tree != NULL) 
160                    {
161                      io->cstr_tree        = Read_Tree_File_Phylip(io->fp_in_constraint_tree);               
162                      io->cstr_tree->rates = RATES_Make_Rate_Struct(io->cstr_tree->n_otu);
163                      RATES_Init_Rate_Struct(io->cstr_tree->rates,io->rates,io->cstr_tree->n_otu);
164                    }
165
166                  if(!tree) continue;
167
168                  if(!tree->n_root) 
169                    {
170                      PhyML_Printf("\n. Sorry, PhyTime requires a rooted tree as input.");
171                      Exit("\n");     
172                    }
173
174                  time(&t_beg);
175                  time(&(tree->t_beg));
176                  tree->rates = RATES_Make_Rate_Struct(tree->n_otu);
177                  RATES_Init_Rate_Struct(tree->rates,io->rates,tree->n_otu);
178
179                  Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
180                  Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);                 
181
182                  RATES_Fill_Lca_Table(tree);
183
184                  tree->mod         = mod;
185                  tree->io          = io;
186                  tree->data        = cdata;
187                  tree->n_pattern   = tree->data->crunch_len/tree->io->state_len;
188
189                  Set_Both_Sides(YES,tree);
190
191                  Prepare_Tree_For_Lk(tree);
192
193                  /* Read node age priors */
194                  Read_Clade_Priors(io->clade_list_file,tree);
195
196                  /* Set upper and lower bounds for all node ages */
197                  TIMES_Set_All_Node_Priors(tree);
198
199                  /* Count the number of time slices */
200                  TIMES_Get_Number_Of_Time_Slices(tree);
201                 
202                  TIMES_Label_Edges_With_Calibration_Intervals(tree);
203                  tree->write_br_lens = NO;
204                  PhyML_Printf("\n");
205                  PhyML_Printf("\n. Input tree with calibration information ('almost' compatible with MCMCtree).\n");
206                  PhyML_Printf("\n%s\n",Write_Tree(tree,YES));
207                  tree->write_br_lens = YES;
208
209
210                  /* Get_Edge_Binary_Coding_Number(tree); */
211                  /* Exit("\n"); */
212
213                  /* Print_CSeq_Select(stdout,NO,tree->data,tree); */
214                  /* Exit("\n"); */
215
216                  /* TIMES_Set_Root_Given_Tip_Dates(tree); */
217                  /* int i; */
218                  /* char *s; */
219                  /* FILE *fp; */
220                  /* For(i,2*tree->n_otu-2) tree->rates->cur_l[i] = 1.; */
221                  /* s = Write_Tree(tree,NO); */
222                  /* fp = fopen("rooted_tree","w"); */
223                  /* fprintf(fp,"%s\n",s); */
224                  /* fclose(fp); */
225                  /* Exit("\n"); */
226
227
228                  /* Work with log of branch lengths? */
229                  if(tree->mod->log_l == YES) Log_Br_Len(tree);
230
231                  if(io->mcmc->use_data == YES)
232                    {
233                      /* Force the exact likelihood score */
234                      user_lk_approx = tree->io->lk_approx;
235                      tree->io->lk_approx = EXACT;
236                     
237                      /* printf("\n. %s",Write_Tree(tree,NO)); */
238                      /* Lk(NULL,tree); */
239                      /* printf("\n. %f",tree->c_lnL); */
240                      /* Exit("\n"); */
241
242                      /* MLE for branch lengths */
243                      PhyML_Printf("\n");
244                      Round_Optimize(tree,tree->data,ROUND_MAX);
245                     
246                      /* Set vector of mean branch lengths for the Normal approximation
247                         of the likelihood */
248                      RATES_Set_Mean_L(tree);
249                     
250                      /* Estimate the matrix of covariance for the Normal approximation of
251                         the likelihood */
252                      PhyML_Printf("\n");
253                      PhyML_Printf("\n. Computing Hessian...");
254                      tree->rates->bl_from_rt = 0;
255                      Free(tree->rates->cov_l);
256                      tree->rates->cov_l = Hessian_Seo(tree);
257                      /* tree->rates->cov_l = Hessian_Log(tree); */
258                      For(i,(2*tree->n_otu-3)*(2*tree->n_otu-3)) tree->rates->cov_l[i] *= -1.0;
259                      if(!Iter_Matinv(tree->rates->cov_l,2*tree->n_otu-3,2*tree->n_otu-3,YES)) Exit("\n");
260                      tree->rates->covdet = Matrix_Det(tree->rates->cov_l,2*tree->n_otu-3,YES);
261                      For(i,(2*tree->n_otu-3)*(2*tree->n_otu-3)) tree->rates->invcov[i] = tree->rates->cov_l[i];
262                      if(!Iter_Matinv(tree->rates->invcov,2*tree->n_otu-3,2*tree->n_otu-3,YES)) Exit("\n");
263                      tree->rates->grad_l = Gradient(tree);
264                     
265                     
266                      /* Pre-calculation of conditional variances to speed up calculations */
267                      RATES_Bl_To_Ml(tree);
268                      RATES_Get_Conditional_Variances(tree);
269                      RATES_Get_All_Reg_Coeff(tree);
270                      RATES_Get_Trip_Conditional_Variances(tree);
271                      RATES_Get_All_Trip_Reg_Coeff(tree);
272                     
273                      Lk(NULL,tree);
274                      PhyML_Printf("\n");
275                      PhyML_Printf("\n. p(data|model) [exact ] ~ %.2f",tree->c_lnL);
276                     
277                      tree->io->lk_approx = NORMAL;
278                      For(i,2*tree->n_otu-3) tree->rates->u_cur_l[i] = tree->rates->mean_l[i] ;
279                      tree->c_lnL = Lk(NULL,tree);
280                      PhyML_Printf("\n. p(data|model) [approx] ~ %.2f",tree->c_lnL);
281
282                      tree->io->lk_approx = user_lk_approx;
283                    }
284                 
285                     
286                  tree->rates->model = io->rates->model;                 
287                  PhyML_Printf("\n. Selected model '%s'",RATES_Get_Model_Name(io->rates->model));
288                  if(tree->rates->model == GUINDON) tree->mod->gamma_mgf_bl = YES;
289                 
290                  tree->rates->bl_from_rt = YES;
291                 
292                  if(tree->io->cstr_tree) Find_Surviving_Edges_In_Small_Tree(tree,tree->io->cstr_tree);
293
294                  time(&t_beg);
295                  tree->mcmc = MCMC_Make_MCMC_Struct();
296                  MCMC_Copy_MCMC_Struct(tree->io->mcmc,tree->mcmc,"phytime");
297                  MCMC_Complete_MCMC(tree->mcmc,tree);
298                  tree->mcmc->is_burnin = NO;
299                  MCMC(tree);
300                  MCMC_Close_MCMC(tree->mcmc);
301                  MCMC_Free_MCMC(tree->mcmc);
302                  PhyML_Printf("\n");
303
304/*                tree->mcmc = MCMC_Make_MCMC_Struct(); */
305/*                MCMC_Copy_MCMC_Struct(tree->io->mcmc,tree->mcmc,"burnin"); */
306/*                MCMC_Complete_MCMC(tree->mcmc,tree); */
307/*                tree->mcmc->adjust_tuning = YES; */
308/*                tree->mcmc->is_burnin     = YES; */
309/*                tree->mcmc->chain_len = tree->io->mcmc->chain_len_burnin; */
310/*                MCMC(tree); */
311/*                MCMC_Close_MCMC(tree->mcmc); */
312
313                 
314/*                new_mcmc = MCMC_Make_MCMC_Struct(tree); */
315/*                MCMC_Complete_MCMC(new_mcmc,tree); */
316/*                MCMC_Copy_MCMC_Struct(tree->mcmc,new_mcmc,"phytime"); */
317/*                MCMC_Free_MCMC(tree->mcmc); */
318
319/*                tree->mcmc                  = new_mcmc; */
320/*                tree->mcmc->chain_len       = tree->io->mcmc->chain_len; */
321/*                tree->mcmc->randomize       = NO; */
322/*                tree->mcmc->adjust_tuning   = NO; */
323/*                tree->mcmc->is_burnin       = NO; */
324
325/*                time(&t_beg); */
326/*                MCMC(tree); */
327/*                MCMC_Close_MCMC(tree->mcmc); */
328/*                MCMC_Free_MCMC(tree->mcmc); */
329/*                PhyML_Printf("\n"); */
330
331
332                  Free_Tree_Pars(tree);
333                  Free_Tree_Lk(tree);
334                  Free_Tree(tree);
335                }
336
337              break;
338            }
339          Free_Cseq(cdata);
340        }
341    }
342
343  Free_Model(mod);
344
345  if(io->fp_in_align)   fclose(io->fp_in_align);
346  if(io->fp_in_tree)    fclose(io->fp_in_tree);
347  if(io->fp_out_lk)     fclose(io->fp_out_lk);
348  if(io->fp_out_tree)   fclose(io->fp_out_tree);
349  if(io->fp_out_trees)  fclose(io->fp_out_trees);
350  if(io->fp_out_stats)  fclose(io->fp_out_stats);
351
352  Free(most_likely_tree);
353  Free_Input(io);
354
355  time(&t_end);
356  Print_Time_Info(t_beg,t_end);
357
358#ifdef MPI
359  MPI_Finalize();
360#endif
361
362  return 0;
363}
364
365#endif
366
367//////////////////////////////////////////////////////////////
368//////////////////////////////////////////////////////////////
369
370
371void TIMES_Least_Square_Node_Times(t_edge *e_root, t_tree *tree)
372{
373
374  /* Solve A.x = b, where x are the t_node time estimated
375     under the least square criterion.
376
377     A is a n x n matrix, with n being the number of
378     nodes in a rooted tree (i.e. 2*n_otu-1).
379
380   */
381
382  phydbl *A, *b, *x;
383  int n;
384  int i,j;
385  t_node *root;
386
387  n = 2*tree->n_otu-1;
388 
389  A = (phydbl *)mCalloc(n*n,sizeof(phydbl));
390  b = (phydbl *)mCalloc(n,  sizeof(phydbl));
391  x = (phydbl *)mCalloc(n,  sizeof(phydbl));
392   
393  if(!tree->n_root && e_root) Add_Root(e_root,tree);
394  else if(!e_root)            Add_Root(tree->a_edges[0],tree);
395 
396  root = tree->n_root;
397
398  TIMES_Least_Square_Node_Times_Pre(root,root->v[2],A,b,n,tree);
399  TIMES_Least_Square_Node_Times_Pre(root,root->v[1],A,b,n,tree);
400 
401  b[root->num] = tree->e_root->l->v/2.;
402 
403  A[root->num * n + root->num]       = 1.0;
404  A[root->num * n + root->v[2]->num] = -.5;
405  A[root->num * n + root->v[1]->num] = -.5;
406   
407  if(!Matinv(A, n, n,YES))
408    {
409      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
410      Exit("\n");     
411    }
412
413  For(i,n) x[i] = .0;
414  For(i,n) For(j,n) x[i] += A[i*n+j] * b[j];
415
416  For(i,n-1) { tree->rates->nd_t[tree->a_nodes[i]->num] = -x[i]; }
417  tree->rates->nd_t[root->num] = -x[n-1];
418  tree->n_root->l[2] = tree->rates->nd_t[root->v[2]->num] - tree->rates->nd_t[root->num];
419  tree->n_root->l[1] = tree->rates->nd_t[root->v[1]->num] - tree->rates->nd_t[root->num];
420  ////////////////////////////////////////
421  return;
422  ////////////////////////////////////////
423
424  /* Rescale the t_node times such that the time at the root
425     is -100. This constraint implies that the clock rate
426     is fixed to the actual tree length divided by the tree
427     length measured in term of differences of t_node times */
428
429  phydbl scale_f,time_tree_length,tree_length;
430
431  scale_f = -100./tree->rates->nd_t[root->num];
432  For(i,2*tree->n_otu-1) tree->rates->nd_t[i] *= scale_f;
433  For(i,2*tree->n_otu-1) if(tree->rates->nd_t[i] > .0) tree->rates->nd_t[i] = .0;
434
435  time_tree_length = 0.0;
436  For(i,2*tree->n_otu-3)
437    if(tree->a_edges[i] != tree->e_root)
438      time_tree_length +=
439        FABS(tree->rates->nd_t[tree->a_edges[i]->left->num] -
440             tree->rates->nd_t[tree->a_edges[i]->rght->num]);
441  time_tree_length += FABS(tree->rates->nd_t[root->num] - tree->rates->nd_t[root->v[2]->num]);
442  time_tree_length += FABS(tree->rates->nd_t[root->num] - tree->rates->nd_t[root->v[1]->num]);
443 
444  tree_length = 0.0;
445  For(i,2*tree->n_otu-3) tree_length += tree->a_edges[i]->l->v;
446
447  tree->rates->clock_r = tree_length / time_tree_length;
448
449  Free(A);
450  Free(b);
451  Free(x);
452
453}
454
455//////////////////////////////////////////////////////////////
456//////////////////////////////////////////////////////////////
457
458
459void TIMES_Least_Square_Node_Times_Pre(t_node *a, t_node *d, phydbl *A, phydbl *b, int n, t_tree *tree)
460{
461  if(d->tax)
462    {
463      A[d->num * n + d->num] = 1.;
464     
465      /* Set the time stamp at tip nodes to 0.0 */
466/*       PhyML_Printf("\n. Tip t_node date set to 0"); */
467      b[d->num] = 0.0;
468      return;
469    }
470  else
471    {
472      int i;
473     
474 
475      For(i,3)
476        if((d->v[i] != a) && (d->b[i] != tree->e_root))
477          TIMES_Least_Square_Node_Times_Pre(d,d->v[i],A,b,n,tree);
478     
479      A[d->num * n + d->num] = 1.;
480      b[d->num] = .0;
481      For(i,3)
482        {
483          A[d->num * n + d->v[i]->num] = -1./3.;
484          if(d->v[i] != a) b[d->num] += d->b[i]->l->v;
485          else             b[d->num] -= d->b[i]->l->v;
486        }
487      b[d->num] /= 3.;
488    }
489}
490
491//////////////////////////////////////////////////////////////
492//////////////////////////////////////////////////////////////
493
494
495/* Adjust t_node times in order to have correct time stamp ranking with
496 respect to the tree topology */
497
498void TIMES_Adjust_Node_Times(t_tree *tree)
499{
500  TIMES_Adjust_Node_Times_Pre(tree->n_root->v[2],tree->n_root->v[1],tree);
501  TIMES_Adjust_Node_Times_Pre(tree->n_root->v[1],tree->n_root->v[2],tree);
502
503  if(tree->rates->nd_t[tree->n_root->num] > MIN(tree->rates->nd_t[tree->n_root->v[2]->num],
504                                                tree->rates->nd_t[tree->n_root->v[1]->num]))
505    {
506      tree->rates->nd_t[tree->n_root->num] = MIN(tree->rates->nd_t[tree->n_root->v[2]->num],
507                                                 tree->rates->nd_t[tree->n_root->v[1]->num]);
508    }
509}
510
511//////////////////////////////////////////////////////////////
512//////////////////////////////////////////////////////////////
513
514
515void TIMES_Adjust_Node_Times_Pre(t_node *a, t_node *d, t_tree *tree)
516{
517  if(d->tax) return;
518  else
519    {
520      int i;
521      phydbl min_height;
522
523      For(i,3)
524        if((d->v[i] != a) && (d->b[i] != tree->e_root))
525          {
526            TIMES_Adjust_Node_Times_Pre(d,d->v[i],tree);
527          }
528
529      min_height = 0.0;
530      For(i,3)
531        {
532          if((d->v[i] != a) && (d->b[i] != tree->e_root))
533            {
534              if(tree->rates->nd_t[d->v[i]->num] < min_height)
535                {
536                  min_height = tree->rates->nd_t[d->v[i]->num];
537                }
538            }
539        }
540
541      if(tree->rates->nd_t[d->num] > min_height) tree->rates->nd_t[d->num] = min_height;
542
543      if(tree->rates->nd_t[d->num] < -100.) tree->rates->nd_t[d->num] = -100.;
544
545    }
546}
547
548//////////////////////////////////////////////////////////////
549//////////////////////////////////////////////////////////////
550
551
552  /* Multiply each time stamp at each internal
553     t_node by  'tree->time_stamp_mult'.
554   */
555
556void TIMES_Mult_Time_Stamps(t_tree *tree)
557{
558  int i;
559  For(i,2*tree->n_otu-2) tree->rates->nd_t[tree->a_nodes[i]->num] *= FABS(tree->mod->s_opt->tree_size_mult);
560  tree->rates->nd_t[tree->n_root->num] *= FABS(tree->mod->s_opt->tree_size_mult);
561}
562
563//////////////////////////////////////////////////////////////
564//////////////////////////////////////////////////////////////
565
566
567void TIMES_Print_Node_Times(t_node *a, t_node *d, t_tree *tree)
568{
569  t_edge *b;
570  int i;
571 
572  b = NULL;
573  For(i,3) if((d->v[i]) && (d->v[i] == a)) {b = d->b[i]; break;}
574
575  PhyML_Printf("\n. (%3d %3d) a->t = %12f d->t = %12f (#=%12f) b->l->v = %12f [%12f;%12f]",
576               a->num,d->num,
577               tree->rates->nd_t[a->num],
578               tree->rates->nd_t[d->num],
579               tree->rates->nd_t[a->num]-tree->rates->nd_t[d->num],
580               (b)?(b->l->v):(-1.0),
581               tree->rates->t_prior_min[d->num],
582               tree->rates->t_prior_max[d->num]);
583  if(d->tax) return;
584  else
585    {
586      int i;
587      For(i,3)
588        if((d->v[i] != a) && (d->b[i] != tree->e_root))
589          TIMES_Print_Node_Times(d,d->v[i],tree);
590    }
591}
592
593//////////////////////////////////////////////////////////////
594//////////////////////////////////////////////////////////////
595
596
597void TIMES_Set_All_Node_Priors(t_tree *tree)
598{
599  int i;
600  phydbl min_prior;
601
602  /* Set all t_prior_max values */
603  TIMES_Set_All_Node_Priors_Bottom_Up(tree->n_root,tree->n_root->v[2],tree);
604  TIMES_Set_All_Node_Priors_Bottom_Up(tree->n_root,tree->n_root->v[1],tree);
605
606  tree->rates->t_prior_max[tree->n_root->num] = 
607    MIN(tree->rates->t_prior_max[tree->n_root->num],
608        MIN(tree->rates->t_prior_max[tree->n_root->v[2]->num],
609            tree->rates->t_prior_max[tree->n_root->v[1]->num]));
610
611
612  /* Set all t_prior_min values */
613  if(!tree->rates->t_has_prior[tree->n_root->num])
614    {
615      min_prior = 1.E+10;
616      For(i,2*tree->n_otu-2)
617        {
618          if(tree->rates->t_has_prior[i])
619            {
620              if(tree->rates->t_prior_min[i] < min_prior)
621                min_prior = tree->rates->t_prior_min[i];
622            }
623        }
624      tree->rates->t_prior_min[tree->n_root->num] = 2.0 * min_prior;
625      /* tree->rates->t_prior_min[tree->n_root->num] = 10. * min_prior; */
626    }
627 
628  if(tree->rates->t_prior_min[tree->n_root->num] > 0.0)
629    {
630      PhyML_Printf("\n== Failed to set the lower bound for the root node.");
631      PhyML_Printf("\n== Make sure at least one of the calibration interval");
632      PhyML_Printf("\n== provides a lower bound.");
633      Exit("\n");
634    }
635
636
637  TIMES_Set_All_Node_Priors_Top_Down(tree->n_root,tree->n_root->v[2],tree);
638  TIMES_Set_All_Node_Priors_Top_Down(tree->n_root,tree->n_root->v[1],tree);
639
640  Get_Node_Ranks(tree);
641  TIMES_Set_Floor(tree);
642}
643
644//////////////////////////////////////////////////////////////
645//////////////////////////////////////////////////////////////
646
647
648void TIMES_Set_All_Node_Priors_Bottom_Up(t_node *a, t_node *d, t_tree *tree)
649{
650  int i;
651  phydbl t_sup;
652 
653  if(d->tax) return;
654  else 
655    {
656      t_node *v1, *v2; /* the two sons of d */
657
658      For(i,3)
659        {
660          if((d->v[i] != a) && (d->b[i] != tree->e_root))
661            {
662              TIMES_Set_All_Node_Priors_Bottom_Up(d,d->v[i],tree);           
663            }
664        }
665     
666      v1 = v2 = NULL;
667      For(i,3) if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
668        {
669          if(!v1) v1 = d->v[i]; 
670          else    v2 = d->v[i];
671        }
672     
673      if(tree->rates->t_has_prior[d->num] == YES)
674        {
675          t_sup = MIN(tree->rates->t_prior_max[d->num],
676                      MIN(tree->rates->t_prior_max[v1->num],
677                          tree->rates->t_prior_max[v2->num]));
678
679          tree->rates->t_prior_max[d->num] = t_sup;
680
681          if(tree->rates->t_prior_max[d->num] < tree->rates->t_prior_min[d->num])
682            {
683              PhyML_Printf("\n. prior_min=%f prior_max=%f",tree->rates->t_prior_min[d->num],tree->rates->t_prior_max[d->num]);
684              PhyML_Printf("\n. Inconsistency in the prior settings detected at t_node %d",d->num);
685              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
686              Warn_And_Exit("\n");
687            }
688        }
689      else
690        {
691          tree->rates->t_prior_max[d->num] = 
692            MIN(tree->rates->t_prior_max[v1->num],
693                tree->rates->t_prior_max[v2->num]);
694        }
695    }
696}
697
698//////////////////////////////////////////////////////////////
699//////////////////////////////////////////////////////////////
700
701
702void TIMES_Set_All_Node_Priors_Top_Down(t_node *a, t_node *d, t_tree *tree)
703{
704  if(d->tax) return;
705  else
706    {
707      int i;     
708     
709      if(tree->rates->t_has_prior[d->num] == YES)
710        {
711          tree->rates->t_prior_min[d->num] = MAX(tree->rates->t_prior_min[d->num],tree->rates->t_prior_min[a->num]);
712         
713          if(tree->rates->t_prior_max[d->num] < tree->rates->t_prior_min[d->num])
714            {
715              PhyML_Printf("\n. prior_min=%f prior_max=%f",tree->rates->t_prior_min[d->num],tree->rates->t_prior_max[d->num]);
716              PhyML_Printf("\n. Inconsistency in the prior settings detected at t_node %d",d->num);
717              PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
718              Warn_And_Exit("\n");
719            }
720        }
721      else
722        {
723          tree->rates->t_prior_min[d->num] = tree->rates->t_prior_min[a->num];
724        }
725           
726      For(i,3)
727        {
728          if((d->v[i] != a) && (d->b[i] != tree->e_root))
729            {
730              TIMES_Set_All_Node_Priors_Top_Down(d,d->v[i],tree);
731            }
732        }
733    }
734}
735
736//////////////////////////////////////////////////////////////
737//////////////////////////////////////////////////////////////
738
739
740void TIMES_Set_Floor(t_tree *tree)
741{
742  TIMES_Set_Floor_Post(tree->n_root,tree->n_root->v[2],tree);
743  TIMES_Set_Floor_Post(tree->n_root,tree->n_root->v[1],tree);
744  tree->rates->t_floor[tree->n_root->num] = MIN(tree->rates->t_floor[tree->n_root->v[2]->num],
745                                                tree->rates->t_floor[tree->n_root->v[1]->num]);
746}
747
748//////////////////////////////////////////////////////////////
749//////////////////////////////////////////////////////////////
750
751
752void TIMES_Set_Floor_Post(t_node *a, t_node *d, t_tree *tree)
753{
754  if(d->tax)
755    {
756      tree->rates->t_floor[d->num] = tree->rates->nd_t[d->num];
757      d->rank_max = d->rank;
758      return;
759    }
760  else
761    {
762      int i;
763      t_node *v1,*v2;
764
765      v1 = v2 = NULL;
766      For(i,3)
767        {
768          if(d->v[i] != a && d->b[i] != tree->e_root)
769            {
770              TIMES_Set_Floor_Post(d,d->v[i],tree);
771              if(!v1) v1 = d->v[i];
772              else    v2 = d->v[i];
773            }
774        }
775      tree->rates->t_floor[d->num] = MIN(tree->rates->t_floor[v1->num],
776                                         tree->rates->t_floor[v2->num]);
777
778      if(tree->rates->t_floor[v1->num] < tree->rates->t_floor[v2->num])
779        {
780          d->rank_max = v1->rank_max;
781        }
782      else if(tree->rates->t_floor[v2->num] < tree->rates->t_floor[v1->num])
783        {
784          d->rank_max = v2->rank_max;
785        }
786      else
787        {
788          d->rank_max = MAX(v1->rank_max,v2->rank_max);
789        }
790    }
791}
792
793//////////////////////////////////////////////////////////////
794//////////////////////////////////////////////////////////////
795
796/* Does it work for serial samples? */
797phydbl TIMES_Log_Conditional_Uniform_Density(t_tree *tree)
798{
799  phydbl min,max;
800  phydbl dens;
801  int i;
802
803  min = tree->rates->nd_t[tree->n_root->num];
804
805  dens = 0.0;
806  For(i,2*tree->n_otu-1)
807    {
808      if((tree->a_nodes[i]->tax == NO) && (tree->a_nodes[i] != tree->n_root))
809        {
810          max = tree->rates->t_floor[i];
811
812          dens += LOG(Dorder_Unif(tree->rates->nd_t[i],
813                                  tree->a_nodes[i]->rank-1,
814                                  tree->a_nodes[i]->rank_max-2,
815                                  min,max));
816        }
817    }
818  return dens;
819}
820
821//////////////////////////////////////////////////////////////
822//////////////////////////////////////////////////////////////
823// Returns the marginal density of tree height assuming the
824// Yule model of speciation.
825phydbl TIMES_Lk_Yule_Root_Marginal(t_tree *tree)
826{
827  int n;
828  int j;
829  t_node *nd;
830  phydbl *t,*ts;
831  phydbl lbda;
832  phydbl T;
833
834  lbda = tree->rates->birth_rate;
835  t    = tree->rates->nd_t;
836  ts   = tree->rates->time_slice_lims;
837  T    = ts[0] - t[tree->n_root->num];
838
839  n = 0;
840  nd = NULL;
841  For(j,2*tree->n_otu-2) 
842    {
843      nd = tree->a_nodes[j];
844
845      if((t[nd->num] > ts[0] && t[nd->anc->num] < ts[0]) || // lineage that is crossing ts[0]
846         (nd->tax == YES && Are_Equal(t[nd->num],ts[0],1.E-6) == YES)) // tip that is lying on ts[0]
847        n++;
848    }
849
850  return LnGamma(n+1) + LOG(lbda) - 2.*lbda*T + (n-2.)*LOG(1. - EXP(-lbda*T));
851}
852
853//////////////////////////////////////////////////////////////
854//////////////////////////////////////////////////////////////
855// Returns the joint density of internal node heights assuming
856// the Yule model of speciation.
857phydbl TIMES_Lk_Yule_Joint(t_tree *tree)
858{
859  int i,j;
860  phydbl loglk;
861  phydbl *t;
862  phydbl dt;
863  int n; // number of lineages at a given time point
864  phydbl lbda;
865  t_node *nd;
866  phydbl *ts;
867  int *tr;
868  phydbl top_t;
869  short int *interrupted;
870  phydbl sumdt;
871
872  interrupted = (short int *)mCalloc(tree->n_otu,sizeof(short int));
873
874  t = tree->rates->nd_t;
875  ts = tree->rates->time_slice_lims;
876  tr = tree->rates->t_ranked;
877  lbda = tree->rates->birth_rate;
878
879  TIMES_Update_Node_Ordering(tree);
880
881  For(j,tree->n_otu) interrupted[j] = NO;
882
883  loglk = .0;
884
885  sumdt = .0;
886  n = 1;
887  For(i,2*tree->n_otu-2) // t[tr[0]] is the oldest node, t[tr[1]], the second oldest and so on...
888    {
889
890      For(j,tree->n_otu)
891        if((t[j] < t[tr[i]]) && (interrupted[j] == NO)) 
892          {
893            interrupted[j] = YES;
894            n--; // How many lineages have stopped above t[tr[i]]?
895          }
896     
897      top_t = t[tr[i+1]];
898      dt = top_t - t[tr[i]];
899      sumdt += dt;
900
901      /* printf("\n. %d node up=%d [%f] node do=%d [%f] dt=%f",i,tr[i],t[tr[i]],tr[i+1],t[tr[i+1]],dt); */
902
903      if(n<1)
904        {
905          PhyML_Printf("\n. i=%d tr[i]=%f",i,t[tr[i]]);
906          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
907          Exit("\n");
908        }
909
910      if(dt > 1.E-10) loglk += LOG((n+1)*lbda) - (n+1)*lbda*dt;
911      n++;     
912    }
913
914  /* printf("\n. sumdt = %f th=%f",sumdt,tree->rates->nd_t[tree->n_root->num]); */
915  /* printf("\n0 loglk = %f",loglk); */
916
917  For(i,tree->rates->n_time_slices-1)
918    {
919      n = 0;
920      dt = 0.;
921      For(j,2*tree->n_otu-2)
922        {
923          nd = tree->a_nodes[j];
924          if(t[nd->num] > ts[i] && t[nd->anc->num] < ts[i]) // How many lineages are crossing this time slice limit?
925            {
926              n++;
927              if(t[nd->num] < dt) dt = t[nd->num]; // take the oldest node younger than the time slice
928            }
929        }
930      dt -= ts[i];
931      loglk += LOG(n*lbda) - n*lbda*dt;
932    }
933
934  /* printf("\n1 loglk = %f",loglk); */
935
936  Free(interrupted);
937
938  return loglk;
939}
940
941//////////////////////////////////////////////////////////////
942//////////////////////////////////////////////////////////////
943// Returns the conditional density of internal node heights
944// given the tree height under the Yule model. Uses the order
945// statistics 'simplification' as described in Yang and Rannala, 2005.
946phydbl TIMES_Lk_Yule_Order(t_tree *tree)
947{
948
949  int j;
950  phydbl *t,*tf;
951  t_node *n;
952  phydbl loglk;
953  phydbl loglbda;
954  phydbl lbda;
955  phydbl *tp_min,*tp_max;
956  phydbl lower_bound,upper_bound;
957
958  tp_min = tree->rates->t_prior_min;
959  tp_max = tree->rates->t_prior_max;
960  tf = tree->rates->t_floor;
961  t  = tree->rates->nd_t;
962  n = NULL;
963  loglbda = LOG(tree->rates->birth_rate);
964  lbda = tree->rates->birth_rate;
965  lower_bound = -1.;
966  upper_bound = -1.;
967
968  /*! Adapted from  Equation (6) in T. Stadler's Systematic Biology, 2012 paper with
969      sampling fraction set to 1 and death rate set to 0. Dropped the 1/(n-1) scaling
970      factor. */
971
972  /* loglk = 0.0; */
973  /* For(j,2*tree->n_otu-2) */
974  /*   { */
975  /*     n = tree->a_nodes[j]; */
976  /*     lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j])); */
977  /*     upper_bound = MIN(FABS(t[tree->n_root->num]),FABS(tp_min[j])); */
978
979  /*     if(n->tax == NO) */
980  /*       { */
981  /*         loglk  += (loglbda - lbda * FABS(t[j])); */
982  /*         /\* loglk -= LOG(EXP(-lbda*lower_bound) - EXP(-lbda*upper_bound)); // incorporate calibration boundaries here. *\/ */
983  /*       } */
984  /*   } */
985
986 
987  /*! Adapted from  Equation (7) in T. Stadler's Systematic Biology, 2012 paper with
988      sampling fraction set to 1 and death rate set to 0. */
989
990  // Check that each node is within its calibration-derived interval
991  For(j,2*tree->n_otu-1) if(t[j] < tp_min[j] || t[j] > tp_max[j]) return(-INFINITY);
992
993  loglk = 0.0;
994  For(j,2*tree->n_otu-2)
995    {
996      n = tree->a_nodes[j];
997      lower_bound = MAX(FABS(tf[j]),FABS(tp_max[j]));
998      upper_bound = FABS(tp_min[j]);
999     
1000      if(n->tax == NO)
1001        {
1002          loglk  += (loglbda - lbda * FABS(t[j]));
1003          loglk -= LOG(EXP(-lbda*lower_bound) - EXP(-lbda*upper_bound)); // incorporate calibration boundaries here.
1004   
1005        }
1006    }
1007
1008  lower_bound = MAX(FABS(tf[tree->n_root->num]),FABS(tp_max[tree->n_root->num]));
1009  upper_bound = FABS(tp_min[tree->n_root->num]);
1010  loglk += LOG(2) + loglbda - 2.*lbda * FABS(t[tree->n_root->num]);
1011  loglk -= LOG(EXP(-2.*lbda*lower_bound) - EXP(-2.*lbda*upper_bound));
1012
1013  return(loglk);
1014}
1015
1016//////////////////////////////////////////////////////////////
1017//////////////////////////////////////////////////////////////
1018
1019phydbl TIMES_Lk_Times(t_tree *tree)
1020{
1021
1022  #ifdef PHYTIME
1023  tree->rates->c_lnL_times =  TIMES_Lk_Yule_Order(tree);
1024  #elif SERGEII
1025  tree->rates->c_lnL_times = TIMES_Calib_Cond_Prob(tree);
1026  //tree->rates->c_lnL_times =  TIMES_Lk_Yule_Order(tree);
1027  #endif
1028
1029
1030  if(isinf(tree->rates->c_lnL_times))
1031    {
1032      tree->rates->c_lnL_times = -INFINITY;
1033    }
1034
1035  return(tree->rates->c_lnL_times);
1036}
1037
1038//////////////////////////////////////////////////////////////
1039//////////////////////////////////////////////////////////////
1040
1041
1042void TIMES_Lk_Times_Trav(t_node *a, t_node *d, phydbl lim_inf, phydbl lim_sup, phydbl *logdens, t_tree *tree)
1043{
1044  int i;
1045 
1046  if(!d->tax)
1047    {
1048      /* if(tree->rates->nd_t[d->num] > lim_sup) */
1049      /*        { */
1050      /*          lim_inf = lim_sup; */
1051      /*          lim_sup = 0.0; */
1052      /*          For(i,2*tree->n_otu-2) */
1053      /*            if((tree->rates->t_floor[i] < lim_sup) && (tree->rates->t_floor[i] > tree->rates->nd_t[d->num])) */
1054      /*              lim_sup = tree->rates->t_floor[i]; */
1055      /*        } */
1056     
1057      /* if(tree->rates->nd_t[d->num] < lim_inf || tree->rates->nd_t[d->num] > lim_sup) */
1058      /*        { */
1059      /*          PhyML_Printf("\n. nd_t = %f lim_inf = %f lim_sup = %f",tree->rates->nd_t[d->num],lim_inf,lim_sup); */
1060      /*          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
1061      /*          Exit("\n"); */
1062      /*        } */
1063 
1064      lim_inf = tree->rates->nd_t[tree->n_root->num];
1065      lim_sup = tree->rates->t_floor[d->num];
1066     
1067      *logdens = *logdens + LOG(lim_sup - lim_inf);   
1068    }
1069 
1070  if(d->tax == YES) return;
1071  else
1072    {     
1073      For(i,3)
1074        {
1075          if(d->v[i] != a && d->b[i] != tree->e_root)
1076            {
1077              TIMES_Lk_Times_Trav(d,d->v[i],lim_inf,lim_sup,logdens,tree);
1078            }
1079        }
1080    }
1081}
1082
1083//////////////////////////////////////////////////////////////
1084//////////////////////////////////////////////////////////////
1085
1086
1087phydbl TIMES_Log_Number_Of_Ranked_Labelled_Histories(t_node *root, int per_slice, t_tree *tree)
1088{
1089  int i;
1090  phydbl logn;
1091  t_node *v1,*v2;
1092  int n1,n2;
1093 
1094  TIMES_Update_Curr_Slice(tree);
1095
1096  logn = .0;
1097  v1 = v2 = NULL;
1098  if(root == tree->n_root)
1099    {
1100      TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[2],per_slice,&logn,tree);
1101      TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[1],per_slice,&logn,tree);
1102      v1 = root->v[2];
1103      v2 = root->v[1];
1104    }
1105  else
1106    {
1107      For(i,3)
1108        {
1109          if(root->v[i] != root->anc && root->b[i] != tree->e_root)
1110            {
1111              TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(root,root->v[i],per_slice,&logn,tree);
1112              if(!v1) v1 = root->v[i];
1113              else    v2 = root->v[i];
1114            }
1115        }
1116    }
1117
1118 
1119  if(per_slice == NO)
1120    {
1121      n1 = tree->rates->n_tips_below[v1->num];
1122      n2 = tree->rates->n_tips_below[v2->num];
1123    }
1124  else
1125    {
1126      if(tree->rates->curr_slice[v1->num] == tree->rates->curr_slice[root->num])
1127        n1 = tree->rates->n_tips_below[v1->num];
1128      else
1129        n1 = 1;
1130     
1131      if(tree->rates->curr_slice[v2->num] == tree->rates->curr_slice[root->num])
1132        n2 = tree->rates->n_tips_below[v2->num];
1133      else
1134        n2 = 1;
1135    }
1136
1137  tree->rates->n_tips_below[root->num] = n1+n2;
1138
1139  logn += Factln(n1+n2-2) - Factln(n1-1) - Factln(n2-1);
1140
1141  return(logn);
1142}
1143
1144//////////////////////////////////////////////////////////////
1145//////////////////////////////////////////////////////////////
1146
1147
1148void TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(t_node *a, t_node *d, int per_slice, phydbl *logn, t_tree *tree)
1149{
1150  if(d->tax == YES) 
1151    {
1152      tree->rates->n_tips_below[d->num] = 1;
1153      return;
1154    }
1155  else
1156    {
1157      int i,n1,n2;
1158      t_node *v1, *v2;
1159
1160      For(i,3)
1161        {
1162          if(d->v[i] != a && d->b[i] != tree->e_root)
1163            {
1164              TIMES_Log_Number_Of_Ranked_Labelled_Histories_Post(d,d->v[i],per_slice,logn,tree);
1165            }
1166        }
1167
1168      v1 = v2 = NULL;
1169      For(i,3)
1170        {
1171          if(d->v[i] != a && d->b[i] != tree->e_root)
1172            {
1173              if(v1 == NULL) {v1 = d->v[i];}
1174              else           {v2 = d->v[i];}
1175            }
1176        }
1177
1178
1179      if(per_slice == NO)
1180        {
1181          n1 = tree->rates->n_tips_below[v1->num];
1182          n2 = tree->rates->n_tips_below[v2->num];
1183        }
1184      else
1185        {
1186          if(tree->rates->curr_slice[v1->num] == tree->rates->curr_slice[d->num])
1187            n1 = tree->rates->n_tips_below[v1->num];
1188          else
1189            n1 = 1;
1190
1191          if(tree->rates->curr_slice[v2->num] == tree->rates->curr_slice[d->num])
1192            n2 = tree->rates->n_tips_below[v2->num];
1193          else
1194            n2 = 1;
1195        }
1196
1197      tree->rates->n_tips_below[d->num] = n1+n2;
1198
1199      (*logn) += Factln(n1+n2-2) - Factln(n1-1) - Factln(n2-1);
1200    }
1201}
1202
1203//////////////////////////////////////////////////////////////
1204//////////////////////////////////////////////////////////////
1205
1206
1207void TIMES_Update_Curr_Slice(t_tree *tree)
1208{
1209  int i,j;
1210
1211  For(i,2*tree->n_otu-1)
1212    {
1213      For(j,tree->rates->n_time_slices)
1214        {
1215          if(!(tree->rates->nd_t[i] > tree->rates->time_slice_lims[j])) break;
1216        }
1217      tree->rates->curr_slice[i] = j;
1218
1219      /* PhyML_Printf("\n. Node %3d [%12f] is in slice %3d.",i,tree->rates->nd_t[i],j); */
1220    }
1221}
1222
1223//////////////////////////////////////////////////////////////
1224//////////////////////////////////////////////////////////////
1225
1226
1227phydbl TIMES_Lk_Uniform_Core(t_tree *tree)
1228{ 
1229  phydbl logn;
1230
1231  logn = TIMES_Log_Number_Of_Ranked_Labelled_Histories(tree->n_root,YES,tree);
1232
1233  tree->rates->c_lnL_times = 0.0;
1234  TIMES_Lk_Uniform_Post(tree->n_root,tree->n_root->v[2],tree);
1235  TIMES_Lk_Uniform_Post(tree->n_root,tree->n_root->v[1],tree);
1236
1237  /* printf("\n. ^ %f %f %f", */
1238  /*     (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.), */
1239  /*     LOG(tree->rates->time_slice_lims[tree->rates->curr_slice[tree->n_root->num]] - */
1240  /*         tree->rates->nd_t[tree->n_root->num]), */
1241  /*     (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.) * */
1242  /*     LOG(tree->rates->time_slice_lims[tree->rates->curr_slice[tree->n_root->num]] - */
1243  /*         tree->rates->nd_t[tree->n_root->num])); */
1244
1245  tree->rates->c_lnL_times +=
1246    Factln(tree->rates->n_tips_below[tree->n_root->num]-2.) -
1247    (phydbl)(tree->rates->n_tips_below[tree->n_root->num]-2.) *
1248    LOG(tree->rates->time_slice_lims[tree->rates->curr_slice[tree->n_root->num]] -
1249        tree->rates->nd_t[tree->n_root->num]);
1250 
1251  tree->rates->c_lnL_times -= logn;
1252 
1253  return(tree->rates->c_lnL_times);
1254}
1255
1256//////////////////////////////////////////////////////////////
1257//////////////////////////////////////////////////////////////
1258
1259
1260void TIMES_Get_Number_Of_Time_Slices(t_tree *tree)
1261{
1262  int i;
1263
1264  tree->rates->n_time_slices=0;
1265  TIMES_Get_Number_Of_Time_Slices_Post(tree->n_root,tree->n_root->v[2],tree);
1266  TIMES_Get_Number_Of_Time_Slices_Post(tree->n_root,tree->n_root->v[1],tree);
1267  Qksort(tree->rates->time_slice_lims,NULL,0,tree->rates->n_time_slices-1);
1268
1269  if(tree->rates->n_time_slices > 1)
1270    {
1271      PhyML_Printf("\n");
1272      PhyML_Printf("\n. Sequences were collected at %d different time points.",tree->rates->n_time_slices);
1273      For(i,tree->rates->n_time_slices) printf("\n+ [%3d] time point @ %12f ",i+1,tree->rates->time_slice_lims[i]);
1274    }
1275}
1276
1277//////////////////////////////////////////////////////////////
1278//////////////////////////////////////////////////////////////
1279
1280
1281void TIMES_Get_Number_Of_Time_Slices_Post(t_node *a, t_node *d, t_tree *tree)
1282{
1283  int i;
1284
1285  if(d->tax == YES)
1286    {
1287      For(i,tree->rates->n_time_slices) 
1288        if(Are_Equal(tree->rates->t_floor[d->num],tree->rates->time_slice_lims[i],1.E-6)) break;
1289
1290      if(i == tree->rates->n_time_slices) 
1291        {
1292          tree->rates->time_slice_lims[i] = tree->rates->t_floor[d->num];
1293          tree->rates->n_time_slices++;
1294        }
1295      return;
1296    }
1297  else
1298    {
1299      For(i,3)
1300        if(d->v[i] != a && d->b[i] != tree->e_root)
1301          TIMES_Get_Number_Of_Time_Slices_Post(d,d->v[i],tree);
1302    }
1303}
1304
1305//////////////////////////////////////////////////////////////
1306//////////////////////////////////////////////////////////////
1307
1308
1309void TIMES_Get_N_Slice_Spans(t_tree *tree)
1310{
1311  int i,j;
1312
1313  For(i,2*tree->n_otu-2)
1314    {
1315      if(tree->a_nodes[i]->tax == NO)
1316        {
1317          For(j,tree->rates->n_time_slices)
1318            {
1319              if(Are_Equal(tree->rates->t_floor[i],tree->rates->time_slice_lims[j],1.E-6))
1320                {
1321                  tree->rates->n_time_slice_spans[i] = j+1;
1322                  /* PhyML_Printf("\n. Node %3d spans %3d slices [%12f].", */
1323                  /*           i+1, */
1324                  /*           tree->rates->n_slice_spans[i], */
1325                  /*           tree->rates->t_floor[i]); */
1326                  break;
1327                }
1328            }
1329        }
1330    }
1331}
1332
1333//////////////////////////////////////////////////////////////
1334//////////////////////////////////////////////////////////////
1335
1336
1337void TIMES_Lk_Uniform_Post(t_node *a, t_node *d, t_tree *tree)
1338{
1339  if(d->tax == YES) return;
1340  else
1341    {
1342      int i;
1343
1344      For(i,3)
1345        {
1346          if(d->v[i] != a && d->b[i] != tree->e_root)
1347            {
1348              TIMES_Lk_Uniform_Post(d,d->v[i],tree);
1349            }
1350        }
1351     
1352      if(tree->rates->curr_slice[a->num] != tree->rates->curr_slice[d->num])
1353        {
1354          tree->rates->c_lnL_times += 
1355            Factln(tree->rates->n_tips_below[d->num]-1.) - 
1356            (phydbl)(tree->rates->n_tips_below[d->num]-1.) *
1357            LOG(tree->rates->time_slice_lims[tree->rates->curr_slice[d->num]] -
1358                tree->rates->nd_t[d->num]);
1359        }
1360    }
1361}
1362
1363//////////////////////////////////////////////////////////////
1364//////////////////////////////////////////////////////////////
1365
1366
1367/* Set the root position so that most of the taxa in the outgroup
1368   correspond to the most ancient time point.
1369*/
1370void TIMES_Set_Root_Given_Tip_Dates(t_tree *tree)
1371{
1372  int i,j;
1373  t_node *left,*rght;
1374  int n_left_in, n_left_out;
1375  int n_rght_in, n_rght_out;
1376  t_edge *b,*best;
1377  phydbl eps,score,max_score;
1378 
1379  Free_Bip(tree);
1380  Alloc_Bip(tree);
1381  Get_Bip(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
1382 
1383  left = rght = NULL;
1384  b = best = NULL;
1385  n_left_in = n_left_out = -1;
1386  n_rght_in = n_rght_out = -1;
1387  eps = 1.E-6;
1388  score = max_score = -1.;
1389
1390  For(i,2*tree->n_otu-3)
1391    {
1392      left = tree->a_edges[i]->left;
1393      rght = tree->a_edges[i]->rght;
1394      b    = tree->a_edges[i];
1395
1396      n_left_in = 0;
1397      For(j,left->bip_size[b->l_r]) 
1398        if(FABS(tree->rates->nd_t[left->bip_node[b->l_r][j]->num] - tree->rates->time_slice_lims[0]) < eps)
1399          n_left_in++;
1400     
1401      n_left_out = left->bip_size[b->l_r]-n_left_in;
1402     
1403      n_rght_in = 0;
1404      For(j,rght->bip_size[b->r_l]) 
1405        if(FABS(tree->rates->nd_t[rght->bip_node[b->r_l][j]->num] - tree->rates->time_slice_lims[0]) < eps)
1406          n_rght_in++;
1407
1408      n_rght_out = rght->bip_size[b->r_l]-n_rght_in;
1409
1410
1411      /* score = POW((phydbl)(n_left_in)/(phydbl)(n_left_in+n_left_out)- */
1412      /*                  (phydbl)(n_rght_in)/(phydbl)(n_rght_in+n_rght_out),2); */
1413      /* score = (phydbl)(n_left_in * n_rght_out + eps)/(n_left_out * n_rght_in + eps); */
1414      /* score = (phydbl)(n_left_in * n_rght_out + eps); */
1415      score = FABS((phydbl)((n_left_in+1.) * (n_rght_out+1.)) - (phydbl)((n_left_out+1.) * (n_rght_in+1.)));
1416     
1417      if(score > max_score)
1418        {
1419          max_score = score;
1420          best = b;
1421        }
1422    }
1423 
1424  Add_Root(best,tree);
1425}
1426
1427//////////////////////////////////////////////////////////////
1428//////////////////////////////////////////////////////////////
1429
1430
1431void Get_Survival_Duration(t_tree *tree)
1432{
1433  Get_Survival_Duration_Post(tree->n_root,tree->n_root->v[2],tree);
1434  Get_Survival_Duration_Post(tree->n_root,tree->n_root->v[1],tree);
1435}
1436
1437//////////////////////////////////////////////////////////////
1438//////////////////////////////////////////////////////////////
1439
1440
1441void Get_Survival_Duration_Post(t_node *a, t_node *d, t_tree *tree)
1442{
1443  if(d->tax)
1444    {
1445      tree->rates->survival_dur[d->num] = tree->rates->nd_t[d->num];
1446      return;
1447    }
1448  else
1449    {
1450      int i;
1451      t_node *v1, *v2;
1452
1453      For(i,3)
1454        if(d->v[i] != a && d->b[i] != tree->e_root)
1455          Get_Survival_Duration_Post(d,d->v[i],tree);
1456     
1457      v1 = v2 = NULL;
1458      For(i,3)
1459        {
1460          if(d->v[i] != a && d->b[i] != tree->e_root)
1461            {
1462              if(!v1) v1 = d->v[i];
1463              else    v2 = d->v[i];
1464            }
1465        }
1466
1467      tree->rates->survival_dur[d->num] = MAX(tree->rates->survival_dur[v1->num],
1468                                              tree->rates->survival_dur[v2->num]);
1469    }
1470}
1471
1472
1473//////////////////////////////////////////////////////////////
1474//////////////////////////////////////////////////////////////
1475
1476/* Update the ranking of node heights. Use bubble sort algorithm */
1477
1478void TIMES_Update_Node_Ordering(t_tree *tree)
1479{
1480  int buff;
1481  int i;
1482  phydbl *t;
1483  int swap = NO;
1484
1485  t = tree->rates->nd_t;
1486
1487  do
1488    {
1489      swap = NO;
1490      For(i,2*tree->n_otu-2)
1491        {
1492          if(t[tree->rates->t_ranked[i]] > t[tree->rates->t_ranked[i+1]]) // Sort in ascending order
1493            {
1494              swap = YES;
1495              buff                       = tree->rates->t_ranked[i];
1496              tree->rates->t_ranked[i]   = tree->rates->t_ranked[i+1];
1497              tree->rates->t_ranked[i+1] = buff;
1498            }       
1499        }
1500    }
1501  while(swap == YES);
1502
1503  /* For(i,2*tree->n_otu-1) */
1504  /*   { */
1505  /*     printf("\n. ..... %f",t[tree->rates->t_ranked[i]]); */
1506  /*   } */
1507}
1508
1509//////////////////////////////////////////////////////////////
1510//////////////////////////////////////////////////////////////
1511
1512void TIMES_Label_Edges_With_Calibration_Intervals(t_tree *tree)
1513{
1514  char *s;
1515  int i;
1516
1517  s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1518 
1519  tree->write_labels = YES;
1520
1521  For(i,2*tree->n_otu-2)
1522    {
1523      if(tree->a_nodes[i]->tax == NO)
1524        {
1525          if(tree->rates->t_has_prior[i] == YES && tree->a_nodes[i]->b[0] != tree->e_root)
1526            {
1527              tree->a_nodes[i]->b[0]->n_labels = 1;
1528              Make_New_Edge_Label(tree->a_nodes[i]->b[0]);
1529              sprintf(s,"'>%f<%f'",FABS(tree->rates->t_prior_max[i])/100.,FABS(tree->rates->t_prior_min[i])/100.);
1530              tree->a_nodes[i]->b[0]->labels[0] = (char *)mCalloc(strlen(s)+1,sizeof(char));
1531              strcpy(tree->a_nodes[i]->b[0]->labels[0],s);
1532            }
1533        }
1534    }
1535 
1536  Free(s);
1537
1538}
1539
1540//////////////////////////////////////////////////////////////
1541//////////////////////////////////////////////////////////////
1542
1543void TIMES_Set_Calibration(t_tree *tree)
1544{
1545  t_cal *cal;
1546  int i;
1547
1548  For(i,2*tree->n_otu-1)
1549    {
1550      tree->rates->t_has_prior[i] = NO;
1551      tree->rates->t_prior_min[i] = BIG;
1552      tree->rates->t_prior_max[i] = BIG; 
1553   }
1554
1555  cal = tree->rates->calib;
1556  while(cal)
1557    {
1558      /* if(cal->is_active == YES) */
1559      /*   { */
1560          /* tree->rates->t_has_prior[cal->node_num] = YES; */
1561          /* tree->rates->t_prior_min[cal->node_num] = cal->lower; */
1562          /* tree->rates->t_prior_max[cal->node_num] = cal->upper;           */
1563        /* } */
1564      cal = cal->next;
1565    }
1566
1567  TIMES_Set_All_Node_Priors(tree);
1568}
1569
1570
1571
1572//////////////////////////////////////////////////////////////
1573//////////////////////////////////////////////////////////////
1574
1575
1576void TIMES_Record_Prior_Times(t_tree *tree)
1577{
1578  int i;
1579  For(i,2*tree->n_otu-1) 
1580    {
1581      tree->rates->t_prior_min_buff[i] = tree->rates->t_prior_min[i];
1582      tree->rates->t_prior_max_buff[i] = tree->rates->t_prior_max[i];
1583    }
1584}
1585
1586//////////////////////////////////////////////////////////////
1587//////////////////////////////////////////////////////////////
1588
1589
1590void TIMES_Reset_Prior_Times(t_tree *tree)
1591{
1592  int i;
1593  For(i,2*tree->n_otu-1) 
1594    {
1595      tree->rates->t_prior_min[i] = tree->rates->t_prior_min_buff[i];
1596      tree->rates->t_prior_max[i] = tree->rates->t_prior_max_buff[i];
1597     }
1598}
1599
1600//////////////////////////////////////////////////////////////
1601//////////////////////////////////////////////////////////////
1602//////////////////////////////////////////////////////////////
1603//////////////////////////////////////////////////////////////
1604//////////////////////////////////////////////////////////////
1605//////////////////////////////////////////////////////////////
1606//////////////////////////////////////////////////////////////
1607//////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.