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

/*

PhyML:  a program that  computes maximum likelihood phyLOGenies from
DNA or AA homoLOGous sequences.

Copyright (C) Stephane Guindon. Oct 2003 onward.

All parts of the source except where indicated are distributed under
the GNU public licence. See http://www.opensource.org for details.

*/

/* Routines for molecular clock trees and molecular dating */


#include "rates.h"

#ifdef RWRAPPER
#include <R.h>
#endif



//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Lk_Rates(t_tree *tree)
{

  tree->rates->c_lnL_rates  = .0;

  RATES_Lk_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
  RATES_Lk_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);

  if(isnan(tree->rates->c_lnL_rates) || isinf(tree->rates->c_lnL_rates))
    {
      PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  return tree->rates->c_lnL_rates;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Lk_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{
  int i;
  phydbl log_dens,mu_a,mu_d,r_a,r_d,dt_a,dt_d;
  int n_a,n_d;

  log_dens = -1.;

  if(d->anc != a)
    {
      PhyML_Printf("\n. d=%d d->anc=%d a=%d root=%d",d->num,d->anc->num,a->num,tree->n_root->num);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  dt_a = -1.;
  if(a != tree->n_root) dt_a = tree->rates->nd_t[a->num] - tree->rates->nd_t[a->anc->num];

  mu_a = tree->rates->br_r[a->num];
  r_a  = tree->rates->nd_r[a->num];
  n_a  = tree->rates->n_jps[a->num];
  
  dt_d = FABS(tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num]);
  mu_d = tree->rates->br_r[d->num];
  r_d  = tree->rates->nd_r[d->num];
  n_d  = tree->rates->n_jps[d->num];

  log_dens = RATES_Lk_Rates_Core(mu_a,mu_d,r_a,r_d,n_a,n_d,dt_a,dt_d,tree);
  tree->rates->c_lnL_rates += log_dens;

  /* printf("\n. a=%3d d=%3d r(a)=%10f r(d)=%10f %f",a->num,d->num,mu_a,mu_d,log_dens); */

  if(isnan(tree->rates->c_lnL_rates))
    {
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      MCMC_Print_Param(tree->mcmc,tree);
      Exit("\n");
    }

  tree->rates->triplet[a->num] += log_dens;


  if(d->tax) return;
  else
    {
      For(i,3)
        {
          if((d->v[i] != a) && (d->b[i] != tree->e_root))
            {
              RATES_Lk_Rates_Pre(d,d->v[i],d->b[i],tree);
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Lk_Change_One_Rate(t_node *d, phydbl new_rate, t_tree *tree)
{
  tree->rates->br_r[d->num] = new_rate;
  RATES_Update_Triplet(d,tree);
  RATES_Update_Triplet(d->anc,tree);
  return(tree->rates->c_lnL_rates);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Lk_Change_One_Time(t_node *n, phydbl new_t, t_tree *tree)
{  
  if(n == tree->n_root)
    {
      PhyML_Printf("\n. Moving the time of the root t_node is not permitted.");
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }
  else
    {
      int i;
      
      tree->rates->nd_t[n->num] = new_t;

      RATES_Update_Triplet(n,tree);
      
      For(i,3)
        {
          if(n->b[i] != tree->e_root) RATES_Update_Triplet(n->v[i],tree);
          else RATES_Update_Triplet(tree->n_root,tree);
        }
    }
  return(tree->rates->c_lnL_rates);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Triplet(t_node *n, t_tree *tree)
{
  phydbl curr_triplet,new_triplet;
  phydbl dt0,dt1,dt2;
  phydbl mu1_mu0,mu2_mu0;
  phydbl mu0,mu1,mu2;
  phydbl r0,r1,r2;
  int n0,n1,n2;
  int i;
  t_node *v1,*v2;

  if(n->tax) return;

  curr_triplet = tree->rates->triplet[n->num];

  dt0 = dt1 = dt2 = -100.0;
  r0 = r1 = r2 = 0.0;

  if(n == tree->n_root)
    {
      phydbl log_dens;
      
      log_dens = 0.0;

      dt0 = tree->rates->nd_t[tree->n_root->v[2]->num] - tree->rates->nd_t[tree->n_root->num];
      dt1 = tree->rates->nd_t[tree->n_root->v[1]->num] - tree->rates->nd_t[tree->n_root->num];
      
      mu0 = tree->rates->br_r[tree->n_root->v[2]->num];
      mu1 = tree->rates->br_r[tree->n_root->v[1]->num];

      r0 = tree->rates->nd_r[tree->n_root->v[2]->num];
      r1 = tree->rates->nd_r[tree->n_root->v[1]->num];
      
      n0  = tree->rates->n_jps[tree->n_root->v[2]->num];
      n1  = tree->rates->n_jps[tree->n_root->v[1]->num];

      switch(tree->rates->model)
        {
        case COMPOUND_COR : case COMPOUND_NOCOR : 
          {
            log_dens  = RATES_Dmu(mu0,n0,dt0,tree->rates->nu,1./tree->rates->nu,tree->rates->lexp,0,1);
            log_dens *= RATES_Dmu(mu1,n1,dt1,tree->rates->nu,1./tree->rates->nu,tree->rates->lexp,0,1);
            log_dens  = LOG(log_dens);
            break;
          }
        case EXPONENTIAL : 
          {
            log_dens = Dexp(mu0,tree->rates->lexp) * Dexp(mu1,tree->rates->lexp);
            log_dens = LOG(log_dens);
            break;
          }
        case GAMMA :
          {
            log_dens = Dgamma(mu0,tree->rates->nu,1./tree->rates->nu) * Dgamma(mu1,tree->rates->nu,1./tree->rates->nu);
            log_dens = LOG(log_dens);
            break;
          }
        case THORNE :
          {
            int err;
            phydbl mean0,sd0;
            phydbl mean1,sd1;

            
            sd0 = SQRT(tree->rates->nu*dt0);
            mean0 = 1.0;

            sd1 = SQRT(tree->rates->nu*dt1);
            mean1 = 1.0;

            log_dens = 
              Log_Dnorm_Trunc(mu0,mean0,sd0,tree->rates->min_rate,tree->rates->max_rate,&err) + 
              Log_Dnorm_Trunc(mu1,mean1,sd1,tree->rates->min_rate,tree->rates->max_rate,&err); 

            break;
          }
        case GUINDON :
          {
            Exit("\n. Not implemented yet.\n");
            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
            break;
          }
        default :
          {
            Exit("\n. Model not implemented yet.\n");
            break;
          }
        }
      new_triplet = log_dens;

      if(isnan(log_dens) || isinf(FABS(log_dens)))
        {
          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
          MCMC_Print_Param(tree->mcmc,tree);
          Exit("\n");
        }
    }
  else
    {
      mu0 = mu1 = mu2 = -1.;
      n0 = n1 = n2 = -1;

      mu0 = tree->rates->br_r[n->num];
      dt0 = FABS(tree->rates->nd_t[n->num] - tree->rates->nd_t[n->anc->num]);
      n0  = tree->rates->n_jps[n->num];
      r0 = tree->rates->nd_r[n->num];

      v1 = v2 = NULL;
      For(i,3)
        {
          if((n->v[i] != n->anc) && (n->b[i] != tree->e_root))
            {
              if(!v1)
                {
                  v1  = n->v[i]; 
                  mu1 = tree->rates->br_r[v1->num];
                  dt1 = FABS(tree->rates->nd_t[v1->num] - tree->rates->nd_t[n->num]);
                  n1  = tree->rates->n_jps[v1->num];
                  r1  = tree->rates->nd_r[v1->num];
                }
              else
                {
                  v2  = n->v[i]; 
                  mu2 = tree->rates->br_r[v2->num];
                  dt2 = FABS(tree->rates->nd_t[v2->num] - tree->rates->nd_t[n->num]);
                  n2  = tree->rates->n_jps[v2->num];
                  r2  = tree->rates->nd_r[v2->num];
                }
            }
        }
 
      mu1_mu0 = RATES_Lk_Rates_Core(mu0,mu1,r0,r1,n0,n1,dt0,dt1,tree);
      mu2_mu0 = RATES_Lk_Rates_Core(mu0,mu2,r0,r1,n0,n2,dt0,dt2,tree);
      
      new_triplet = mu1_mu0 + mu2_mu0;
    }

  tree->rates->c_lnL_rates = tree->rates->c_lnL_rates + new_triplet - curr_triplet;
  tree->rates->triplet[n->num] = new_triplet;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* Returns LOG(f(br_r_rght;br_r_left)) */
phydbl RATES_Lk_Rates_Core(phydbl br_r_a, phydbl br_r_d, phydbl nd_r_a, phydbl nd_r_d, int n_a, int n_d, phydbl dt_a, phydbl dt_d, t_tree *tree)
{
  phydbl log_dens;
  phydbl mean,sd;
  phydbl min_r, max_r;
  phydbl cr,logcr;
  
  cr        = tree->rates->clock_r;
  logcr     = LOG(cr);
  log_dens  = UNLIKELY;
  mean = sd = -1.;
  min_r     = tree->rates->min_rate;
  max_r     = tree->rates->max_rate;

  dt_d = MAX(0.5,dt_d); // We give only one decimal when printing out node heights. It is therefore a fair approximation

  switch(tree->rates->model)
    {
    case THORNE :
      {
        int err;        

        if(tree->rates->model_log_rates == YES)
          {

            br_r_d += logcr;
            br_r_a += logcr;
            min_r  += logcr;
            max_r  += logcr;

            sd   = SQRT(tree->rates->nu*dt_d);
            mean = br_r_a - .5*sd*sd;

            log_dens = Log_Dnorm_Trunc(br_r_d,mean,sd,min_r,max_r,&err);

            if(err)
              {
                PhyML_Printf("\n. Run: %d",tree->mcmc->run);
                PhyML_Printf("\n. br_r_d = %f br_r_a = %f dt_d = %f logcr = %f",br_r_d,br_r_a,dt_d,logcr);
                PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
                Exit("\n");
              }
          }
        else
          {
            PhyML_Printf("\n. Not implemented yet.");
            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
            Exit("\n");
          }

        /* int err; */
        /* phydbl cr; */
        /* phydbl min_r, max_r; */

        /* cr = tree->rates->clock_r; */
        /* min_r = tree->rates->min_rate * cr; */
        /* max_r = tree->rates->max_rate * cr; */

        /* /\* !!!!!!!!!!!!1 *\/ */
        /* br_r_d *= cr; */
        /* br_r_a *= cr; */

        /* err = NO; */
        
        /* /\* if(dt_d < 5.0) dt_d = 5.0; *\/ */

        /* sd = SQRT(dt_d*tree->rates->nu); */

        /* /\* sd   = SQRT(dt_d*EXP(tree->rates->nu)); *\/ */
        /* /\* mean = LOG(br_r_a) - .5*sd*sd; *\/ */

        /* if(tree->rates->model_log_rates == YES) */
        /*   { */
        /*     mean = br_r_a - .5*sd*sd; */
        /*   } */
        /* else */
        /*   { */
        /*     mean = br_r_a; */
        /*   } */

        /* log_dens = Log_Dnorm_Trunc(br_r_d,mean,sd,min_r,max_r,&err); */

        /* if(err || isnan(log_dens) || isinf(log_dens)) */
        /*   { */
        /*     PhyML_Printf("\n. br_r_a=%f br_r_d=%f dt_d=%f nu=%G",br_r_a,br_r_d,dt_d,tree->rates->nu); */
        /*     PhyML_Printf("\n. Run: %d",tree->mcmc->run); */
        /*     PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
        /*     Exit("\n"); */
        /*   } */

        break;

      }
    case GUINDON :
      {
        int err;        

        min_r = tree->rates->min_rate;
        max_r = tree->rates->max_rate;

        br_r_d += logcr;
        br_r_a += logcr;
        min_r  += logcr;
        max_r  += logcr;

        sd   = SQRT(tree->rates->nu*dt_d);
        mean = br_r_a - .5*sd*sd;

        log_dens = Log_Dnorm_Trunc(br_r_d,mean,sd,min_r,max_r,&err);

        if(err)
          {
            PhyML_Printf("\n. Run: %d",tree->mcmc->run);
            PhyML_Printf("\n. br_r_d=%f mean=%f sd=%f min_r=%f max_r=%f dt_d=%f",br_r_d,mean,sd,min_r,max_r,dt_d);
            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
            Exit("\n");
          }
        break;
      }
    case GAMMA :
      {
        if(tree->rates->model_log_rates == YES)
          log_dens = Dgamma(EXP(br_r_d),1./tree->rates->nu,tree->rates->nu);
        else
          log_dens = Dgamma(br_r_d,1./tree->rates->nu,tree->rates->nu);

        log_dens = LOG(log_dens);
        break;
      }
    case STRICTCLOCK :
      {
        log_dens = 0.0;
        break;
      }

    default : 
      {
        PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
        Warn_And_Exit("");
      }
    }

  if(isnan(log_dens) || isinf(FABS(log_dens)))
    {
      PhyML_Printf("\n. Run=%4d br_r_d=%f br_r_a=%f dt_d=%f dt_a=%f nu=%f log_dens=%G sd=%f mean=%f\n",
                   tree->mcmc->run,
                   br_r_d,br_r_a,dt_d,dt_a,tree->rates->nu,log_dens,
                   sd,mean);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  return log_dens;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Compound_Core(phydbl mu1, phydbl mu2, int n1, int n2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx)
{
  if((n1 > -1) && (n2 > -1))
    {
      return RATES_Compound_Core_Joint(mu1,mu2,n1,n2,dt1,dt2,alpha,beta,lexp,eps,approx);
    }
  else
    {
      return RATES_Compound_Core_Marginal(mu1,mu2,dt1,dt2,alpha,beta,lexp,eps,approx);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Compound_Core_Marginal(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx)
{
  phydbl p0,p1,v0,v1,v2;
  phydbl dmu1;
  int    equ;
  phydbl dens;
  
  v0 = v1 = v2 = 0.0;
  
  /* Probability of 0 and 1 jumps */
  p0   = Dpois(0,lexp*(dt2+dt1));       
  p1   = Dpois(1,lexp*(dt2+dt1));
  
  dmu1 = RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,0,0);
  
  /* Are the two rates equal ? */
  equ = 0;
  if(FABS(mu1-mu2) < eps) equ = 1;
  
  /* No jump */
  if(equ)
    {
      v0 = 1.0*Dgamma(mu1,alpha,beta)/dmu1;
      /*       Rprintf("\n. mu1=%f mu2=%f",mu1,mu2); */
    }
  else
    {
      v0 = 1.E-100;
    }

  /* One jump */
  v1 = RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(mu1,mu2,dt1,dt2,alpha,beta);
  v1 /= dmu1;
    
  /* Two jumps and more (approximation) */
  if(approx == 1)
    {
      v2 =
        RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,0,0)*RATES_Dmu(mu2,-1,dt2,alpha,beta,lexp,0,0) -
        Dpois(0,lexp*dt1) * Dpois(0,lexp*dt2) *
        Dgamma(mu1,alpha,beta) * Dgamma(mu2,alpha,beta);
    }
  else
    {
      v2 = 
        RATES_Dmu_One(mu1,dt1,alpha,beta,lexp) * 
        RATES_Dmu_One(mu2,dt2,alpha,beta,lexp);

      v2 += Dpois(0,lexp*dt1)*Dgamma(mu1,alpha,beta)*RATES_Dmu(mu2,-1,dt2,alpha,beta,lexp,1,0);
      v2 += Dpois(0,lexp*dt2)*Dgamma(mu2,alpha,beta)*RATES_Dmu(mu1,-1,dt1,alpha,beta,lexp,1,0);

    }
/*   PhyML_Printf("\n. %f %f %f %f %f ",mu1,mu2,dt1,dt2,v2); */
  v2 /= dmu1;

  dens = p0*v0 + p1*v1 + v2;
/*   dens = p1*v1 + v2; */
  /*       dens = p1*v1 + v2; */
  /*   dens = v0; */
/*   dens *= dmu1; */
  
  if(dens < SMALL)
    {
      PhyML_Printf("\n. dens=%12G mu1=%12G mu2=%12G dt1=%12G dt2=%12G lexp=%12G alpha=%f v0=%f v1=%f v2=%f p0=%f p1=%f p2=%f",
             dens,
             mu1,mu2,dt1,dt2,
             lexp,
             alpha,
             v0,v1,v2,
             p0,p1,1.-p0-p1);
    }
  
  return dens;

}

/**********************************************************/

phydbl RATES_Compound_Core_Joint(phydbl mu1, phydbl mu2, int n1, int n2, phydbl dt1, phydbl dt2, 
                                 phydbl alpha, phydbl beta, phydbl lexp, phydbl eps, int approx)
{
  phydbl density;
  phydbl dmu1;
  

  if(n1 < 0 || n2 < 0)
    {
      PhyML_Printf("\n. n1=%d n2=%d",n1,n2);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  dmu1 = RATES_Dmu(mu1,n1,dt1,alpha,beta,lexp,0,0);
  
  if((n1 < 0) || (n2 < 0))
    {
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  if((n1 == 0) && (n2 == 0))
    {
      if(FABS(mu1-mu2) < eps) { density = Dgamma(mu1,alpha,beta); }
      else                    { density = 1.E-70; }
    }
  else if((n1 == 0) && (n2 == 1))
    {
      density = 
        Dgamma(mu1,alpha,beta) *
        RATES_Dmu1_Given_V_And_N(mu2,mu1,1,dt2,alpha,beta);
    }
  else if((n1 == 1) && (n2 == 0))
    {
      density = 
        Dgamma(mu2,alpha,beta) *
        RATES_Dmu1_Given_V_And_N(mu1,mu2,1,dt1,alpha,beta);
    }
  else /* independent */
    {
      density = 
        RATES_Dmu(mu1,n1,dt1,alpha,beta,lexp,0,0) * 
        RATES_Dmu(mu2,n2,dt2,alpha,beta,lexp,0,0);
    }

  density /= dmu1;

  density *= Dpois(n2,dt2*lexp);
  
  if(density < 1.E-70) density = 1.E-70;

/*   PhyML_Printf("\n. density = %15G mu1=%3.4f mu2=%3.4f dt1=%3.4f dt2=%3.4f n1=%2d n2=%2d",density,mu1,mu2,dt1,dt2,n1,n2); */
  return density;
}

/**********************************************************/

void RATES_Print_Triplets(t_tree *tree)
{
  int i;
  For(i,2*tree->n_otu-1) PhyML_Printf("\n. Node %3d t=%f",i,tree->rates->triplet[i]);
}


/**********************************************************/

void RATES_Print_Rates(t_tree *tree)
{
  RATES_Print_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
  RATES_Print_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Print_Rates_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{  
  if((d == tree->n_root->v[2] && d->tax) || (d == tree->n_root->v[1] && d->tax))
    PhyML_Printf("\n. a=%3d ++d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
           a->num,d->num,
           tree->rates->br_r[d->num],
           tree->rates->nd_t[a->num],tree->rates->nd_t[d->num],
           tree->rates->ml_l[d->num],
           tree->rates->cur_l[d->num],
           (tree->rates->nd_t[d->num]-tree->rates->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
  
  else if((d == tree->n_root->v[2]) || (d == tree->n_root->v[1]))
    PhyML_Printf("\n. a=%3d __d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
           a->num,d->num,
           tree->rates->br_r[d->num],
           tree->rates->nd_t[a->num],tree->rates->nd_t[d->num],
           tree->rates->ml_l[d->num],
           tree->rates->cur_l[d->num],
           (tree->rates->nd_t[d->num]-tree->rates->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
  else 
    PhyML_Printf("\n. a=%3d   d=%3d rate=%12f t_left=%12f t_rght=%12f ml=%12f l=%12f %12f",
           a->num,d->num,
           tree->rates->br_r[d->num],
           tree->rates->nd_t[a->num],tree->rates->nd_t[d->num],
           tree->rates->ml_l[d->num],
           tree->rates->cur_l[d->num],
           (tree->rates->nd_t[d->num]-tree->rates->nd_t[a->num])*tree->rates->clock_r*tree->rates->br_r[d->num]);
  
  if(d->tax) return;
  else
    {
      int i;

      For(i,3) 
        {
          if((d->v[i] != a) && (d->b[i] != tree->e_root))
            {
              RATES_Print_Rates_Pre(d,d->v[i],d->b[i],tree);
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Average_Rate(t_tree *tree)
{
  int i;
  phydbl sum;
  sum = 0.0;
  For(i,2*tree->n_otu-2) sum += tree->rates->br_r[i];
  return sum/(2*tree->n_otu-2);
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

phydbl RATES_Average_Substitution_Rate(t_tree *tree)
{
  phydbl sum_r,sum_dt;
  phydbl u,t,t_anc;
  int i;

  u = 0.0;
  sum_r  = 0.0;
  sum_dt = 0.0;

  /* For(i,2*tree->n_otu-3)  */
  /*   { */
  /*     t     = tree->rates->nd_t[i]; */
  /*     t_anc = tree->rates->nd_t[tree->a_nodes[i]->anc->num];       */
  /*     u = tree->a_edges[i]->l->v; */
  /*     if(tree->rates->model == GUINDON) u = tree->a_edges[i]->gamma_prior_mean; */
  /*     sum_r += u;       */
  /*     sum_dt += FABS(t-t_anc); */
  /*   } */

  For(i,2*tree->n_otu-3)  
    {
      if(tree->a_edges[i] != tree->e_root)
        {
          t     = tree->rates->nd_t[tree->a_edges[i]->left->num];
          t_anc = tree->rates->nd_t[tree->a_edges[i]->rght->num];
          u = tree->a_edges[i]->l->v;
          sum_r += u;
          sum_dt += FABS(t-t_anc);
        }
      
      sum_dt += FABS(tree->rates->nd_t[tree->n_root->v[2]->num] - tree->rates->nd_t[tree->n_root->num]);
      sum_dt += FABS(tree->rates->nd_t[tree->n_root->v[1]->num] - tree->rates->nd_t[tree->n_root->num]);

      u = tree->e_root->l->v;
      sum_r += u;        
    }
  return(sum_r / sum_dt);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Check_Mean_Rates_True(t_tree *tree)
{
  phydbl sum;
  int i;
  
  sum = 0.0;
  For(i,2*tree->n_otu-2) sum += tree->rates->true_r[i];
  return(sum/(phydbl)(2*tree->n_otu-2));
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


int RATES_Check_Node_Times(t_tree *tree)
{
  int err;
  err = NO;
  RATES_Check_Node_Times_Pre(tree->n_root,tree->n_root->v[2],&err,tree);
  RATES_Check_Node_Times_Pre(tree->n_root,tree->n_root->v[1],&err,tree);
  return err;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Check_Node_Times_Pre(t_node *a, t_node *d, int *err, t_tree *tree)
{
  if((tree->rates->nd_t[d->num] < tree->rates->nd_t[a->num]) || (FABS(tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num]) < 1.E-20))
    {
      PhyML_Printf("\n== a->t=%f d->t=%f",tree->rates->nd_t[a->num],tree->rates->nd_t[d->num]);
      PhyML_Printf("\n== a->t_prior_min=%f a->t_prior_max=%f",tree->rates->t_prior_min[a->num],tree->rates->t_prior_max[a->num]);
      PhyML_Printf("\n== d->t_prior_min=%f d->t_prior_max=%f",tree->rates->t_prior_min[d->num],tree->rates->t_prior_max[d->num]);
      *err = YES;
    }
  if(d->tax) return;
  else
    {
      int i;

      For(i,3)
        if((d->v[i] != a) && (d->b[i] != tree->e_root))   
          RATES_Check_Node_Times_Pre(d,d->v[i],err,tree);
    }
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////







void RATES_Bracket_N_Jumps(int *up, int *down, phydbl param)
{
  phydbl cdf,eps,a,b,c;
  int step;

  step = 10;
  eps = 1.E-10;
  cdf = 0.0;
  c = 1;
  
  while(cdf < 1.-eps)
    {
      c = (int)FLOOR(c * step);
      cdf = Ppois(c,param);      
    }
  
  a = 0.0;
  b = (c-a)/2.;
  step = 0;
  do
    {
      step++;
      cdf = Ppois(b,param);
      if(cdf < eps) a = b;
      else 
        {
          break;
        }
      b = (c-a)/2.;
    }
  while(step < 1000);
  
  if(step == 1000)
    {
      PhyML_Printf("\n. a=%f b=%f c=%f param=%f",a,b,c,param);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }
  *up = c;
  *down = a;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* 
   mu   : average rate of the time period dt
   dt   : time period to be considered
   a    : rate at a given time point is gamma distributed. a is the shape parameter
   b    : rate at a given time point is gamma distributed. b is the scale parameter
   lexp : the number of rate switches is Poisson distributed with parameter lexp * dt
*/ 
/* compute f(mu;dt,a,b,lexp), the probability density of mu. We need to integrate over the
   possible number of jumps (n) during the time interval dt */
phydbl RATES_Dmu(phydbl mu, int n_jumps, phydbl dt, phydbl a, phydbl b, phydbl lexp, int min_n, int jps_dens)
{
  if(n_jumps < 0) /* Marginal, i.e., the number of jumps is not fixed */
    {
      phydbl var,cumpoissprob,dens,mean,poissprob,ab2,gammadens,lexpdt;
      int n,up,down;
      
      var          = 0.0;
      cumpoissprob = 0.0;
      dens         = 0.0;
      n            = 0;
      mean         = a*b;
      ab2          = a*b*b;
      lexpdt       = lexp*dt;  
      
      RATES_Bracket_N_Jumps(&up,&down,lexpdt);
      For(n,MAX(down,min_n)-1) cumpoissprob += Dpois(n,lexpdt);
      
      for(n=MAX(down,min_n);n<up+1;n++)
        {
          poissprob    = Dpois(n,lexpdt); /* probability of having n jumps */      
          var          = (2./(n+2.))*ab2; /* var(mu|n) = var(mu|n=0) * 2 / (n+2) */
          gammadens    = Dgamma_Moments(mu,mean,var);
          dens         += poissprob * gammadens;
          cumpoissprob += poissprob;
          if(cumpoissprob > 1.-1.E-04) break;
        }
      
      if(dens < 1.E-70) dens = 1.E-70;

      return(dens);      
    }
  else /* Joint, i.e., return P(mu | dt, n_jumps) */
    {
      phydbl mean, var, density;


      mean = 1.0;
      var = (2./(n_jumps+2.))*a*b*b;

      if(jps_dens)
        density = Dgamma_Moments(mu,mean,var) * Dpois(n_jumps,dt*lexp);
      else
        density = Dgamma_Moments(mu,mean,var);
      
      if(density < 1.E-70) density = 1.E-70;

      return density;
    }
}
  
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

  
phydbl RATES_Dmu_One(phydbl mu, phydbl dt, phydbl a, phydbl b, phydbl lexp)
{
  phydbl var,cumpoissprob,dens,mean,poissprob,ab2,gammadens,lexpdt;
  int n,up,down;
  
  var          = 0.0;
  cumpoissprob = 0.0;
  dens         = 0.0;
  n            = 0;
  mean         = a*b;
  ab2          = a*b*b;
  lexpdt       = lexp*dt;  
  
  if(dt < 0.0)
    {
      PhyML_Printf("\n. dt=%f",dt);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }
  
  if(lexpdt < SMALL)
    {
      PhyML_Printf("\n. lexpdt=%G lexp=%G dt=%G",lexpdt,lexp,dt);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  if(mu < 1.E-10)
    {
      PhyML_Printf("\n. mu=%G",mu);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");      
    }
  
  RATES_Bracket_N_Jumps(&up,&down,lexpdt);

  For(n,MAX(1,down)-1) cumpoissprob += Dpois(n,lexpdt);

  for(n=MAX(1,down);n<up+1;n++) /* WARNING: we are considering that at least one jump occurs in the interval */
    {
      poissprob    = Dpois(n,lexpdt); /* probability of having n jumps */      
      var          = (n/((n+1)*(n+1)*(n+2)))*(POW(1-a*b,2) + 2/(n+1)*ab2) + 2*n*n*ab2/POW(n+1,3);      
      gammadens    = Dgamma_Moments(mu,mean,var);
      dens         += poissprob * gammadens;
      cumpoissprob += poissprob;
      if(cumpoissprob > 1.-1.E-06) break;
    }

  return(dens);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* Given the times of nodes a (ta) and d (td), the shape of the gamma distribution of instantaneous
   rates, the parameter of the exponential distribution of waiting times between rate jumps and the 
   instantaneous rate at t_node a, this function works out an expected number of (amino-acids or 
   nucleotide) substitutions per site.
*/
void RATES_Expect_Number_Subst(phydbl t_beg, phydbl t_end, phydbl r_beg,  int *n_jumps, phydbl *mean_r, phydbl *r_end, t_rate *rates, t_tree *tree)
{
  phydbl curr_r, curr_t, next_t;

  switch(rates->model)
    {
    case COMPOUND_COR:case COMPOUND_NOCOR:
      {
        /* Compound Poisson */
        if(rates->model == COMPOUND_COR)
          {
            curr_r  = r_beg;
            *mean_r = r_beg;
          }
        else
          {
            curr_r  = Rgamma(rates->nu,1./rates->nu);;
            *mean_r = curr_r;
          }

        curr_t = t_beg + Rexp(rates->lexp); /* Exponentially distributed waiting times */
        next_t = curr_t;
        
        *n_jumps = 0;
        while(curr_t < t_end)
          {
            curr_r = Rgamma(rates->nu,1./rates->nu); /* Gamma distributed random instantaneous rate */
            
            (*n_jumps)++;
            
            next_t = curr_t + Rexp(rates->lexp);
            
            if(next_t < t_end)
              {
                *mean_r = (1./(next_t - t_beg)) * (*mean_r * (curr_t - t_beg) + curr_r * (next_t - curr_t));
              }
            else
              {
                *mean_r = (1./(t_end - t_beg)) * (*mean_r * (curr_t - t_beg) + curr_r * (t_end - curr_t));
              }
            curr_t = next_t;
          }
        
        /*   PhyML_Printf("\n. [%3d %f %f]",*n_jumps,*mean_r,r_beg); */
        
        if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
        if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;

        *r_end = curr_r;
        break;
      }
    case EXPONENTIAL:
      {
        *mean_r = Rexp(rates->nu);

        if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
        if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;

        *r_end  = *mean_r;
        break;
      }
    case GAMMA:
      {
        *mean_r = Rgamma(rates->nu,1./rates->nu);

        if(*mean_r < rates->min_rate) *mean_r = rates->min_rate;
        if(*mean_r > rates->max_rate) *mean_r = rates->max_rate;

        *r_end  = *mean_r;
        break;
      }
    case THORNE:
      {
        phydbl sd,mean;
        int err;
        
        sd = SQRT(rates->nu*FABS(t_beg-t_end));
        mean = r_beg;

        *mean_r = Rnorm_Trunc(mean,sd,rates->min_rate,rates->max_rate,&err);

        if(err) PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run);
        *r_end  = *mean_r;
        break;
      }
    default:
      {
        PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
        Exit("\n. Model not implemented yet.\n");
        break;
      }
    }
}
  
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Mean_Rates_Pre(t_node *a, t_node *d, t_edge *b, phydbl r_a, t_tree *tree)
{
  phydbl a_t,d_t;
  phydbl mean_r;
  int n_jumps;
  phydbl r_d;

  a_t = tree->rates->nd_t[a->num];
  d_t = tree->rates->nd_t[d->num];
      
  RATES_Expect_Number_Subst(a_t,d_t,r_a,&n_jumps,&mean_r,&r_d,tree->rates,tree);
  
  tree->rates->br_r[d->num]   = mean_r;
  tree->rates->true_r[d->num] = mean_r;
  tree->rates->t_jps[d->num]  = n_jumps;


  /* Move to the next branches */
  if(d->tax) return;
  else
    {
      int i;
      
      For(i,3)
        {
          if((d->v[i] != a) && (d->b[i] != tree->e_root))
            {
              RATES_Get_Mean_Rates_Pre(d,d->v[i],d->b[i],r_d,tree);
            }
        }
    }

}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Random_Branch_Lengths(t_tree *tree)
{
  phydbl r0;
   
  r0 = 1.0;

  tree->rates->br_r[tree->n_root->num] = r0;

  RATES_Get_Mean_Rates_Pre(tree->n_root,tree->n_root->v[2],NULL,r0,tree);
  RATES_Get_Mean_Rates_Pre(tree->n_root,tree->n_root->v[1],NULL,r0,tree);

/*   RATES_Normalise_Rates(tree); */

  RATES_Update_Cur_Bl(tree);
  RATES_Initialize_True_Rates(tree);

  tree->n_root_pos = 
    tree->rates->cur_l[tree->n_root->v[2]->num] /
    (tree->rates->cur_l[tree->n_root->v[2]->num] + tree->rates->cur_l[tree->n_root->v[1]->num]);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Set_Node_Times(t_tree *tree)
{
  RATES_Set_Node_Times_Pre(tree->n_root,tree->n_root->v[2],tree);
  RATES_Set_Node_Times_Pre(tree->n_root,tree->n_root->v[1],tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Init_Triplets(t_tree *tree)
{
  int i;
  For(i,2*tree->n_otu-1) tree->rates->triplet[i] = 0.0;
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Set_Node_Times_Pre(t_node *a, t_node *d, t_tree *tree)
{
  if(d->tax) return;
  else
    {
      t_node *v1, *v2; /* the two sons of d */
      phydbl t_sup, t_inf;
      int i;

      v1 = v2 = NULL;
      For(i,3) if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
        {
          if(!v1) v1 = d->v[i]; 
          else    v2 = d->v[i];
        }
          
      t_inf = MIN(tree->rates->nd_t[v1->num],tree->rates->nd_t[v2->num]);
      t_sup = tree->rates->nd_t[a->num];

      if(t_sup > t_inf)
        {
          PhyML_Printf("\n. t_sup = %f t_inf = %f",t_sup,t_inf);
          PhyML_Printf("\n. Run = %d",tree->mcmc->run);
          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
          Warn_And_Exit("");
        }
            
      if(tree->rates->nd_t[d->num] > t_inf)      
        {
          tree->rates->nd_t[d->num] = t_inf;
          PhyML_Printf("\n. Time for t_node %d is larger than should be. Set it to %f",d->num,tree->rates->nd_t[d->num]);
        }
      else if(tree->rates->nd_t[d->num] < t_sup) 
        {
          tree->rates->nd_t[d->num] = t_sup;
          PhyML_Printf("\n. Time for t_node %d is lower than should be. Set it to %f",d->num,tree->rates->nd_t[d->num]);
        }

      For(i,3)
        {
          if((d->v[i] != a) && (d->b[i] != tree->e_root))
            {
              RATES_Set_Node_Times_Pre(d,d->v[i],tree);       
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////



phydbl RATES_Dmu1_And_Mu2_One_Jump_Two_Intervals(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, phydbl alpha, phydbl beta)
{
  phydbl dens;

  if(mu2 < 1.E-10)
    {
      PhyML_Printf("\n. mu2=%G",mu2);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  if(mu1 < 1.E-10)
    {
      PhyML_Printf("\n. mu2=%G",mu1);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  dens =
    ((dt1/(dt1+dt2)) * RATES_Dmu1_Given_V_And_N(mu1,mu2,1,dt1,alpha,beta) * Dgamma(mu2,alpha,beta)) +
    ((dt2/(dt1+dt2)) * RATES_Dmu1_Given_V_And_N(mu2,mu1,1,dt2,alpha,beta) * Dgamma(mu1,alpha,beta));

  return dens;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Dmu1_Given_V_And_N(phydbl mu1, phydbl v, int n, phydbl dt1, phydbl a, phydbl b)
{
  phydbl lbda,dens,h,u;
  phydbl mean,var;
  int n_points,i;
  phydbl ndb;
  phydbl end, beg;

  n_points = 20;

  end = MIN(mu1/v-0.01,0.99);
  beg = 0.01;
  
  dens = 0.0;
  
  if(end > beg)
    {
      mean = a*b;
      var = a*b*b*2./(n+1.);
      ndb = (phydbl)n/dt1;
      
      h = (end - beg) / (phydbl)n_points;
      
      lbda = beg;
      For(i,n_points-1) 
        {
          lbda += h;
          u = (mu1 - lbda*v)/(1.-lbda);
          
          if(u < 1.E-10)
            {
              PhyML_Printf("\n. u = %G",u);
              PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
              Warn_And_Exit("");
            }
          
          dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
        }
      dens *= 2.;
      
      lbda = beg;
      u = (mu1 - lbda*v)/(1.-lbda);
      if(u < 1.E-10)
        {
          PhyML_Printf("\n. mu1 = %f lambda = %f v = %f u = %G beg = %f end = %f",mu1,lbda,v,u,beg,end);
          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
          Warn_And_Exit("");
        }
      
      dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
      
      lbda = end;
      u = (mu1 - lbda*v)/(1.-lbda);
      if(u < 1.E-10)
        {
          PhyML_Printf("\n. u = %G",u);
          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
          Warn_And_Exit("");
        }
      
      dens += Dgamma_Moments(u,mean,var) / (1.-lbda) * ndb * POW(1.-lbda,n-1);
      
      dens *= (h/2.);
      dens *= dt1;
    }

  return(dens);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* Joint density of mu1 and a minimum number of jumps occuring in the interval dt1+dt2 given mu1. 
   1 jump occurs at the junction of the two intervals, which makes mu1 and mu2 independant */
phydbl RATES_Dmu2_And_Mu1_Given_Min_N(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n_min, phydbl a, phydbl b, phydbl lexp)
{
  phydbl density, lexpdt,cumpoiss,poiss;
  int i;
  int up,down;

  density = 0.0;
  lexpdt = lexp * (dt1+dt2);
  cumpoiss = 0.0;

  RATES_Bracket_N_Jumps(&up,&down,lexpdt);

  For(i,MAX(up,n_min)-1)
    {
      poiss = Dpois(i,lexpdt);
      cumpoiss = cumpoiss + poiss;
    }

  for(i=MAX(up,n_min);i<up;i++)
    {
/*       poiss = Dpois(i-1,lexpdt); /\* Complies with the no correlation model *\/ */
      poiss = Dpois(i,lexpdt);
      cumpoiss = cumpoiss + poiss;

      density = density + poiss * RATES_Dmu2_And_Mu1_Given_N(mu1,mu2,dt1,dt2,i-1,a,b,lexp);
/*       density = density + poiss * RATES_Dmu2_And_Mu1_Given_N_Full(mu1,mu2,dt1,dt2,i,a,b,lexp); */
      if(cumpoiss > 1.-1.E-6) break;
    }

  if(density < 0.0)
    {
      PhyML_Printf("\n. density=%f cmpoiss = %f i=%d n_min=%d mu1=%f mu2=%f dt1=%f dt2=%f",
             density,cumpoiss,i,n_min,mu1,mu2,dt1,dt2);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Warn_And_Exit("");
    }

  return(density);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* Joint density of mu1 and mu2 given the number of jumps (n) in the interval dt1+dt2, which are considered
   as independant. Hence, for n jumps occuring in dt1+dt2, the number of jumps occuring in dt1 and dt2 are
   n and 0, or n-1 and 1, or n-2 and 2 ... or 0 and n. This function sums over all these scenarios, with
   weights corresponding to the probability of each partitition.
 */

phydbl RATES_Dmu2_And_Mu1_Given_N(phydbl mu1, phydbl mu2, phydbl dt1, phydbl dt2, int n, phydbl a, phydbl b, phydbl lexp)
  {
    phydbl density,cumpoiss,poiss,abb,ab,LOGnf,LOGdt1,LOGdt2,nLOGdt;
    int i;

    density  = 0.0;
    abb      = a*b*b;
    ab       = a*b;
    cumpoiss = 0.0;
    poiss    = 0.0;
    LOGnf    = LnFact(n);
    LOGdt1   = LOG(dt1);
    LOGdt2   = LOG(dt2);
    nLOGdt   = n*LOG(dt1+dt2);

    For(i,n+1)
      {
        poiss = LOGnf - LnFact(i) - LnFact(n-i) + i*LOGdt1 + (n-i)*LOGdt2 - nLOGdt;
        poiss = EXP(poiss);
        cumpoiss = cumpoiss + poiss;

        if(mu2 < 1.E-10)
          {
            PhyML_Printf("\n. mu2=%f",mu2);
            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
            Warn_And_Exit("");
          }

        if(mu1 < 1.E-10)
          {
            PhyML_Printf("\n. mu1=%f",mu1);
            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
            Warn_And_Exit("");
          }


        density = density + poiss * Dgamma_Moments(mu2,ab,2./((n-i)+2.)*abb) * Dgamma_Moments(mu1,ab,2./(i+2.)*abb);
        if(cumpoiss > 1.-1.E-6) break;
      }

    return(density);
  }

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Initialize_True_Rates(t_tree *tree)
{
  int i;
  For(i,2*tree->n_otu-2) tree->rates->true_r[i] = tree->rates->br_r[i];
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Rates_From_Bl(t_tree *tree)
{
  phydbl dt,cr;
  t_node *left, *rght;
  int i;

  dt = -1.0;
  cr = tree->rates->clock_r;

  if(tree->n_root)
    {
      dt = FABS(tree->rates->nd_t[tree->n_root->num] - tree->rates->nd_t[tree->n_root->v[2]->num]);
      tree->rates->br_r[tree->n_root->v[2]->num] = 0.5 * tree->e_root->l->v / (dt*cr);
      dt = FABS(tree->rates->nd_t[tree->n_root->num] - tree->rates->nd_t[tree->n_root->v[1]->num]);
      tree->rates->br_r[tree->n_root->v[1]->num] = 0.5 * tree->e_root->l->v / (dt*cr);
    }
  

  For(i,2*tree->n_otu-3)
    {
      if(tree->a_edges[i] != tree->e_root)
        {
          left = tree->a_edges[i]->left;
          rght = tree->a_edges[i]->rght;
          dt = FABS(tree->rates->nd_t[left->num] - tree->rates->nd_t[rght->num]);         
          
          if(left->anc == rght) tree->rates->br_r[left->num] = tree->a_edges[i]->l->v / (dt*cr);
          else                  tree->rates->br_r[rght->num] = tree->a_edges[i]->l->v / (dt*cr);
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Lk_Jumps(t_tree *tree)
{
  int i,n_jps;
  phydbl dens,dt,lexp;
  t_node *n;

  n = NULL;
  lexp = tree->rates->lexp;
  n_jps = 0;
  dt = 0.0;
  dens = 0.0;

  For(i,2*tree->n_otu-2)
    {
      n = tree->a_nodes[i];
      dt = FABS(tree->rates->nd_t[n->num]-tree->rates->nd_t[n->anc->num]);
      n_jps = tree->rates->n_jps[n->num];
      dens += LOG(Dpois(n_jps,lexp*dt));
    }

  tree->rates->c_lnL_jps = dens;
  
  return dens;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Posterior_Rates(t_tree *tree)
{
  int node_num; 
  node_num = Rand_Int(0,2*tree->n_otu-3);      
  RATES_Posterior_One_Rate(tree->a_nodes[node_num]->anc,tree->a_nodes[node_num],NO,tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Posterior_Times(t_tree *tree)
{
  int node_num;
  node_num = Rand_Int(tree->n_otu,2*tree->n_otu-3);
  RATES_Posterior_One_Time(tree->a_nodes[node_num]->anc,tree->a_nodes[node_num],NO,tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Posterior_One_Rate(t_node *a, t_node *d, int traversal, t_tree *tree)
{
  phydbl like_mean, like_var;
  phydbl prior_mean, prior_var;
  phydbl post_mean, post_var, post_sd;
  phydbl dt,rd,cr,cel,cvl;
  int dim;
  phydbl l_opp; /* length of the branch connected to the root, opposite to the one connected to d */
  t_edge *b;
  int i;
  t_node *v2,*v3;
  phydbl T0,T1,T2,T3;
  phydbl U0,U1,U2,U3;
  phydbl V1,V2,V3;
  phydbl r_min,r_max;
  int err;
  phydbl ratio,u;
  phydbl cvr, cer;
  phydbl nf;
  phydbl new_lnL_data, cur_lnL_data, new_lnL_rate, cur_lnL_rate;
  phydbl sd1,sd2,sd3;
  phydbl inflate_var;

  if(d == tree->n_root) return;

  dim = 2*tree->n_otu-3;
  err = NO;

  inflate_var = tree->rates->inflate_var;

  b = NULL;
  if(a == tree->n_root) b = tree->e_root;
  else For(i,3) if(d->v[i] == a) { b = d->b[i]; break; }
  
  v2 = v3 = NULL;
  if(!d->tax)
    {
      For(i,3)
        if((d->v[i] != a) && (d->b[i] != tree->e_root))
          {
            if(!v2) { v2 = d->v[i]; }
            else    { v3 = d->v[i]; }
          }
    }

  T3 = T2 = 0.0;
  T0 = tree->rates->nd_t[a->num];
  T1 = tree->rates->nd_t[d->num];
  U0 = tree->rates->br_r[a->num];
  U1 = tree->rates->br_r[d->num];
  U3 = U2 = -1.0;

  if(!d->tax)
    {
      T2  = tree->rates->nd_t[v2->num];
      T3  = tree->rates->nd_t[v3->num];
      U2  = tree->rates->br_r[v2->num];
      U3  = tree->rates->br_r[v3->num];
    }
  

  V1 = tree->rates->nu * FABS(T1 - T0);
  V2 = tree->rates->nu * FABS(T2 - T1);
  V3 = tree->rates->nu * FABS(T3 - T1);

  dt     = T1 - T0;
  rd     = U1;
  cr     = tree->rates->clock_r;
  r_min  = tree->rates->min_rate;
  r_max  = tree->rates->max_rate;
  nf     = tree->rates->norm_fact;
  sd1    = SQRT(V1);
  sd2    = SQRT(V2);
  sd3    = SQRT(V3);


  /* Likelihood */
  cel=0.0;
  For(i,dim) 
    if(i != b->num) 
        cel += tree->rates->reg_coeff[b->num*dim+i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
  cel += tree->rates->mean_l[b->num];
  cvl = tree->rates->cond_var[b->num];

  l_opp  = 0.0;
  if(a == tree->n_root)
    {
      if(d == a->v[0]) l_opp = tree->rates->cur_l[a->v[1]->num];
      else             l_opp = tree->rates->cur_l[a->v[0]->num]; 
      cel -= l_opp;
    }


  if(isnan(cvl) || isnan(cel) || cvl < .0) 
    {
      For(i,dim) if(i != b->num) printf("\n. reg: %f %f %f nu=%f clock=%f",
                                        tree->rates->reg_coeff[b->num*dim+i],
                                        tree->rates->u_cur_l[i],
                                        tree->rates->mean_l[i],
                                        tree->rates->nu,tree->rates->clock_r);
      PhyML_Printf("\n. cel=%f cvl=%f\n",cel,cvl); 
      PhyML_Printf("\n. Warning: invalid expected and/or std. dev. values. Skipping this step.\n"); 
      Exit("\n");
    }

  /* Model rates */
  /* if(tree->mod->log_l == YES) cel = EXP(cel); */
  /* like_mean = cel / (dt*cr*nf); */
  /* like_var  = cvl / POW(dt*cr*nf,2);  */

  /* Model log of rates. Branch lengths are given in log units */
  if(tree->rates->model_log_rates == YES) 
    {      
      like_mean = cel -    LOG(dt*cr*nf);
      like_var  = cvl - 2.*LOG(dt*cr*nf);
    }
  else
    {
      like_mean = cel / (dt*cr*nf);
      like_var  = cvl / POW(dt*cr*nf,2);
    }
  
  /* Prior */
  if(!d->tax)
    {
      cvr = 1./(1./V1 + 1./V2 + 1./V3);
      cer = cvr*(U0/V1 + U2/V2 + U3/V3);
    }
  else
    {
      cvr = V1;
      cer = U0;
    }
  
  if(cvr < 0.0)
    {
      PhyML_Printf("\n. cvr=%f d->tax=%d V1=%f v2=%f V3=%f",cvr,d->tax,V1,V2,V3);
      PhyML_Printf("\n. T0=%f T1=%f T2=%f T3=%f",T0,T1,T2,T3);
      Exit("\n");
    }

  prior_mean = cer;
  prior_var  = cvr;

  /* Posterior */
  post_mean = (prior_mean/prior_var + like_mean/like_var)/(1./prior_var + 1./like_var);
  post_var  = 1./(1./prior_var + 1./like_var);


  /* Sample according to priors */
  if(tree->mcmc->use_data == NO)
    {
      post_mean = prior_mean;
      post_var  = prior_var;
    }

  post_sd = SQRT(post_var);

  rd = Rnorm_Trunc(post_mean,inflate_var*post_sd,r_min,r_max,&err);

  if(err || isnan(rd))
    {
      PhyML_Printf("\n");
      PhyML_Printf("\n. run: %d err=%d d->tax=%d",tree->mcmc->run,err,d->tax);
      PhyML_Printf("\n. rd=%f cvr=%G dt=%G cr=%G",rd,cvr,dt,cr);
      PhyML_Printf("\n. prior_mean = %G prior_var = %G",prior_mean,prior_var);
      PhyML_Printf("\n. like_mean = %G like_var = %G",like_mean,like_var);
      PhyML_Printf("\n. post_mean = %G post_var = %G",post_mean,post_var);
      PhyML_Printf("\n. clock_r = %f",tree->rates->clock_r);
      PhyML_Printf("\n. T0=%f T1=%f T2=%f T3=%f",T0,T1,T2,T3);
      PhyML_Printf("\n. U0=%f U1=%f U2=%f U3=%f",U0,U1,U2,U3);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }    

  /* !!!!!!!!!!!!!! */
/*   u = Uni(); */
/*   rd = U1 * EXP(1.*(u-0.5)); */

  if(rd > r_min && rd < r_max)
    {
      
      cur_lnL_data = tree->c_lnL;
      cur_lnL_rate = tree->rates->c_lnL_rates;
      new_lnL_data = tree->c_lnL;
      new_lnL_rate = tree->rates->c_lnL_rates;

      tree->rates->br_r[d->num] = rd;
      RATES_Update_Norm_Fact(tree);
      RATES_Update_Cur_Bl(tree);
      
      if(tree->mcmc->use_data) new_lnL_data = Lk(b,tree);
        /* new_lnL_rate = RATES_Lk_Rates(tree); */
      new_lnL_rate = 
        cur_lnL_rate - 
        (Log_Dnorm_Trunc(U1,U0,sd1,r_min,r_max,&err)) +
        (Log_Dnorm_Trunc(rd,U0,sd1,r_min,r_max,&err));

      if(!d->tax) 
        {
          new_lnL_rate -= (Log_Dnorm_Trunc(U2,U1,sd2,r_min,r_max,&err) + Log_Dnorm_Trunc(U3,U1,sd3,r_min,r_max,&err));
          new_lnL_rate += (Log_Dnorm_Trunc(U2,rd,sd2,r_min,r_max,&err) + Log_Dnorm_Trunc(U3,rd,sd3,r_min,r_max,&err));
        }

      tree->rates->c_lnL_rates = new_lnL_rate;
      
      /* printf("\n. %f %f sd1=%f U1=%f rd=%f ra=%f a=%d d=%d [%f] [%f %f]", */
      /*             new_lnL_rate,RATES_Lk_Rates(tree),sd1,U1,rd,ra,a->num,d->num,tree->rates->br_r[tree->n_root->num], */
      /*             Log_Dnorm_Trunc(U1,ra,sd1,r_min,r_max,&err), */
      /*             Log_Dnorm_Trunc(rd,ra,sd1,r_min,r_max,&err)); */

      ratio = 0.0;
      /* Proposal ratio */
      ratio += (Log_Dnorm_Trunc(U1,post_mean,inflate_var*post_sd,r_min,r_max,&err) - Log_Dnorm_Trunc(rd,post_mean,inflate_var*post_sd,r_min,r_max,&err));
      /*   ratio += LOG(rd/U1); */
      /* Prior ratio */
      ratio += (new_lnL_rate - cur_lnL_rate);
      /* Likelihood ratio */
      ratio += (new_lnL_data - cur_lnL_data);
      
      
      ratio = EXP(ratio);
      
      /*   printf("\n. R a=%3d T0=%6.1f T1=%6.1f T2=%6.1f T3=%6.1f ratio=%8f pm=%7f U1=%7.2f rd=%7.2f %f %f lr=%f %f ld=%f %f [%f]",a->num,T0,T1,T2,T3,ratio,post_mean,U1,rd, */
      /*         Log_Dnorm_Trunc(U1,post_mean,post_sd,r_min,r_max,&err), */
      /*         Log_Dnorm_Trunc(rd,post_mean,post_sd,r_min,r_max,&err), */
      /*         new_lnL_rate,cur_lnL_rate, */
      /*         new_lnL_data,cur_lnL_data, */
      /*         ((Pnorm(r_max,U1,sd2)-Pnorm(r_min,U1,sd2)) * */
      /*          (Pnorm(r_max,U1,sd3)-Pnorm(r_min,U1,sd3)))/ */
      /*         ((Pnorm(r_max,rd,sd2)-Pnorm(r_min,rd,sd2)) * */
      /*          (Pnorm(r_max,rd,sd3)-Pnorm(r_min,rd,sd3)))); */
      
      
      u = Uni();
      
      if(u > MIN(1.,ratio))
        {
          tree->rates->br_r[d->num] = U1; /* reject */
          tree->rates->c_lnL_rates        = cur_lnL_rate;
          tree->c_lnL               = cur_lnL_data;
          RATES_Update_Cur_Bl(tree);
          Update_PMat_At_Given_Edge(b,tree);
        }
      else 
        {
          tree->mcmc->acc_move[tree->mcmc->num_move_br_r+d->num]++;
        }

      RATES_Update_Norm_Fact(tree);
      tree->mcmc->acc_rate[tree->mcmc->num_move_br_r+d->num] = 
        (tree->mcmc->acc_move[tree->mcmc->num_move_br_r+d->num]+1.E-6)/
        (tree->mcmc->run_move[tree->mcmc->num_move_br_r+d->num]+1.E-6);
    }

  tree->mcmc->run_move[tree->mcmc->num_move_br_r+d->num]++;
      

  if(traversal == YES)
    {
      if(d->tax == YES) return;
      else
        {
          For(i,3)
            if(d->v[i] != a && d->b[i] != tree->e_root)
              {
                if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,d->b[i],d);
                /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) { tree->both_sides = YES; Lk(tree); } */
                RATES_Posterior_One_Rate(d,d->v[i],YES,tree);
              }
        }
      
      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,b,d);
      /* if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) { tree->both_sides = YES; Lk(tree); } */
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Posterior_One_Time(t_node *a, t_node *d, int traversal, t_tree *tree)
{
  /*
           T0 (a)
            |
            | l1,U1,b1
            |
           T1 (d)
           / \
 l2,u2,b2 /   \ L3,u3,b3
         /     \
       t2       t3
      (v2)     (v3)

   */

  phydbl l1,l2,l3;
  phydbl new_l1;
  phydbl El1,El2,El3;
  phydbl u0,r1,r2,r3;
  phydbl t0,t1,t2,t3;
  phydbl t_min, t_max;
/*   phydbl t_max_12,t_max_13; */
  phydbl bl_min, bl_max;
  phydbl t1_new;
  phydbl EX,EY;
  phydbl VX,VY;
  phydbl cr;
  int    i,j;
  phydbl *mu, *cov;
  phydbl *cond_mu, *cond_cov;
  short int *is_1;
  phydbl sig11, sig1X, sig1Y, sigXX, sigXY, sigYY;
  phydbl cov11,cov12,cov13,cov22,cov23,cov33;
  int dim;
  t_edge *b1, *b2, *b3;
  t_node *v2,*v3;
  phydbl *l2XY;
  phydbl l_opp;
  t_edge *buff_b;
  t_node *buff_n;
  int err;
  int num_1, num_2, num_3;
  phydbl nf;
  phydbl u, ratio;
  phydbl new_lnL_data, cur_lnL_data, new_lnL_rate, cur_lnL_rate;
  int num_move;
  phydbl inflate_var;

  dim = 2*tree->n_otu-3;
  num_move = tree->mcmc->num_move_nd_t+d->num-tree->n_otu;
  inflate_var = tree->rates->inflate_var;

  if(d->tax) return;
  
  if(FABS(tree->rates->t_prior_min[d->num] - tree->rates->t_prior_max[d->num]) < 1.E-10) return;

  l2XY     = tree->rates->_2n_vect2;
  mu       = tree->rates->_2n_vect3;
  cov      = tree->rates->_2n2n_vect1;
  cond_mu  = tree->rates->_2n_vect1;
  cond_cov = tree->rates->_2n2n_vect2;
  is_1     = tree->rates->_2n_vect5;
  err      = NO;

  b1 = NULL;
  if(a == tree->n_root) b1 = tree->e_root;
  else For(i,3) if(d->v[i] == a) { b1 = d->b[i]; break; }

  b2 = b3 = NULL;
  v2 = v3 = NULL;
  For(i,3)
    if((d->v[i] != a) && (d->b[i] != tree->e_root))
      {
        if(!v2) { v2 = d->v[i]; b2 = d->b[i]; }
        else    { v3 = d->v[i]; b3 = d->b[i]; }
      }

  t2 = tree->rates->nd_t[v2->num];
  t3 = tree->rates->nd_t[v3->num];

  buff_n = NULL;
  buff_b = NULL;
  if(t3 > t2)
    {
      buff_n = v2;
      v2     = v3;
      v3     = buff_n;

      buff_b = b2;
      b2     = b3;
      b3     = buff_b;
    }  
    
  t0 = tree->rates->nd_t[a->num];
  t1 = tree->rates->nd_t[d->num];
  t2 = tree->rates->nd_t[v2->num];
  t3 = tree->rates->nd_t[v3->num];
  u0 = tree->rates->br_r[a->num];
  r1 = tree->rates->br_r[d->num];
  r2 = tree->rates->br_r[v2->num];
  r3 = tree->rates->br_r[v3->num];
  l1 = tree->rates->cur_l[d->num];
  l2 = tree->rates->cur_l[v2->num];
  l3 = tree->rates->cur_l[v3->num];
  cr = tree->rates->clock_r;
  nf = tree->rates->norm_fact;

  For(i,dim) is_1[i] = 0;
  is_1[b1->num] = 1;
  is_1[b2->num] = 1;
  is_1[b3->num] = 1;

/*   For(i,dim)     cond_mu[i]  = 0.0; */
/*   For(i,dim*dim) cond_cov[i] = 0.0; */
/*   Normal_Conditional(tree->rates->mean_l,tree->rates->cov_l,tree->rates->u_cur_l,dim,is_1,3,cond_mu,cond_cov); */
  
/*   El1 = cond_mu[b1->num]; */
/*   El2 = cond_mu[b2->num]; */
/*   El3 = cond_mu[b3->num]; */

/*   cov11 = cond_cov[b1->num*dim+b1->num]; */
/*   cov12 = cond_cov[b1->num*dim+b2->num]; */
/*   cov13 = cond_cov[b1->num*dim+b3->num]; */
/*   cov23 = cond_cov[b2->num*dim+b3->num]; */
/*   cov22 = cond_cov[b2->num*dim+b2->num]; */
/*   cov33 = cond_cov[b3->num*dim+b3->num]; */


/*   El1 = tree->rates->u_ml_l[b1->num]; */
/*   El2 = tree->rates->u_ml_l[b2->num]; */
/*   El3 = tree->rates->u_ml_l[b3->num]; */

/*   cov11 = tree->rates->cov[b1->num*dim+b1->num]; */
/*   cov12 = tree->rates->cov[b1->num*dim+b2->num]; */
/*   cov13 = tree->rates->cov[b1->num*dim+b3->num]; */
/*   cov23 = tree->rates->cov[b2->num*dim+b3->num]; */
/*   cov22 = tree->rates->cov[b2->num*dim+b2->num]; */
/*   cov33 = tree->rates->cov[b3->num*dim+b3->num]; */

/*   PhyML_Printf("\n- El1=%10f El2=%10f El3=%10f",El1,El2,El3); */
/*   PhyML_Printf("\n- cov11=%10f cov12=%10f cov13=%10f cov23=%10f cov22=%10f cov33=%10f", cov11,cov12,cov13,cov23,cov22,cov33); */

  num_1 = num_2 = num_3 = -1;
  if(b1->num < b2->num && b2->num < b3->num) { num_1 = 0; num_2 = 1; num_3 = 2; }
  if(b1->num < b3->num && b3->num < b2->num) { num_1 = 0; num_2 = 2; num_3 = 1; }
  if(b2->num < b1->num && b1->num < b3->num) { num_1 = 1; num_2 = 0; num_3 = 2; }
  if(b2->num < b3->num && b3->num < b1->num) { num_1 = 2; num_2 = 0; num_3 = 1; }
  if(b3->num < b2->num && b2->num < b1->num) { num_1 = 2; num_2 = 1; num_3 = 0; }
  if(b3->num < b1->num && b1->num < b2->num) { num_1 = 1; num_2 = 2; num_3 = 0; }

  cov11 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_1];
  cov12 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_2];
  cov13 = tree->rates->trip_cond_cov[d->num * 9 + num_1 * 3 + num_3];
  cov23 = tree->rates->trip_cond_cov[d->num * 9 + num_2 * 3 + num_3];
  cov22 = tree->rates->trip_cond_cov[d->num * 9 + num_2 * 3 + num_2];
  cov33 = tree->rates->trip_cond_cov[d->num * 9 + num_3 * 3 + num_3];

  El1=0.0;
  For(i,dim)
    if(i != b1->num && i != b2->num && i != b3->num)
      El1 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_1 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
  El1 += tree->rates->mean_l[b1->num];

  El2=0.0;
  For(i,dim)
    if(i != b1->num && i != b2->num && i != b3->num)
      El2 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_2 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
  El2 += tree->rates->mean_l[b2->num];

  El3=0.0;
  For(i,dim)
    if(i != b1->num && i != b2->num && i != b3->num)
      El3 += tree->rates->trip_reg_coeff[d->num * (6*tree->n_otu-9) + num_3 * dim +i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
  El3 += tree->rates->mean_l[b3->num];


/*   PhyML_Printf("\n+ El1=%10f El2=%10f El3=%10f",El1,El2,El3); */
/*   PhyML_Printf("\n+ cov11=%10f cov12=%10f cov13=%10f cov23=%10f cov22=%10f cov33=%10f", cov11,cov12,cov13,cov23,cov22,cov33); */


  t1_new = +1;

  t_min = MAX(t0,tree->rates->t_prior_min[d->num]);
  t_max = MIN(MIN(t2,t3),tree->rates->t_prior_max[d->num]);
    
  t_min += tree->rates->min_dt;
  t_max -= tree->rates->min_dt;

  if(t_min > t_max) 
    {
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  bl_min = bl_max = -1.0;

  l_opp = 0.0;
  if(a == tree->n_root)
    {
      l_opp = (d == a->v[0])?(tree->rates->cur_l[a->v[1]->num]):(tree->rates->cur_l[a->v[0]->num]);
      El1 -= l_opp;
    }
     
  EX = El1/r1 + El2/r2;
  EY = El1/r1 + El3/r3;
  
  VX = cov11/(r1*r1) + cov22/(r2*r2) + 2.*cov12/(r1*r2);
  VY = cov11/(r1*r1) + cov33/(r3*r3) + 2.*cov13/(r1*r3);

  mu[0] = El1;
  mu[1] = EX;
  mu[2] = EY;
  
  sig11 = cov11;
  sig1X = cov11/r1 + cov12/r2;
  sig1Y = cov11/r1 + cov13/r3;
  sigXX = VX;
  sigYY = VY;
  sigXY = cov11/(r1*r1) + cov13/(r1*r3) + cov12/(r1*r2) + cov23/(r2*r3);

  cov[0*3+0] = sig11; cov[0*3+1] = sig1X; cov[0*3+2] = sig1Y;
  cov[1*3+0] = sig1X; cov[1*3+1] = sigXX; cov[1*3+2] = sigXY;
  cov[2*3+0] = sig1Y; cov[2*3+1] = sigXY; cov[2*3+2] = sigYY;

  l2XY[0] = 0.0; /* does not matter */
  l2XY[1] = (t2-t0)*cr*nf; /* constraint 1  */
  l2XY[2] = (t3-t0)*cr*nf; /* constraint 2  */

  is_1[0] = is_1[1] = is_1[2] = 0;
  is_1[0] = 1;

  Normal_Conditional(mu,cov,l2XY,3,is_1,1,cond_mu,cond_cov);


  if(cond_cov[0*3+0] < 0.0)
    {
      PhyML_Printf("\n. a: %d d: %d",a->num,d->num);
      PhyML_Printf("\n. Conditional mean=%G var=%G",cond_mu[0],cond_cov[0*3+0]);
      PhyML_Printf("\n. t0=%G t1=%f t2=%f t3=%f l1=%G l2=%G l3=%G",t0,t1,t2,t3,l1,l2,l3);
      PhyML_Printf("\n. El1=%G El2=%G El3=%G Nu=%G r1=%G r2=%G r3=%G cr=%G",El1,El2,El3,tree->rates->nu,r1,r2,r3,cr);
      PhyML_Printf("\n. COV11=%f COV12=%f COV13=%f COV22=%f COV23=%f COV33=%f",cov11,cov12,cov13,cov22,cov23,cov33);
      PhyML_Printf("\n. constraint1: %f constraints2: %f",l2XY[1],l2XY[2]);
      PhyML_Printf("\n");
      For(i,3)
        {
          PhyML_Printf(". mu%d=%12lf\t",i,mu[i]);
          For(j,3)
            {
              PhyML_Printf("%12lf ",cov[i*3+j]);
            }
          PhyML_Printf("\n");
        }      
      cond_cov[0*3+0] = 1.E-10;
    }


  bl_min = (t_min - t0) * r1 * cr * nf;
  bl_max = (t_max - t0) * r1 * cr * nf;

  new_l1 = Rnorm_Trunc(cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err);
/*   new_l1 = Rnorm(cond_mu[0],SQRT(cond_cov[0*3+0])); */


  if(new_l1 < bl_min) new_l1 = l1;
  if(new_l1 > bl_max) new_l1 = l1;
      
  t1_new = new_l1/(r1*cr*nf) + t0;


  if(err)
    {
      PhyML_Printf("\n");
      PhyML_Printf("\n. Root ? %s",(tree->n_root==a)?("yes"):("no")); 
      PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run);
      PhyML_Printf("\n. t0=%f t1=%f t2=%f t3=%f",t0,t1,t2,t3);
      PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
      PhyML_Printf("\n. bl_min=%f bl_max=%f",bl_min,bl_max);
      PhyML_Printf("\n. cond_mu[0]=%f cond_cov[0]=%f",cond_mu[0],SQRT(cond_cov[0*3+0]));
      PhyML_Printf("\n. El1=%f El2=%f El3=%f",El1,El2,El3);
      PhyML_Printf("\n. l1=%f l2=%f l3=%f",l1,l2,l3);
      PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
      PhyML_Printf("\n. COV11=%f COV22=%f",cov11,cov22);
      PhyML_Printf("\n. Clock rate = %f",tree->rates->clock_r);
      PhyML_Printf("\n. Setting t1_new to %f",t1);
      t1_new = t1;
      /*          Exit("\n"); */
    }
  /*     } */
  /*   else */
  /*     { */
  /*       bl_min = (t2 - t_max) * r2; */
  /*       bl_max = (t2 - t_min) * r2; */
  
  /*       l2 = Rnorm_Trunc(cond_mu[0],SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err); */
  /*       t1_new = -l2/r2 + t2; */
  
  /*       if(err) */
  /*    { */
  /*      PhyML_Printf("\n"); */
  /*      PhyML_Printf("\n. %s %d %d",__FILE__,__LINE__,tree->mcmc->run); */
  /*      PhyML_Printf("\n. t0=%f t1=%f t2=%f t3=%f",t0,t1,t2,t3); */
  /*      PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max); */
  /*      PhyML_Printf("\n. bl_min=%f bl_max=%f",bl_min,bl_max); */
  /*      PhyML_Printf("\n. cond_mu[0]=%f cond_cov[0]=%f",cond_mu[0],SQRT(cond_cov[0*3+0])); */
  /*      PhyML_Printf("\n. El1=%f El2=%f El3=%f",El1,El2,El3); */
  /*      PhyML_Printf("\n. l1=%f l2=%f l3=%f",l1,l2,l3); */
  /*      PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3); */
  /*      PhyML_Printf("\n. COV11=%f COV22=%f",cov11,cov22); */
  /*      PhyML_Printf("\n. Clock rate = %f",tree->rates->clock_r); */
  /*      PhyML_Printf("\n. Setting t1_new to %f",t1); */
  /*      t1_new = t1; */
  /* /\*          Exit("\n"); *\/ */
  /*    } */
  /*     } */
  
  if(t1_new < t0)
    {
      t1_new = t0+1.E-4;
      PhyML_Printf("\n");
      PhyML_Printf("\n. a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
      PhyML_Printf("\n. t0 = %f t1_new = %f t1 = %f",t0,t1_new,t1);
      PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
      PhyML_Printf("\n. l1 = %f",l1);
      PhyML_Printf("\n. bl_min = %f bl_max = %f",bl_min,bl_max);
      PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1-t0,t2-t1);
      PhyML_Printf("\n. l1 = %f l2 = %f cov11=%f cov22=%f cov33=%f",l1,l2,cov11,cov22,cov33);
      PhyML_Printf("\n. clock=%G",tree->rates->clock_r);
      PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      /*       Exit("\n"); */
    }
  if(t1_new > MIN(t2,t3))
    {
      PhyML_Printf("\n");
      PhyML_Printf("\n. a is root -> %s",(a == tree->n_root)?("YES"):("NO"));
      PhyML_Printf("\n. t0 = %f t1_new = %f t1 = %f t2 = %f t3 = %f MIN(t2,t3)=%f",t0,t1_new,t1,t2,t3,MIN(t2,t3));
      PhyML_Printf("\n. t_min=%f t_max=%f",t_min,t_max);
      PhyML_Printf("\n. l2 = %f",l2);
      PhyML_Printf("\n. bl_min = %f bl_max = %f",bl_min,bl_max);
      PhyML_Printf("\n. (t1-t0)=%f (t2-t1)=%f",t1-t0,t2-t1);
      PhyML_Printf("\n. l1 = %f l2 = %f cov11=%f cov22=%f cov33=%f",l1,l2,cov11,cov22,cov33);
      PhyML_Printf("\n. clock=%G",tree->rates->clock_r);
      PhyML_Printf("\n. u0=%f r1=%f r2=%f r3=%f",u0,r1,r2,r3);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      /*       Exit("\n"); */
    }
  
  if(isnan(t1_new))
    {
      PhyML_Printf("\n. run=%d",tree->mcmc->run);
      PhyML_Printf("\n. mean=%G var=%G",cond_mu[0],cond_cov[0*3+0]);
      PhyML_Printf("\n. t1=%f l1=%G r1=%G t0=%G",t1,l1,r1,t0);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      /*       Exit("\n"); */
    }
  
  if(tree->mcmc->use_data == YES)
    tree->rates->nd_t[d->num] = t1_new;
  else
    {
      u = Uni();
      tree->rates->nd_t[d->num] = u*(t_max-t_min) + t_min;
    }

  RATES_Update_Norm_Fact(tree);
  RATES_Update_Cur_Bl(tree);
  
  cur_lnL_data = tree->c_lnL;
  cur_lnL_rate = tree->rates->c_lnL_rates;
  new_lnL_data = tree->c_lnL;
  new_lnL_rate = tree->rates->c_lnL_rates;
  
  if(tree->mcmc->use_data) 
    {
      /* !!!!!!!!!!!!!!!!! */
      /* tree->both_sides = NO; */
      /* new_lnL_data = Lk(tree); */

      if(tree->io->lk_approx == EXACT)
        {
          Update_PMat_At_Given_Edge(b1,tree);
          Update_PMat_At_Given_Edge(b2,tree);
          Update_PMat_At_Given_Edge(b3,tree);
          Update_P_Lk(tree,b1,d);
        }
      new_lnL_data = Lk(b1,tree);
    }
  
  new_lnL_rate = RATES_Lk_Rates(tree);

  ratio = 0.0;

  /* Proposal ratio */
  if(tree->mcmc->use_data)
    ratio += (Log_Dnorm_Trunc(l1,    cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err) - 
              Log_Dnorm_Trunc(new_l1,cond_mu[0],inflate_var*SQRT(cond_cov[0*3+0]),bl_min,bl_max,&err));
    
  /* Prior ratio */
  ratio += (new_lnL_rate - cur_lnL_rate);

  /* Likelihood ratio */
  ratio += (new_lnL_data - cur_lnL_data);

  /*   printf("\n* d:%d Ratio=%f l1=%f new_l1=%f mean=%f ml=%f sd=%f [%f %f]", */
  /*     d->num, */
  /*     EXP(ratio), */
  /*     l1,new_l1, */
  /*     cond_mu[0], */
  /*     tree->rates->mean_l[b1->num], */
  /*     SQRT(cond_cov[0*3+0]), */
  /*     Log_Dnorm(l1,cond_mu[0],SQRT(cond_cov[0*3+0]),&err),Log_Dnorm(new_l1,cond_mu[0],SQRT(cond_cov[0*3+0]),&err)); */
  
  ratio = EXP(ratio);


  u = Uni();
  if(u > MIN(1.,ratio))
    {
      tree->rates->nd_t[d->num] = t1; /* reject */
      tree->rates->c_lnL_rates        = cur_lnL_rate;
      tree->c_lnL               = cur_lnL_data;
      RATES_Update_Cur_Bl(tree);
      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) 
        {
          Update_PMat_At_Given_Edge(b1,tree);
          Update_PMat_At_Given_Edge(b2,tree);
          Update_PMat_At_Given_Edge(b3,tree);
          Update_P_Lk(tree,b1,d);
        }
    }
  else
    {
      tree->mcmc->acc_move[num_move]++;
    }
  
  RATES_Update_Norm_Fact(tree);

  tree->mcmc->run_move[num_move]++;
  
  if(traversal == YES)
    {
      if(d->tax == YES) return;
      else
        {
          For(i,3)
            if(d->v[i] != a && d->b[i] != tree->e_root)
              {
                if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,d->b[i],d);
                RATES_Posterior_One_Time(d,d->v[i],YES,tree);
              }
        }
      if(tree->io->lk_approx == EXACT && tree->mcmc->use_data) Update_P_Lk(tree,b1,d);
    }
  
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Posterior_Time_Root(t_tree *tree)
{
  /* 
        Root, T0
        / \
  l1   /   \ l2
      /     \
     v1,t1   v2,t2

  */

  phydbl t1,t2;
  phydbl u0,u1,u2;
  phydbl cel,cvl;
  int i,dim;
  t_edge *b;
  t_node *root;
  phydbl t0,t0_min, t0_max;
  phydbl new_t0;
/*   phydbl t_max_01, t_max_02; */
  int err;
  phydbl cr;
  phydbl bl_min, bl_max;
  phydbl new_l;
  phydbl new_lnL_data, cur_lnL_data, cur_lnL_rate;
  phydbl u,ratio;

  dim    = 2*tree->n_otu-3;
  b      = tree->e_root;
  root   = tree->n_root;
  t0     = tree->rates->nd_t[root->num];
  t1     = tree->rates->nd_t[root->v[2]->num];
  t2     = tree->rates->nd_t[root->v[1]->num];
  u1     = tree->rates->br_r[root->v[2]->num];
  u2     = tree->rates->br_r[root->v[1]->num];
  u0     = 1.0;
  cr     = tree->rates->clock_r;
  t0_min = -BIG;
  t0_max = MIN(t1,t2);
  
  /*   t0_max = MIN(t1 - (1./tree->rates->nu)*POW((u0-u1)/tree->rates->z_max,2), */
  /*           t2 - (1./tree->rates->nu)*POW((u0-u2)/tree->rates->z_max,2)); */
  
  /*   if(u1 > u0) t_max_01 = t1 - (1./tree->rates->nu)*POW((u1-u0)/(PointNormal(tree->rates->p_max*Pnorm(.0,.0,1.))),2); */
  /*   else        t_max_01 = t1 - (1./tree->rates->nu)*POW((u1-u0)/(PointNormal(tree->rates->p_max*(1.-Pnorm(.0,.0,1.)))),2); */
  
  /*   if(u2 > u0) t_max_02 = t2 - (1./tree->rates->nu)*POW((u2-u0)/(PointNormal(tree->rates->p_max*Pnorm(.0,.0,1.))),2); */
  /*   else        t_max_02 = t2 - (1./tree->rates->nu)*POW((u2-u0)/(PointNormal(tree->rates->p_max*(1.-Pnorm(.0,.0,1.)))),2); */
  
  
  /*   t_max_01 = RATES_Find_Max_Dt_Bisec(u1,u0,tree->rates->t_prior_min[root->num],t1,tree->rates->nu,tree->rates->p_max,(u1 < u0)?YES:NO); */
  /*   t_max_02 = RATES_Find_Max_Dt_Bisec(u2,u0,tree->rates->t_prior_min[root->num],t2,tree->rates->nu,tree->rates->p_max,(u2 < u0)?YES:NO); */
  /*   t0_max = MIN(t_max_01,t_max_02); */
    
  /*   RATES_Min_Max_Interval(u0,u1,u2,r3,t0,t2,t3,&t_min,&t_max,tree->rates->nu,tree->rates->p_max,tree); */

  t0_min = MAX(t0_min,tree->rates->t_prior_min[root->num]);
  t0_max = MIN(t0_max,tree->rates->t_prior_max[root->num]);

  t0_max -= tree->rates->min_dt;

  if(t0_min > t0_max) return;

  tree->rates->t_prior[root->num] = Uni()*(t0_max - t0_min) + t0_min;

  u0 *= cr;
  u1 *= cr;
  u2 *= cr;

  if(FABS(tree->rates->t_prior_min[root->num] - tree->rates->t_prior_max[root->num]) < 1.E-10) return;

  cel=0.0;
  For(i,dim) if(i != b->num) cel += tree->rates->reg_coeff[b->num*dim+i] * (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]);
  cel += tree->rates->mean_l[b->num];

  cvl = tree->rates->cond_var[b->num];
  
  bl_min = u1 * (t1 - t0_max) + u2 * (t2 - t0_max); 
  bl_max = u1 * (t1 - t0_min) + u2 * (t2 - t0_min);

  if(bl_min > bl_max) return;

  new_l = Rnorm_Trunc(cel,SQRT(cvl),bl_min,bl_max,&err);

  new_t0 = (u1*t1 + u2*t2 - new_l)/(u1+u2);

  if(t0 > t1 || t0 > t2)
    {
      PhyML_Printf("\n");
      PhyML_Printf("\n. Run = %d",tree->mcmc->run);
      PhyML_Printf("\n. t0=%f t1=%f t2=%f",t0,t1,t2);
      PhyML_Printf("\n. t0_min=%f t0_max=%f",t0_min,t0_max);
      PhyML_Printf("\n. new_l=%f cel=%f",new_l,cel);
      PhyML_Printf("\n. u0=%f u1=%f u2=%f",u0/cr,u1/cr,u2/cr);
      PhyML_Printf("\n. Nu = %f Clock = %f",tree->rates->nu,tree->rates->clock_r);
      PhyML_Printf("\n. Setting t0 to %f",tree->rates->nd_t[root->num]);
      return;
    }

  if(t0 < t0_min || t0 > t0_max)
    {
      PhyML_Printf("\n");
      PhyML_Printf("\n. Run = %d",tree->mcmc->run);
      PhyML_Printf("\n. t0=%f t1=%f t2=%f",t0,t1,t2);
      PhyML_Printf("\n. t0_min=%f t0_max=%f",t0_min,t0_max);
      PhyML_Printf("\n. u0=%f u1=%f u2=%f",u0/cr,u1/cr,u2/cr);
      PhyML_Printf("\n. Nu = %f Clock = %f",tree->rates->nu,tree->rates->clock_r);
      PhyML_Printf("\n. Setting t0 to %f",tree->rates->nd_t[root->num]);
      t0 = tree->rates->nd_t[root->num];      
    }

  /* Sample according to prior */
/*   tree->rates->nd_t[root->num] = tree->rates->t_prior[root->num]; */

  /* Sample according to posterior */
  tree->rates->nd_t[root->num] = new_t0;

  RATES_Update_Norm_Fact(tree);
  RATES_Update_Cur_Bl(tree);

  cur_lnL_data = tree->c_lnL;
  cur_lnL_rate = tree->rates->c_lnL_rates;
  new_lnL_data = tree->c_lnL;
  
  if(tree->mcmc->use_data) new_lnL_data = Lk(NULL,tree);
    
  ratio = 0.0;
  /* Prior ratio */
  ratio += .0;
  /* Likelihood ratio */
  ratio += new_lnL_data - cur_lnL_data;

  ratio = EXP(ratio);
  u = Uni();
  if(u > MIN(1.,ratio))
    {
      tree->rates->nd_t[root->num] = t0; /* reject */
      tree->rates->c_lnL_rates     = cur_lnL_rate;
      tree->c_lnL                  = cur_lnL_data;
    }
  else
    {
      tree->mcmc->acc_move[tree->mcmc->num_move_nd_t]++;
    }

  RATES_Update_Norm_Fact(tree);
  RATES_Update_Cur_Bl(tree);

  tree->mcmc->run_move[tree->mcmc->num_move_nd_t]++;
  tree->mcmc->acc_rate[tree->mcmc->num_move_nd_t] = 
    (tree->mcmc->acc_move[tree->mcmc->num_move_nd_t]+1.E-6)/ 
    (tree->mcmc->run_move[tree->mcmc->num_move_nd_t]+1.E+6);

}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Cur_Bl(t_tree *tree)
{
  RATES_Update_Norm_Fact(tree);

  RATES_Update_Cur_Bl_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
  RATES_Update_Cur_Bl_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
  
  if(tree->mod && tree->mod->log_l == YES)
    {
      tree->e_root->l->v = 
        EXP(tree->rates->cur_l[tree->n_root->v[2]->num]) +
        EXP(tree->rates->cur_l[tree->n_root->v[1]->num]) ;
      tree->e_root->l->v = LOG(tree->e_root->l->v);
    }
  else
    {
      tree->e_root->l->v = 
        tree->rates->cur_l[tree->n_root->v[2]->num] +
        tree->rates->cur_l[tree->n_root->v[1]->num] ;
    }

  tree->rates->u_cur_l[tree->e_root->num] = tree->e_root->l->v;
  tree->n_root_pos = tree->rates->cur_l[tree->n_root->v[2]->num] / tree->e_root->l->v;

  if(tree->rates->model == GUINDON)
    {
      phydbl t0,t1,t2;
      t_node *n0, *n1;
      
      n0 = tree->n_root->v[2];
      n1 = tree->n_root->v[1];
      t1 = tree->rates->nd_t[tree->n_root->v[2]->num];
      t2 = tree->rates->nd_t[tree->n_root->v[1]->num];
      t0 = tree->rates->nd_t[tree->n_root->num];

      tree->e_root->l->v = 
        (t1-t0)/(t1+t2-2.*t0)*tree->rates->cur_gamma_prior_mean[n0->num] +
        (t2-t0)/(t1+t2-2.*t0)*tree->rates->cur_gamma_prior_mean[n1->num];
      
      tree->e_root->l_var->v = 
        POW((t1-t0)/(t1+t2-2.*t0),2)*tree->rates->cur_gamma_prior_var[n0->num] +
        POW((t2-t0)/(t1+t2-2.*t0),2)*tree->rates->cur_gamma_prior_var[n1->num];

      /* printf("\n. ROOT: %f %f %f %f", */
      /*             tree->rates->cur_gamma_prior_mean[n0->num], */
      /*             tree->rates->cur_gamma_prior_var[n0->num], */
      /*             tree->rates->cur_gamma_prior_mean[n1->num], */
      /*             tree->rates->cur_gamma_prior_var[n1->num]); */
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Cur_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{
  phydbl dt,rr,cr,ra,rd,ta,td,nu;

  tree->rates->br_do_updt[d->num] = YES;

  if(tree->rates->br_do_updt[d->num] == YES)
    {
      tree->rates->br_do_updt[d->num] = NO;

      dt = tree->rates->nd_t[d->num] - tree->rates->nd_t[a->num];
      cr = tree->rates->clock_r;
      rd = tree->rates->br_r[d->num];
      ra = tree->rates->br_r[a->num];
      td = tree->rates->nd_t[d->num];
      ta = tree->rates->nd_t[a->num];
      nu = tree->rates->nu;


      if(tree->rates->model_log_rates == YES)
        {
          /* Artihmetic average */
          rr = (EXP(ra) + EXP(rd))/2.;
        }
      else
        {
          rr = (ra+rd)/2.;
        }

      tree->rates->cur_l[d->num] = dt*rr*cr;
      
      if(tree->rates->model == GUINDON)
        {
          phydbl m,v;

          Integrated_Geometric_Brownian_Bridge_Moments(dt,ra,rd,nu,&m,&v);
          
          m *= cr*dt; // the actual rate average is m * cr. We multiply by dt in order to derive the value for the branch length
          v *= (cr*cr)*(dt*dt);

          tree->rates->cur_gamma_prior_mean[d->num] = m;
          tree->rates->cur_gamma_prior_var[d->num]  = v;

          tree->rates->cur_l[d->num] = tree->rates->cur_gamma_prior_mean[d->num]; // Required for having proper branch lengths in Write_Tree function
        }
      
      if(tree->mod && tree->mod->log_l == YES) tree->rates->cur_l[d->num] = LOG(tree->rates->cur_l[d->num]);
      
      if(b)
        {
          b->l->v                      = tree->rates->cur_l[d->num];
          tree->rates->u_cur_l[b->num] = tree->rates->cur_l[d->num];
          b->l_var->v                  = tree->rates->cur_gamma_prior_var[d->num];
        }
      
      if(b && (isnan(b->l->v) || isnan(b->l_var->v)))
        {
          PhyML_Printf("\n. dt=%G rr=%G cr=%G ra=%G rd=%G nu=%G %f %f ",dt,rr,cr,ra,rd,nu,b->l_var->v,b->l->v);   
          PhyML_Printf("\n. ta=%G td=%G ra*cr=%G rd*cr=%G sd=%G",
                       ta,td,ra*cr,rd*cr,
                       SQRT(dt*nu)*cr);
          PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
          Exit("\n");
        }
    }

  if(d->tax) return;
  else
    {
      int i;
      For(i,3) 
        if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
          RATES_Update_Cur_Bl_Pre(d,d->v[i],d->b[i],tree);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Bl_To_Bl(t_tree *tree)
{
  RATES_Bl_To_Bl_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
  RATES_Bl_To_Bl_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
  /* tree->rates->cur_l[tree->n_root->v[2]->num] = tree->a_edges[tree->e_root->num]->l->v * tree->n_root_pos; */
  /* tree->rates->cur_l[tree->n_root->v[1]->num] = tree->a_edges[tree->e_root->num]->l->v * (1. - tree->n_root_pos); */
  tree->rates->cur_l[tree->n_root->v[2]->num] = tree->a_edges[tree->e_root->num]->l->v * 0.5;
  tree->rates->cur_l[tree->n_root->v[1]->num] = tree->a_edges[tree->e_root->num]->l->v * (1. - 0.5);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Bl_To_Bl_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{

  if(b) 
    {
      tree->rates->cur_l[d->num] = b->l->v;
    }

  if(d->tax) return;
  else
    {
      int i;

      For(i,3) 
        if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
          RATES_Bl_To_Bl_Pre(d,d->v[i],d->b[i],tree);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Bl_To_Ml(t_tree *tree)
{
  RATES_Bl_To_Ml_Pre(tree->n_root,tree->n_root->v[2],NULL,tree);
  RATES_Bl_To_Ml_Pre(tree->n_root,tree->n_root->v[1],NULL,tree);
  tree->rates->u_ml_l[tree->e_root->num] = tree->a_edges[tree->e_root->num]->l->v;
  tree->rates->ml_l[tree->n_root->v[2]->num] = tree->rates->u_ml_l[tree->e_root->num] * tree->n_root_pos;
  tree->rates->ml_l[tree->n_root->v[1]->num] = tree->rates->u_ml_l[tree->e_root->num] * (1. - tree->n_root_pos);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Bl_To_Ml_Pre(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{

  if(b) 
    {
      tree->rates->u_ml_l[b->num] = b->l->v;
      tree->rates->ml_l[d->num]   = b->l->v;
    }

  if(d->tax) return;
  else
    {
      int i;

      For(i,3) 
        if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
          RATES_Bl_To_Ml_Pre(d,d->v[i],d->b[i],tree);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Cov_Matrix_Rooted(phydbl *unroot_cov, t_tree *tree)
{
  int i,dim;

  dim = 2*tree->n_otu-3;

  RATES_Get_Cov_Matrix_Rooted_Pre(tree->n_root,tree->n_root->v[2],NULL,unroot_cov,tree);
  RATES_Get_Cov_Matrix_Rooted_Pre(tree->n_root,tree->n_root->v[1],NULL,unroot_cov,tree);

  For(i,dim+1) tree->rates->cov_l[i*(dim+1)+tree->n_root->v[2]->num] /= 2.;
  For(i,dim+1) tree->rates->cov_l[i*(dim+1)+tree->n_root->v[1]->num] /= 2.;
  For(i,dim+1) tree->rates->cov_l[tree->n_root->v[2]->num*(dim+1)+i] /= 2.;
  For(i,dim+1) tree->rates->cov_l[tree->n_root->v[1]->num*(dim+1)+i] /= 2.;

}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Cov_Matrix_Rooted_Pre(t_node *a, t_node *d, t_edge *b, phydbl *cov, t_tree *tree)
{
  int i, dim;
  t_node *n;

  dim       = 2*tree->n_otu-3;
  n         = NULL;

  For(i,dim) 
    { 
      if(tree->a_edges[i] != tree->e_root)
        {
          n = 
            (tree->a_edges[i]->left->anc == tree->a_edges[i]->rght)?
            (tree->a_edges[i]->left):
            (tree->a_edges[i]->rght);

          if(b)
            {
              tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
            }
          else
            {
              tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];
            }
        }
      else
        {
          n = tree->e_root->left;
          if(b)
            tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
          else
            tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];

          n = tree->e_root->rght;
          if(b)
            tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[b->num*dim + i];
          else
            tree->rates->cov_l[d->num*(dim+1) + n->num] = cov[tree->e_root->num*dim + i];
        }
    }


  if(d->tax) return;
  else
    {
      For(i,3)
        if((d->v[i] != a) && (d->b[i] != tree->e_root)) 
          RATES_Get_Cov_Matrix_Rooted_Pre(d,d->v[i],d->b[i],cov,tree);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Covariance_Mu(t_tree *tree)
{
  int i,j;
  phydbl dt,var;
  int dim;
  int lca_num;

  dim = 2*tree->n_otu-2;

  For(i,dim*dim) tree->rates->cov_r[i] = 0.0;

  dt =  tree->rates->nd_t[tree->n_root->v[2]->num] - tree->rates->nd_t[tree->n_root->num];
  var = dt * tree->rates->nu;
  tree->rates->cov_r[tree->n_root->v[2]->num*dim+tree->n_root->v[2]->num] = var;


  dt = tree->rates->nd_t[tree->n_root->v[1]->num] - tree->rates->nd_t[tree->n_root->num];
  var = dt * tree->rates->nu;
  tree->rates->cov_r[tree->n_root->v[1]->num*dim+tree->n_root->v[1]->num] = var;

  RATES_Variance_Mu_Pre(tree->n_root,tree->n_root->v[2],tree);
  RATES_Variance_Mu_Pre(tree->n_root,tree->n_root->v[1],tree);

  For(i,dim)
    {
      for(j=i+1;j<dim;j++)
        {
          lca_num = tree->rates->lca[i*(dim+1)+j]->num;
          if(lca_num < dim) 
            {
              tree->rates->cov_r[i*dim+j] = tree->rates->cov_r[lca_num*dim+lca_num];        
              tree->rates->cov_r[j*dim+i] = tree->rates->cov_r[i*dim+j];
            }
          else if(lca_num == dim)
            {
              tree->rates->cov_r[i*dim+j] = 0.0;            
              tree->rates->cov_r[j*dim+i] = 0.0;
            }
          else
            {
              PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
              Exit("\n");
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Variance_Mu_Pre(t_node *a, t_node *d, t_tree *tree)
{
  int dim;
  phydbl var0;
  phydbl dt1,var1;
  phydbl dt2,var2;
  int i;
  int dir1, dir2;

  dim = 2*tree->n_otu-2;

  if(d->tax) return;

  dir1 = dir2 = -1;
  For(i,3)
    {
      if((d->v[i] != a) && (d->b[i] != tree->e_root))
        {
          if(dir1 < 0) dir1 = i;
          else         dir2 = i;
        }
    }


  var0 = tree->rates->cov_r[d->num*dim+d->num];

  dt1  = tree->rates->nd_t[d->v[dir1]->num] - tree->rates->nd_t[d->num];
  var1 = tree->rates->nu*dt1;

  dt2  = tree->rates->nd_t[d->v[dir2]->num] - tree->rates->nd_t[d->num];
  var2 = tree->rates->nu*dt2;

  tree->rates->cov_r[d->v[dir1]->num*dim+d->v[dir1]->num] = var0+var1;
  tree->rates->cov_r[d->v[dir2]->num*dim+d->v[dir2]->num] = var0+var2;


  For(i,3)
    {
      if((d->v[i] != a) && (d->b[i] != tree->e_root))
        {
          RATES_Variance_Mu_Pre(d,d->v[i],tree);
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Fill_Lca_Table(t_tree *tree)
{
  int i,j;
  int dim;
  
  dim = 2*tree->n_otu-1;

  For(i,dim)
    {
      for(j=i;j<dim;j++)
        {
          tree->rates->lca[i*dim+j] = Find_Lca_Pair_Of_Nodes(tree->a_nodes[i],tree->a_nodes[j],tree);
          tree->rates->lca[j*dim+i] = tree->rates->lca[i*dim+j];
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

/* Get V(L_{i} | L_{-i}) for all i */
void RATES_Get_Conditional_Variances(t_tree *tree)
{
  int i,j;
  short int *is_1;
  phydbl *a;
  int dim;
  t_edge *b;
  phydbl *cond_mu, *cond_cov;

  dim      = 2*tree->n_otu-3;
  a        = tree->rates->_2n_vect1;
  is_1     = tree->rates->_2n_vect5;
  b        = NULL;
  cond_mu  = tree->rates->_2n_vect2;
  cond_cov = tree->rates->_2n2n_vect1;

  For(i,dim) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);

  For(i,dim)
    {
      b = tree->a_edges[i];

      For(j,dim) is_1[j] = 0;
      is_1[b->num]       = 1;

      For(j,dim*dim) cond_cov[j] = 0.0;
      For(j,dim)     cond_mu[j]  = 0.0;
      Normal_Conditional(tree->rates->mean_l,tree->rates->cov_l,a,dim,is_1,1,cond_mu,cond_cov);
      
      tree->rates->cond_var[b->num] = cond_cov[b->num*dim+b->num];
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_All_Reg_Coeff(t_tree *tree)
{
  int i,j;
  short int *is_1;
  phydbl *a;
  int dim;
  t_edge *b;

  dim  = 2*tree->n_otu-3;
  a    = tree->rates->_2n_vect1;
  is_1 = tree->rates->_2n_vect5;
  b    = NULL;

  For(i,dim) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);

  For(i,dim)
    {
      b = tree->a_edges[i];

      For(j,dim) is_1[j] = 0;
      is_1[b->num]       = 1;
      
      Get_Reg_Coeff(tree->rates->mean_l,tree->rates->cov_l,a,dim,is_1,1,tree->rates->reg_coeff+b->num*dim);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


/* Get V(L_{i} | L_{-i}) for all i */
void RATES_Get_Trip_Conditional_Variances(t_tree *tree)
{
  int i,j;
  short int *is_1;
  phydbl *a;
  phydbl *cond_mu, *cond_cov;
  t_node *n;
  int n_otu;

  a        = tree->rates->_2n_vect1;
  is_1     = tree->rates->_2n_vect5;
  cond_mu  = tree->rates->_2n_vect2;
  cond_cov = tree->rates->_2n2n_vect1;
  n        = NULL;
  n_otu    = tree->n_otu;

  For(i,2*n_otu-3) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);

  For(i,2*n_otu-2)
    {
      n = tree->a_nodes[i];
      if(!n->tax)
        {
          For(j,2*n_otu-3) is_1[j] = 0;
          is_1[n->b[0]->num] = 1;
          is_1[n->b[1]->num] = 1;
          is_1[n->b[2]->num] = 1;
          
          Normal_Conditional_Unsorted(tree->rates->mean_l,tree->rates->cov_l,a,2*n_otu-3,is_1,3,cond_mu,cond_cov);
          
          For(j,9) tree->rates->trip_cond_cov[n->num*9+j] = cond_cov[j];
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_All_Trip_Reg_Coeff(t_tree *tree)
{
  int i,j;
  short int *is_1;
  phydbl *a;
  t_node *n;
  int n_otu;

  a     = tree->rates->_2n_vect1;
  is_1  = tree->rates->_2n_vect5;
  n     = NULL;
  n_otu = tree->n_otu;

  For(i,2*n_otu-3) a[i] = tree->rates->mean_l[i] * (Uni() * 0.2 + 0.9);

  For(i,2*n_otu-2)
    {
      n = tree->a_nodes[i];
      if(!n->tax)
        {         
          For(j,2*n_otu-3) is_1[j] = 0;
          is_1[n->b[0]->num] = 1;
          is_1[n->b[1]->num] = 1;
          is_1[n->b[2]->num] = 1;
          
          Get_Reg_Coeff(tree->rates->mean_l,tree->rates->cov_l,a,2*n_otu-3,is_1,3,tree->rates->trip_reg_coeff+n->num*(6*n_otu-9));
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Check_Lk_Rates(t_tree *tree, int *err)
{
  int i;
  phydbl u,u_anc;
  phydbl t,t_anc;
  
  *err = 0;

  For(i,2*tree->n_otu-2)
    {
      u     = tree->rates->br_r[i];
      u_anc = tree->rates->br_r[tree->a_nodes[i]->anc->num];
      t     = tree->rates->nd_t[i];
      t_anc = tree->rates->nd_t[tree->a_nodes[i]->anc->num];

      if(t_anc > t)
        {
          PhyML_Printf("\n. %d %d u=%f u_anc=%f t=%f t_anc=%f",i,tree->a_nodes[i]->anc->num,u,u_anc,t,t_anc);
          PhyML_Printf("\n. %d %d %d",tree->n_root->num,tree->n_root->v[2]->num,tree->n_root->v[1]->num);
          *err = 1;
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Expected_Tree_Length(t_tree *tree)
{
  int n;
  phydbl mean;

  n    = 0;
  mean = 0.0;
  RATES_Expected_Tree_Length_Pre(tree->n_root,tree->n_root->v[2],1.0,&mean,&n,tree);
  RATES_Expected_Tree_Length_Pre(tree->n_root,tree->n_root->v[1],1.0,&mean,&n,tree);

  if(n != 2*tree->n_otu-2)
    {
      PhyML_Printf("\n. n=%d 2n-2=%d",n,2*tree->n_otu-2);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  return mean;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Expected_Tree_Length_Pre(t_node *a, t_node *d, phydbl eranc, phydbl *mean, int *n, t_tree *tree)
{
  phydbl erdes;
  int i;
  phydbl loc_mean;
  int loc_n;
  
  
/*  erdes = u_anc - */
/*     sd*(Dnorm((tree->rates->min_rate-u_anc)/sd,.0,1.) - Dnorm((tree->rates->max_rate-u_anc)/sd,.0,1.))/ */
/*        (Pnorm((tree->rates->max_rate-u_anc)/sd,.0,1.) - Pnorm((tree->rates->min_rate-u_anc)/sd,.0,1.)); */

  /* erdes = Norm_Trunc_Mean(eranc,sd,tree->rates->min_rate,tree->rates->max_rate); */
  erdes = 1.0;

  loc_mean = *mean;
  loc_n = *n;

  loc_mean *= loc_n;
  loc_mean += erdes;
  loc_mean /= (loc_n + 1);

  *mean = loc_mean;
  *n    = loc_n + 1;


  if(d->tax) return;
  else
    For(i,3)
      if(d->v[i] != a && d->b[i] != tree->e_root)
        RATES_Expected_Tree_Length_Pre(d,d->v[i],erdes,mean,n,tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Norm_Fact(t_tree *tree)
{
  int i;
  phydbl r,t,t_anc;
  phydbl num,denom;
  
  num = denom = 0.0;

  For(i,2*tree->n_otu-2)
    {
      r     = tree->rates->br_r[i];
      t     = tree->rates->nd_t[i];
      t_anc = tree->rates->nd_t[tree->a_nodes[i]->anc->num];

      num   += (t-t_anc);
      denom += (t-t_anc) * r;

    }
  tree->rates->norm_fact = num/denom;

  /* !!!!!!!!!!!!!!!!!! */
  tree->rates->norm_fact = 1.0;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Normalise_Rates(t_tree *tree)
{
  phydbl expr,curr;
  int i;
  
  curr = RATES_Average_Substitution_Rate(tree);
  curr /= tree->rates->clock_r;

  /* Set expected mean rate to one such that clock_r is
     easy to interpret regarding the mean values of br_r */
  expr = 1.0;
  
  For(i,2*tree->n_otu-2) 
    {
      tree->rates->br_r[i] *= expr/curr;
     
      if(tree->rates->br_r[i] > tree->rates->max_rate)
        tree->rates->br_r[i] = tree->rates->max_rate;
      
      if(tree->rates->br_r[i] < tree->rates->min_rate)
        tree->rates->br_r[i] = tree->rates->min_rate;
    }
      
  tree->rates->clock_r *= curr/expr;
  /* Branch lengths therefore do not change */

/*   if(tree->rates->clock_r < tree->rates->min_clock) */
/*     { */
/*       PhyML_Printf("\n. Curr mean rates: %G",curr); */
/*       PhyML_Printf("\n. Set clock rate to: %G",tree->rates->clock_r); */
/*       tree->rates->clock_r = tree->rates->min_clock; */
/*     } */
/*   if(tree->rates->clock_r > tree->rates->max_clock) */
/*     { */
/*       PhyML_Printf("\n. Curr mean rates: %G",curr); */
/*       PhyML_Printf("\n. Set clock rate to: %G",tree->rates->clock_r); */
/*       tree->rates->clock_r = tree->rates->max_clock; */
/*     } */

  RATES_Update_Norm_Fact(tree);
  RATES_Update_Cur_Bl(tree);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Find_Max_Dt_Bisec(phydbl r, phydbl r_mean, phydbl ta, phydbl tc, phydbl nu, phydbl threshp, int inf)
{
  phydbl cdfr,cdf0;
  phydbl sd;
  phydbl trunc_cdf;
  phydbl ori_tc, ori_ta;
  phydbl tb;

  ori_tc = tc;
  ori_ta = ta;

/*   PhyML_Printf("\n Max %s r=%f r_mean%f",inf?"inf":"sup",r,r_mean); */
  do
    {
      tb = ta + (tc - ta)/2.;


      sd   = SQRT(nu*(ori_tc - tb));
      cdfr = Pnorm(r ,r_mean,sd);
      cdf0 = Pnorm(.0,r_mean,sd);
      
      if(inf)
        trunc_cdf = (cdfr - cdf0)/(1. - cdf0);
      else
        trunc_cdf = (1. - cdfr)/(1. - cdf0);
        
/*       PhyML_Printf("\n. ta=%15f tb=%15f tc=%15f cdf = %15f",ta,tb,tc,trunc_cdf); */

      if(trunc_cdf > threshp)
        {
          ta = tb;
        }
      else
        {
          tc = tb;
        }

    }while((tc - ta)/(ori_tc - ori_ta) > 0.001);

  if(tb < ori_ta)
    {
      PhyML_Printf("\n. tb < ta r=%f r_mean=%f",r,r_mean);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }
  if(tb > ori_tc)
    {
      PhyML_Printf("\n. tb > tc r=%f r_mean=%f ori_ta=%f ori_tc=%f tb=%f",r,r_mean,ori_ta,ori_tc,tb);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  return tb;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Find_Min_Dt_Bisec(phydbl r, phydbl r_mean, phydbl ta, phydbl tc, phydbl nu, phydbl threshp, int inf)
{
  phydbl cdfr,cdf0;
  phydbl sd;
  phydbl trunc_cdf;
  phydbl ori_tc, ori_ta;
  phydbl tb;

  ori_tc = tc;
  ori_ta = ta;

/*   PhyML_Printf("\n Min %s r=%f r_mean=%f",inf?"inf":"sup",r,r_mean); */
  do
    {

      tb = ta + (tc - ta)/2.;


      sd   = SQRT(nu*(tb - ori_ta));
      cdfr = Pnorm(r ,r_mean,sd);
      cdf0 = Pnorm(.0,r_mean,sd);
      
      if(inf)
        trunc_cdf = (cdfr - cdf0)/(1. - cdf0);
      else
        trunc_cdf = (1. - cdfr)/(1. - cdf0);
        
/*       PhyML_Printf("\n. ta=%15f tb=%15f tc=%15f cdf = %15f",ta,tb,tc,trunc_cdf); */

      if(trunc_cdf > threshp)
        {
          tc = tb;
        }
      else
        {
          ta = tb;
        }

    }while((tc - ta)/(ori_tc - ori_ta) > 0.001);

  if(tb < ori_ta)
    {
      PhyML_Printf("\n. tb < ta r=%f r_mean=%f",r,r_mean);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }
  if(tb > ori_tc)
    {
      PhyML_Printf("\n. tb > tc r=%f r_mean=%f ori_ta=%f ori_tc=%f tb=%f",r,r_mean,ori_ta,ori_tc,tb);
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  return tb;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Min_Max_Interval(phydbl u0, phydbl u1, phydbl u2, phydbl u3, phydbl t0, phydbl t2, phydbl t3, 
                            phydbl *t_min, phydbl *t_max, phydbl nu, phydbl p_thresh, t_tree *tree)
{
  int n_eval;
  int i;
  phydbl sum_cdf;
  phydbl *cdf;
  phydbl t1;
  phydbl mint2t3;
  
  mint2t3 = MIN(t2,t3);
  n_eval = 5;

  cdf = (phydbl *)mCalloc(n_eval,sizeof(phydbl));

  sum_cdf = .0;
  For(i,n_eval)
    {
      t1 = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
      cdf[i] = 
        Dnorm_Trunc(u1,u0,SQRT(nu*(t1-t0)),tree->rates->min_rate,tree->rates->max_rate) *
        Dnorm_Trunc(u2,u1,SQRT(nu*(t2-t1)),tree->rates->min_rate,tree->rates->max_rate) *
        Dnorm_Trunc(u3,u1,SQRT(nu*(t3-t1)),tree->rates->min_rate,tree->rates->max_rate) *
        (mint2t3 - t0) / (n_eval + 1);
      if(i) cdf[i] += cdf[i-1];
    }
  sum_cdf = cdf[i-1];

  For(i,n_eval) cdf[i] /= sum_cdf;

/*   PhyML_Printf("\n"); */
/*   For(i,n_eval) PhyML_Printf("\n* %f %f",cdf[i],sum_cdf); */

  For(i,n_eval) 
    if(cdf[i] > p_thresh)
      {
        *t_min = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
        break;
      }

  For(i,n_eval) 
    if(cdf[i] > (1. - p_thresh))
      {
        *t_max = t0 + (i + 1)*(mint2t3 - t0)/(n_eval + 1);
        break;
      }

  Free(cdf);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Get_Correction_Factor(phydbl mode, phydbl sd, int *err, t_tree *tree)
{

  phydbl K,X0,X1,X2,Y0,Y1,Y2;
  phydbl eps=0.01;
  phydbl A,B;
  phydbl slope,inter;
  phydbl denom;

  *err = NO;


  /* DO NOTHING */
  return 0.0;



  A = tree->rates->min_rate / sd;
  B = tree->rates->max_rate / sd;
  K = mode / sd;

  X0 = 0.;
/*   Y0 = Dnorm(X0-K,.0,1.)/(1. - Pnorm(X0-K,.0,1.)) - X0; */
  Y0 = (Dnorm(A-K+X0,.0,1.)-Dnorm(B-K+X0,.0,1.))/(Pnorm(B-K+X0,.0,1.)-Pnorm(A-K+X0,.0,1.)) - X0;

  X1 = .1;
/*   Y1 = Dnorm(X1-K,.0,1.)/(1. - Pnorm(X1-K,.0,1.)) - X1; */
  Y1 = (Dnorm(A-K+X1,.0,1.)-Dnorm(B-K+X1,.0,1.))/(Pnorm(B-K+X1,.0,1.)-Pnorm(A-K+X1,.0,1.)) - X1;

/*   printf("\n. ^^ mean=%f sd=%f",mode,sd); */

/*   printf("\n. X0=%f Y0=%f X1=%f Y1=%f",X0,Y0,X1,Y1); */

  do
    {
      
      slope = (Y1-Y0)/(X1-X0);
      inter = Y0 - X0*(Y1-Y0)/(X1-X0);

      X2 = -inter/slope;
      
      denom = (Pnorm(B-K+X2,.0,1.)-Pnorm(A-K+X2,.0,1.));

      if(denom < 1.E-10) 
        {
/*        printf("\n. X2 = %f Y2=%f num=%f denom=%f Y1=%f Y0=%f X1=%f X0=%f mode=%f sd=%f",X2,Y2,num,denom,Y1,Y0,X1,X0,mode,sd);           */
          *err = YES;
          break;
        }

/*       Y2 = Dnorm(X2-K,.0,1.)/(1. - Pnorm(X2-K,.0,1.)) - X2; */
      Y2 = (Dnorm(A-K+X2,.0,1.)-Dnorm(B-K+X2,.0,1.))/(Pnorm(B-K+X2,.0,1.)-Pnorm(A-K+X2,.0,1.)) - X2;
      
     /*  printf("\n. X2 = %f Y2=%f num=%f denom=%f Y1=%f Y0=%f X1=%f X0=%f",X2,Y2,num,denom,Y1,Y0,X1,X0); */

      if(X2 > X1)
        {
          X0 = X1;
          X1 = X2;
          Y0 = Y1;
          Y1 = Y2;
        }
      else
        {
          X1 = X0;
          X0 = X2;
          Y1 = Y0;
          Y0 = Y2;
        }      
    }while(fabs(Y2) > eps);

/*   printf("\n. shift = %f X2=%f Y2 = %f",X2*sd,X2,Y2); */

  return X2 * sd;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl Sample_Average_Rate(t_node *a, t_node *d, t_tree *tree)
{
  return(-1.);
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Mean_Br_Len(int iter, t_tree *tree)
{
  int i,dim;
  phydbl *mean;

  if(tree->rates->update_mean_l == NO) return;

  dim = 2*tree->n_otu-3;
  mean = tree->rates->mean_l;

  For(i,dim)
    {      
     mean[i] *= (phydbl)iter;
     mean[i] += tree->a_edges[i]->l->v;
     mean[i] /= (phydbl)(iter+1);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Update_Cov_Br_Len(int iter, t_tree *tree)
{
  int i,j,dim;
  phydbl *mean,*cov;

  if(tree->rates->update_cov_l == NO) return;

  dim = 2*tree->n_otu-3;
  mean = tree->rates->mean_l;
  cov  = tree->rates->cov_l;

  For(i,dim)
    {     
      For(j,dim)
        {
          cov[i*dim+j] += mean[i]*mean[j];
          cov[i*dim+j] *= (phydbl)tree->mcmc->run;
          cov[i*dim+j] += tree->a_edges[i]->l->v * tree->a_edges[j]->l->v;
          cov[i*dim+j] /= (phydbl)(tree->mcmc->run+1);
          cov[i*dim+j] -= mean[i]*mean[j];

          if(i == j && cov[i*dim+j] < MIN_VAR_BL) cov[i*dim+j] = MIN_VAR_BL;
        }      
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Set_Mean_L(t_tree *tree)
{
  int i;
  For(i,2*tree->n_otu-3) 
    {
      tree->rates->mean_l[i] = tree->a_edges[i]->l->v;
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Record_Times(t_tree *tree)
{
  int i;

  if(tree->rates->nd_t_recorded == YES)
    {
      PhyML_Printf("\n. Overwriting recorded times is forbidden.\n");
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  For(i,2*tree->n_otu-1) tree->rates->buff_t[i] = tree->rates->nd_t[i];
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Reset_Times(t_tree *tree)
{
  int i;
  tree->rates->nd_t_recorded = NO;
  For(i,2*tree->n_otu-1) tree->rates->nd_t[i] = tree->rates->buff_t[i];
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Record_Rates(t_tree *tree)
{
  int i;

  if(tree->rates->br_r_recorded == YES)
    {
      PhyML_Printf("\n. Overwriting recorded rates is forbidden.\n");
      PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
      Exit("\n");
    }

  For(i,2*tree->n_otu-2) tree->rates->buff_r[i] = tree->rates->br_r[i];
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Reset_Rates(t_tree *tree)
{
  int i;
  tree->rates->br_r_recorded = NO;
  For(i,2*tree->n_otu-2) tree->rates->br_r[i] = tree->rates->buff_r[i];
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Set_Clock_And_Nu_Max(t_tree *tree)
{
  phydbl dt,nu;
  phydbl min_t;
  int i;
  phydbl step;
  phydbl l_max;
  phydbl max_clock;
  phydbl r_max;
  phydbl tune;
  phydbl pa,pb;

  if(tree->rates->model == THORNE || 
     tree->rates->model == GUINDON)
    {
      tune = 1.05;
      
      if(tree->rates->model_log_rates == NO)
        {
          r_max = LOG(tree->rates->max_rate);
        }
      else
        {
          r_max = tree->rates->max_rate;
        }
      
      l_max = tree->mod->l_max;
      
      min_t = .0;
      For(i,2*tree->n_otu-1) if(tree->rates->t_prior_min[i] < min_t) min_t = tree->rates->t_prior_min[i];
      
      dt = FABS(min_t);
      max_clock = l_max / dt; 
      
      nu   = 1.E-10;
      step = 1.E-1;
      do
        {
          do
            {
              nu += step;
              pa = Dnorm(0.0,  0.0,SQRT(nu*dt)); 
              pb = Dnorm(r_max,0.0,SQRT(nu*dt));
            }while(pa/pb > tune);
          nu -= step;
          step /= 10.;
        }while(step > 1.E-10);
      
      tree->rates->max_nu    = nu;
      /* tree->rates->max_nu    = 1.0; */
      tree->rates->max_clock = max_clock;
      
      PhyML_Printf("\n. Clock rate parameter upper bound set to %f expected subst./site/time unit",tree->rates->max_clock);
      PhyML_Printf("\n. Autocorrelation parameter upper bound set to %f",tree->rates->max_nu);
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

void RATES_Set_Birth_Rate_Boundaries(t_tree *tree)
{
  phydbl lbda;
  phydbl p_above_min,p_below_max;
  phydbl min,max;
  int assign = YES;

  min = -0.01*tree->rates->t_prior_max[tree->n_root->num];
  max = -100.*tree->rates->t_prior_min[tree->n_root->num];

  for(lbda = 0.0001; lbda < 10; lbda+=0.0001)
    {
      p_above_min = 1. - POW(1.-EXP(-lbda*min),tree->n_otu);
      p_below_max = POW(1.-EXP(-lbda*max),tree->n_otu);
 
      if(p_above_min < 1.E-10) 
        { 
          tree->rates->birth_rate_max = lbda;
          break;
        }
      if(p_below_max > 1.E-10 && assign==YES)
        {
          assign = NO;
          tree->rates->birth_rate_min = lbda;
        }
    }
  
  /* tree->rates->birth_rate_min = 1.E-6; */
  /* tree->rates->birth_rate_max = 1.; */
  PhyML_Printf("\n. Birth rate lower bound set to %f.",tree->rates->birth_rate_min);
  PhyML_Printf("\n. Birth rate upper bound set to %f.",tree->rates->birth_rate_max);

}


//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Write_Mean_R_On_Edge_Label(t_node *a, t_node *d, t_edge *b, t_tree *tree)
{

  if(b) 
    {
      if(!b->labels) 
        {
          Make_New_Edge_Label(b);       
          b->n_labels++;
        }
      sprintf(b->labels[0],"%f",tree->rates->mean_r[d->num] / (phydbl)(tree->mcmc->run/tree->mcmc->sample_interval+1.));
    }

  if(d->tax) return;
  else
    {
      int i;
      For(i,3)
        {
          if((d->v[i] != a) && (d->b[i] != tree->e_root))
            {
              RATES_Write_Mean_R_On_Edge_Label(d,d->v[i],d->b[i],tree);
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////



//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


phydbl RATES_Get_Mean_Rate_In_Subtree(t_node *root, t_tree *tree)
{
  phydbl sum;
  int n;

  sum = 0.0;
  n   = 0;

  if(root->tax == NO)
    {
      if(root == tree->n_root)
        {
          RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[2],&sum,&n,tree);
          RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[1],&sum,&n,tree);
        }
      else
        {
          int i;
          For(i,3)
            {
              if(root->v[i] != root->anc && root->b[i] != tree->e_root)
                {
                  RATES_Get_Mean_Rate_In_Subtree_Pre(root,root->v[i],&sum,&n,tree);
                }
            }
        }
      return sum/(phydbl)n;
    }
  else
    {
      return 0.0;
    }
  
  
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Mean_Rate_In_Subtree_Pre(t_node *a, t_node *d, phydbl *sum, int *n, t_tree *tree)
{
  (*sum) += EXP(tree->rates->nd_r[d->num]);
  (*n)   += 1;

  if(d->tax == YES)  return;
  else
    {
      int i;
      For(i,3)
        {
          if(d->v[i] != a && d->b[i] != tree->e_root)
            {
              RATES_Get_Mean_Rate_In_Subtree_Pre(d,d->v[i],sum,n,tree);
            }
        }
    }
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


char *RATES_Get_Model_Name(int model)
{
  char *s;

  s = (char *)mCalloc(T_MAX_NAME,sizeof(char));

  switch(model)
    {
    case GUINDON     : {strcpy(s,"gbs"); break;}
    case THORNE      : {strcpy(s,"gbd"); break;}
    case GAMMA       : {strcpy(s,"gamma"); break;}
    case STRICTCLOCK : {strcpy(s,"clock"); break;}
    default : 
      {
        PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
        Exit("\n");
        break;
     }
    }
  return s;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


void RATES_Get_Survival_Ranks(t_tree *tree)
{
  int i,j;
  phydbl rank;


  For(i,2*tree->n_otu-2)
    {
      rank = 0.;
      For(j,2*tree->n_otu-2)
        {
          if(tree->rates->br_r[i] > tree->rates->br_r[j]) rank += 1.0;
        }

      tree->rates->survival_rank[i] = rank;
    }


}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////