/*

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.

*/

#include "utilities.h"
#include "optimiz.h"
#include "lk.h"
#include "free.h"

/* double UNLIKELY; */
/* double ROUND_MAX; */
/* double MIN_DIFF_LK; */

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

void Optimize_Single_Param_Generic(arbre *tree, double *param, 
				   double start, 
				   double lim_inf, double lim_sup,
				   int n_max_iter)
{
  double ax,bx,cx;
/*   double fa,fb,fc; */
  double lk_init;

  lk_init = tree->tot_loglk;

  tree->mod->s_opt->opt_bl = 
  tree->both_sides         = 0;

  ax =  lim_inf;
  if((*param < lim_inf) ||
     (*param > lim_sup)) bx = (lim_sup-lim_inf)/2.;
  else bx = start;
  cx = lim_sup;
  
/*   Generic_Brak(param, */
/* 	       &ax,&bx,&cx, */
/* 	       &fa,&fb,&fc, */
/* 	       lim_inf, lim_sup, */
/* 	       tree); */
  
  Generic_Brent(param,
		ax,bx,cx,1.e-10,
		param,tree,n_max_iter);

  if(tree->tot_loglk < lk_init-MIN_DIFF_LK) 
    {
      printf("\n %.10f < %.10f --> diff=%.10f\n",
	     tree->tot_loglk,lk_init,
	     tree->tot_loglk-lk_init);
      Exit("\n. Optimisation failed !\n");
    }
}

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

int Generic_Brak(double *param,
		 double *ax, double *bx, double *cx, 
		 double *fa, double *fb, double *fc,
		 double lim_inf, double lim_sup,
		 arbre *tree)
{
   double ulim,u,r,q,fu,dum;

   u = 0.0;
   *param = *ax;

   if(*param > lim_sup) *param = lim_sup;
   if(*param < lim_inf) *param = lim_inf;
   *fa=-Return_Lk(tree);
   *param = *bx;
   if(*param > lim_sup) *param = lim_sup;
   if(*param < lim_inf) *param = lim_inf;
   *fb=-Return_Lk(tree);
   if (*fb > *fa) {
      SHFT(dum,*ax,*bx,dum)
      SHFT(dum,*fb,*fa,dum)
   }
   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
   *param = fabs(*cx);
   if(*param > lim_sup) *param = lim_sup;
   if(*param < lim_inf) *param = lim_inf;
   *fc=-Return_Lk(tree); 
   while (*fb > *fc) 
     {
        
       if(*ax > lim_sup) *ax = lim_sup;
       if(*ax < lim_inf) *ax = lim_inf;
       if(*bx > lim_sup) *bx = lim_sup;
       if(*bx < lim_inf) *bx = lim_inf;
       if(*cx > lim_sup) *cx = lim_sup;
       if(*cx < lim_inf) *cx = lim_inf;
       if(u   > lim_sup) u   = lim_sup;
       if(u   < lim_inf) u   = lim_inf;

       r=(*bx-*ax)*(*fb-*fc);
       q=(*bx-*cx)*(*fb-*fa);
       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
       
       if ((*bx-u)*(u-*cx) > lim_inf) 
	 {
	   *param = fabs(u);
	   if(*param > lim_sup) {*param = u = lim_sup;}
	   if(*param < lim_inf) {*param = u = lim_inf;}
	   fu=-Return_Lk(tree);
	   if (fu < *fc) 
	     {
	       *ax=(*bx);
	       *bx=u;
	       *fa=(*fb);
	       *fb=fu;
	       (*ax)=fabs(*ax);
	       (*bx)=fabs(*bx);
	       (*cx)=fabs(*cx);
	       return(0);
	     } 
	   else if (fu > *fb) 
	     {
	       *cx=u;
	       *fc=fu;	
	       (*ax)=fabs(*ax);
	       (*bx)=fabs(*bx);
	       (*cx)=fabs(*cx);
	       return(0);
	     }
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   *param = fabs(u);
	   if(*param > lim_sup) {*param = u = lim_sup;}
	   if(*param < lim_inf) {*param = u = lim_inf;}
	   fu=-Return_Lk(tree);
	 } 
       else if ((*cx-u)*(u-ulim) > lim_inf) 
	 {
	   *param = fabs(u);
	   if(*param > lim_sup) {*param = u = lim_sup;}
	   if(*param < lim_inf) {*param = u = lim_inf;}
	   fu=-Return_Lk(tree);
	   if (fu < *fc) 
	     {
	       SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
	       *param = fabs(u); 
	       SHFT(*fb,*fc,fu,-Return_Lk(tree))
	     }
	 } 
       else if ((u-ulim)*(ulim-*cx) >= lim_inf) 
	 {
	   u=ulim;
	   *param = fabs(u);
	   if(*param > lim_sup) {*param = u = lim_sup;}
	   if(*param < lim_inf) {*param = u = lim_inf;}
	   fu=-Return_Lk(tree);
	 } 
       else 
	 {
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   *param = fabs(u);
	   if(*param > lim_sup) {*param = u = lim_sup;}
	   if(*param < lim_inf) {*param = u = lim_inf;}
	   fu=-Return_Lk(tree);
	 }
       SHFT(*ax,*bx,*cx,u)
       SHFT(*fa,*fb,*fc,fu)


     }
   (*ax)=fabs(*ax);
   (*bx)=fabs(*bx);
   (*cx)=fabs(*cx);
   return(0);
}

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

double Generic_Brent(double *param, 
		     double ax, double bx, double cx, double tol, 
		     double *xmin, arbre *tree, int n_iter_max)
{
  int iter;
  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
  double e=0.0;
  double init_loglk, max_loglk;
  double bestx;

  d=0.0;
  a=((ax < cx) ? ax : cx);
  b=((ax > cx) ? ax : cx);
  x=w=v=bx;
  *param=bx;
  Lk(tree,tree->data);
  fw=fv=fx=-tree->tot_loglk;
  init_loglk = tree->tot_loglk;
  max_loglk = UNLIKELY;
  bestx = bx;

  for(iter=1;iter<=BRENT_ITMAX;iter++) 
    {
      xm=0.5*(a+b);
      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
      if(fabs(x-xm) <= (tol2-0.5*(b-a))) 
	{
	  if(tree->tot_loglk < init_loglk - MIN_DIFF_LK)
              {
                  printf("\n. WARNING : Brent failed\n");
                  *param = bestx;
                  Lk(tree,tree->data);
              }
	  *xmin=x;
	  return -fx;
	}
      
      if(fabs(e) > tol1) 
	{
	  r=(x-w)*(fx-fv);
	  q=(x-v)*(fx-fw);
	  p=(x-v)*q-(x-w)*r;
	  q=2.0*(q-r);
	  if(q > 0.0) p = -p;
	  q=fabs(q);
	  etemp=e;
	  e=d;
	  if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
              {
                  d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
/*                   printf("Golden section step\n"); */
              }
	  else
              {
                  d=p/q;
                  u=x+d;
                  if (u-a < tol2 || b-u < tol2)
                      d=SIGN(tol1,xm-x);
/*                   printf("Parabolic step\n"); */
              }
        }
      else
          {
              d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
/*               printf("Golden section step (default)\n"); */
          }

      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
      *param=u;
      Lk(tree,tree->data);
      fu=-tree->tot_loglk;
      
      if(tree->tot_loglk > max_loglk)
          {
              max_loglk = tree->tot_loglk;
              bestx = u;
          }

/*       printf("param=%f loglk=%f\n",*param,tree->tot_loglk); */

      if(fu <= fx) 
	{
	  if(iter > n_iter_max) 
	    {
	      if(tree->tot_loglk < init_loglk - MIN_DIFF_LK)
		  printf("\n. WARNING : Brent failed\n");
		  
	      return tree->tot_loglk;
	    }
	  if(u >= x) a=x; else b=x;
	  SHFT(v,w,x,u)
	    SHFT(fv,fw,fx,fu)
	    } 
      else
	{
	  if (u < x) a=u; else b=u;
	  if (fu <= fw || w == x) 
	    {
	      v=w;
	      w=u;
	      fv=fw;
            fw=fu;
	    } 
	  else if (fu <= fv || v == x || v == w) {
            v=u;
            fv=fu;
	  }
	}
    }
  Exit("\n. Too many iterations in BRENT !");
  return(-1);
  /* Not Reached ??  *xmin=x;   */
  /* Not Reached ??  return fx; */
}

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

double RRparam_GTR_Golden(double ax, double bx, double cx, double tol, 
			  double *xmin, arbre *tree, allseq *alldata, double *param, int n_iter_max)
{
   double f1,f2,x0,x1,x2,x3;
   int n_iter;


   x0=ax;
   x3=cx;
   if (fabs(cx-bx) > fabs(bx-ax)) 
     {
       x1=bx;
       x2=bx+GOLDEN_C*(cx-bx);
     } 
   else 
     {
       x2=bx;
       x1=bx-GOLDEN_C*(bx-ax);
     }
   (*param)=x1;

   Lk(tree,alldata);
   f1=-tree->tot_loglk;
   (*param)=x2;

   Lk(tree,alldata);
   f2=-tree->tot_loglk;

   n_iter = 0;
   while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) 
     {

       if (f2 < f1) 
	 {
	   SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
	   (*param)=x2;
	   Lk(tree,alldata);
	   SHFT2(f1,f2,-tree->tot_loglk)
	 } 
       else 
	 {
	   SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
	   (*param)=x1;
	   Lk(tree,alldata);
	   SHFT2(f2,f1,-tree->tot_loglk)
	 }
       
       if(n_iter++ > n_iter_max) break;
       
/*        printf("p=%E %f\n",(*param),tree->tot_loglk); */
     }
   if (f1 < f2) 
    {
       *xmin=x1;
       return f1;
     } 
   else 
     {
       *xmin=x2;
       return f2;
     }
}

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

double Br_Len_Golden(double ax, double bx, double cx, double tol, 
		     double *xmin, edge *b_fcus, arbre *tree)
{
   double f1,f2,x0,x1,x2,x3;

   x0=ax;
   x3=cx;
   if (fabs(cx-bx) > fabs(bx-ax)) 
     {
       x1=bx;
       x2=bx+GOLDEN_C*(cx-bx);
     } 
   else 
     {
       x2=bx;
       x1=bx-GOLDEN_C*(bx-ax);
     }
   
   b_fcus->l=x1;
   f1 = -Lk_At_Given_Edge(tree,b_fcus);
   b_fcus->l=x2;
   f2 = -Lk_At_Given_Edge(tree,b_fcus);
   while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) 
     {
       if (f2 < f1) 
	 {
	   SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
	   b_fcus->l=x2;
	   SHFT2(f1,f2,-Lk_At_Given_Edge(tree,b_fcus))
	 } 
       else 
	 {
	   SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
	   b_fcus->l=x1;
	   SHFT2(f2,f1,-Lk_At_Given_Edge(tree,b_fcus))
	 }
     }
   if (f1 < f2) 
     {
       *xmin=fabs(x1);
       return -f1;
     } 
   else 
     {
       *xmin=fabs(x2);
       return -f2;
     }
}

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

int Br_Len_Brak(double *ax, double *bx, double *cx, 
		double *fa, double *fb, double *fc, 
		edge *b_fcus, arbre *tree)
{
   double ulim,u,r,q,fu,dum;

   b_fcus->l = *ax;
   *fa=-Lk_At_Given_Edge(tree,b_fcus);
   b_fcus->l = *bx;
   *fb=-Lk_At_Given_Edge(tree,b_fcus);
   if (*fb > *fa) {
      SHFT(dum,*ax,*bx,dum)
      SHFT(dum,*fb,*fa,dum)
   }
   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
   b_fcus->l = fabs(*cx);
   *fc=-Lk_At_Given_Edge(tree,b_fcus);
   while (*fb > *fc) 
     {

       r=(*bx-*ax)*(*fb-*fc);
       q=(*bx-*cx)*(*fb-*fa);
       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
       
       if ((*bx-u)*(u-*cx) > 0.0) 
	 {
	   b_fcus->l = fabs(u);
	   fu=-Lk_At_Given_Edge(tree,b_fcus);
	   if (fu < *fc) 
	     {
	       *ax=(*bx);
	       *bx=u;
	       *fa=(*fb);
	       *fb=fu;
	       (*ax)=fabs(*ax);
	       (*bx)=fabs(*bx);
	       (*cx)=fabs(*cx);
	       return(0);
	     } 
	   else if (fu > *fb) 
	     {
	       *cx=u;
	       *fc=fu;	
	       (*ax)=fabs(*ax);
	       (*bx)=fabs(*bx);
	       (*cx)=fabs(*cx);
	       return(0);
	     }
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   b_fcus->l = fabs(u);
	   fu=-Lk_At_Given_Edge(tree,b_fcus);
	 } 
       else if ((*cx-u)*(u-ulim) > 0.0) 
	 {
	   b_fcus->l = fabs(u);
	   fu=-Lk_At_Given_Edge(tree,b_fcus);
	   if (fu < *fc) 
	     {
	       SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
	       b_fcus->l = fabs(u); 
	       SHFT(*fb,*fc,fu,-Lk_At_Given_Edge(tree,b_fcus))
	     }
	 } 
       else if ((u-ulim)*(ulim-*cx) >= 0.0) 
	 {
	   u=ulim;
	   b_fcus->l = fabs(u);
	   fu=-Lk_At_Given_Edge(tree,b_fcus);
	 } 
       else 
	 {
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   b_fcus->l = fabs(u);
	   fu=-Lk_At_Given_Edge(tree,b_fcus);
	 }
       SHFT(*ax,*bx,*cx,u)
       SHFT(*fa,*fb,*fc,fu)
      }
   (*ax)=fabs(*ax);
   (*bx)=fabs(*bx);
   (*cx)=fabs(*cx);
   return(0);
}

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

double Br_Len_Brent(double ax, double bx, double cx, double tol,
		    double *xmin, edge *b_fcus, arbre *tree, int n_iter_max)
{
  int iter;
  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
  double e=0.0;
  
  d=0.0;
  a=((ax < cx) ? ax : cx);
  b=((ax > cx) ? ax : cx);
  x=w=v=bx;
  b_fcus->l = fabs(bx);
  fw=fv=fx=-Lk_At_Given_Edge(tree,b_fcus);
  
  for(iter=1;iter<=BRENT_ITMAX;iter++)
    {
      xm=0.5*(a+b);
      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
      if(fabs(x-xm) <= (tol2-0.5*(b-a)))
	{
	  *xmin=x;
	  Lk_At_Given_Edge(tree,b_fcus);
	  return -fx;
	}
      
      if(fabs(e) > tol1)
	{
	  r=(x-w)*(fx-fv);
	  q=(x-v)*(fx-fw);
	  p=(x-v)*q-(x-w)*r;
	  q=2.0*(q-r);
	  if(q > 0.0) p = -p;
	  q=fabs(q);
	  etemp=e;
	  e=d;
	  if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
	    d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
	  else{
	    d=p/q;
	    u=x+d;
	    if (u-a < tol2 || b-u < tol2)
	      d=SIGN(tol1,xm-x);
	  }
	}
      else
	{
	  d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
	}
      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
      if(u<BL_MIN) u = BL_MIN;
      b_fcus->l=fabs(u);
      fu=-Lk_At_Given_Edge(tree,b_fcus);
/*       printf("edge %d l=%f lnL=%f\n",b_fcus->num,b_fcus->l,fu); */
      if(fu <= fx)
	{
	  if(iter > n_iter_max) 
	    {
	      printf("\n. WARNING : too many iterations in Brent\n");
	      b_fcus->l = fabs(bx);
	      Lk_At_Given_Edge(tree,b_fcus);
	      return tree->tot_loglk;
	    }

	  if(u >= x) a=x; else b=x;
	  SHFT(v,w,x,u)
	    SHFT(fv,fw,fx,fu)
	    }
      else
	{
	  if (u < x) a=u; else b=u;
	  if (fu <= fw || w == x)
	    {
	      v=w;
	      w=u;
	      fv=fw;
            fw=fu;
	    }
	  else if (fu <= fv || v == x || v == w) {
            v=u;
            fv=fu;
	  }
	}
    }
  printf("Too many iterations in BRENT");
  return(-1);
  /* Not Reached ??  *xmin=x;   */
  /* Not Reached ??  return fx; */
}

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

int Dist_Seq_Brak(double *ax, double *bx, double *cx, 
		  double *fa, double *fb, double *fc, 
		  allseq *data, int numseq1, int numseq2, 
		  model *mod)
{
   double ulim,u,r,q,fu,dum;
   double dist;
   double lk,dlk,d2lk;

   dist = *ax;
   *fa=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
   dist = *bx;
   *fb=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
   if (*fb > *fa) {
      SHFT(dum,*ax,*bx,dum)
      SHFT(dum,*fb,*fa,dum)
   }
   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
   dist = fabs(*cx);
   *fc=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
   while (*fb > *fc) 
     {
       r=(*bx-*ax)*(*fb-*fc);
       q=(*bx-*cx)*(*fb-*fa);
       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
       
       if ((*bx-u)*(u-*cx) > 0.0) 
	 {
	   dist = fabs(u);
	   fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
	   if (fu < *fc) 
	     {
	       *ax=(*bx);
	       *bx=u;
	       *fa=(*fb);
	       *fb=fu;
	       return(0);
	     } 
	   else if (fu > *fb) 
	     {
	       *cx=u;
	       *fc=fu;
	       return(0);
	     }
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   dist = fabs(u);
	   fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
	 } 
       else if ((*cx-u)*(u-ulim) > 0.0) 
	 {
	   dist = fabs(u);
	   fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
	   if (fu < *fc) 
	     {
	       SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
	       dist = fabs(u); 
	       SHFT(*fb,*fc,fu,-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk))
	     }
	 } 
       else if ((u-ulim)*(ulim-*cx) >= 0.0) 
	 {
	   u=ulim;
	   dist = fabs(u);
	   fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
	 } 
       else 
	 {
	   u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
	   dist = fabs(u);
	   fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
	 }
       SHFT(*ax,*bx,*cx,u)
       SHFT(*fa,*fb,*fc,fu)
      }
   return(0);
}

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

double Dist_Seq_Brent(double ax, double bx, double cx, double tol, 
		      double *xmin, allseq *data, 
		      int numseq1, int numseq2, model *mod)
{
  int iter;
  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
  double e=0.0;
  double dist;
  double lk,dlk,d2lk;
  
  d=0.0;
  a=((ax < cx) ? ax : cx);
  b=((ax > cx) ? ax : cx);
  x=w=v=bx;
  dist = fabs(bx);
  fw=fv=fx=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
  for(iter=1;iter<=BRENT_ITMAX;iter++) 
    {
      xm=0.5*(a+b);
      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
      if(fabs(x-xm) <= (tol2-0.5*(b-a))) 
	{
	  *xmin=x;
	  return -fx;
	}
      
      if(fabs(e) > tol1) 
	{
	  r=(x-w)*(fx-fv);
	  q=(x-v)*(fx-fw);
	  p=(x-v)*q-(x-w)*r;
	  q=2.0*(q-r);
	  if(q > 0.0) p = -p;
	  q=fabs(q);
	  etemp=e;
	  e=d;
	  if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
	    d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
	  else{
	    d=p/q;
	    u=x+d;
	    if (u-a < tol2 || b-u < tol2)
	      d=SIGN(tol1,xm-x);
	  }
	} else
	  {
	    d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
	  }
      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
      dist=fabs(u);
      fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
      if(fu <= fx) {
	if(u >= x) a=x; else b=x;
	SHFT(v,w,x,u)
	  SHFT(fv,fw,fx,fu)
	  } 
      else
	{
	  if (u < x) a=u; else b=u;
	  if (fu <= fw || w == x) 
	    {
	      v=w;
	      w=u;
	      fv=fw;
	      fw=fu;
	    } 
	  else if (fu <= fv || v == x || v == w) {
            v=u;
            fv=fu;
	  }
	}
    }
  printf("\n . BRENT method failed, trying Newton-Raphson");
  return(+1.0);
  /* Not Reached ??  *xmin=x;   */
  /* Not Reached ??  return fx; */
}

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

double Optimize_Dist(model *mod, double init, allseq *twoseqs)
{
  double d_infa,d_max,d_infb;
  double lk_infa, lk_max, lk_infb, lk;

  d_infa = 100.*BL_MIN;
  d_max  = init;
  d_infb = 3.*init;
  if(init <= BL_MIN) {d_infa = -BL_START; d_max = .0; d_infb = BL_START;}
  lk_infa = lk_max = lk_infb = .0;

  Dist_Seq_Brak(&d_infa, &d_max, &d_infb,
		&lk_infa,&lk_max,&lk_infb,
		twoseqs,0,1,mod);

  lk = (double)Dist_Seq_Brent(d_infa,d_max,d_infb,
			      1.e-5,&d_max,twoseqs,0,1,mod);
  if(lk > .0) return -1.0;
  else        return d_max;

}

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

void Round_Optimize(arbre *tree, allseq *data)
{
  int n_round,each;
  double lk_old, lk_new;
  node *root;

  lk_new = tree->tot_loglk;
  lk_old = UNLIKELY;
  n_round = 0;
  each = 1;
  root = tree->noeud[0];


  tree->mod->s_opt->opt_bl = 0;
  tree->both_sides         = 1;
  Lk(tree,data);


  while(n_round < ROUND_MAX)
    {

        (!((n_round+2)%2))?(root=tree->noeud[0]):(root=tree->noeud[tree->n_otu-1]);
        Optimize_Br_Len_Serie(root,
                              root->v[0],
                              root->b[0],
                              tree,
                              data,
                              5);
        

        tree->mod->s_opt->opt_bl = 0;
        tree->both_sides         = 1;
        Lk(tree,data);
        
        if(!each)
            {
                each = 1;
                if(tree->mod->s_opt->print) printf("\n");
                Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
                tree->mod->s_opt->opt_bl = 0;
                tree->both_sides         = 1;
                Lk(tree,data);
            }
        
        lk_new = tree->tot_loglk;
        
        if(tree->mod->s_opt->print)
            {
                if(lk_old < UNLIKELY+1)
                    printf("\n. Log(lk) :               * -> %15.6f ",lk_new);
                else
                    printf("\n. Log(lk) : %15.6f -> %15.6f ",lk_old,lk_new);
                fflush(NULL);
            }
        
        if(lk_new < lk_old - MIN_DIFF_LK) Exit("\n. Optimisation failed ! (Round_Optimize)\n");
        if(fabs(lk_new - lk_old) < MIN_DIFF_LK)  break;
        else lk_old  = lk_new;
        n_round++;
        each--;
    }
  
  if(tree->mod->s_opt->print) printf("\n");
  Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
}

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

void Optimize_Br_Len_Serie(node *a, node *d, edge *b_fcus, 
			   arbre *tree,allseq *alldata, int n_passes)
{
  int i;
  double l_infa,l_max,l_infb;
  double lk_init;
  
  lk_init = tree->tot_loglk;

  l_infa = 10.*b_fcus->l;
  l_max  = b_fcus->l;
  l_infb = BL_MIN;
 

  Br_Len_Brent(l_infa,l_max,l_infb,
               1.e-5,
               &(b_fcus->l),
               b_fcus,tree,1000);

  /* Golden method is generally slower than Brent */
  /*   Br_Len_Golden(l_infa,l_max,l_infb, */
  /* 		1.e-5, */
  /* 		&(b_fcus->l), */
  /* 		b_fcus,tree); */
/*   printf("Edge %d -> %20f\n",b_fcus->num,tree->tot_loglk); */
  if(tree->tot_loglk < lk_init - MIN_DIFF_LK)
    {
      printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
      printf("%f -- %f \n",lk_init,tree->tot_loglk);
      Exit("\n. Err. in Optimize_Br_Len_Serie\n");
    }

/*   printf("Edge %3d -> %f %f\n", */
/* 	 b_fcus->num, */
/* 	 tree->tot_loglk, */
/* 	 b_fcus->l); fflush(NULL); */
    
  
  if(d->tax) return;
  else For(i,3) if(d->v[i] != a)
    {
      Update_P_Lk(tree,d->b[i],d);
      Optimize_Br_Len_Serie(d,d->v[i],d->b[i],tree,alldata,n_passes);
    }
  For(i,3) if((d->v[i] == a) && !(d->v[i]->tax)) Update_P_Lk(tree,d->b[i],d);
}

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

void Moving_Backward(arbre *tree, allseq *alldata, 
		     double **step, double lk_base, double *l_branch_init)
{
  double lk_new;
  int i,n_iter_max;


  tree->mod->s_opt->opt_bl = tree->both_sides = 0;
  n_iter_max = 10;
  lk_new = lk_base-MIN_DIFF_LK-1.;
  

  For(i,2*tree->n_otu-3) 
    {
      tree->t_edges[i]->l = l_branch_init[i];
      tree->tot_dloglk[i] = .0;
      tree->tot_d2loglk[i] = -.1;
    }

  while(lk_new < lk_base-MIN_DIFF_LK)
    {
      if(!n_iter_max) 
	{
	  For(i,2*tree->n_otu-3) tree->t_edges[i]->l = l_branch_init[i];
	  tree->mod->s_opt->opt_bl = tree->both_sides = 1;
	  return;
	}

      For(i,2*tree->n_otu-3)
	{
	  (*step)[i] *= .5;
	  tree->t_edges[i]->l += (*step)[i];
	  if(tree->t_edges[i]->l < BL_MIN) tree->t_edges[i]->l = BL_MIN;
	}

      Lk(tree,alldata);
      lk_new = tree->tot_loglk;
      n_iter_max--;
    }
  tree->mod->s_opt->opt_bl = tree->both_sides = 1;
}


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

double Optimize_One_Dist(allseq *data, int numseq1, int numseq2, double init_dist, model *mod)
{

  int n_iter, n_iter_mov;
  double lk_new,lk_old,lk_best;
  double step;
  double dlk,d2lk;
  double dist_new,dist_best;

  step = .0;
  dlk  = d2lk = .0;
  lk_old = lk_new = lk_best = .0;
  dist_new = dist_best = init_dist;
  
  n_iter = 0;
  for(;;)
    {
      
      if(dist_new < BL_MIN) dist_new = BL_MIN;

      Lk_Given_Two_Seq(data,numseq1,numseq2,dist_new,mod,&lk_new,&dlk,&d2lk);

      if(((fabs(lk_new - lk_old) < MIN_DIFF_LK) && (lk_new >= lk_old)))
	break;
      else
	{
	  if(((fabs(lk_new - lk_old) > MIN_DIFF_LK) 
	      && (lk_new > lk_old)) 
	      || (!n_iter))
	    {
	      dist_best = dist_new;
	      lk_best   = lk_new;
	      lk_old    = lk_new;
	    }
	  else
	    {
	      n_iter_mov = 20;
	      do
		{
		  dist_new = dist_best;
		  step /= 2.;
		  dist_new += step;
		  Lk_Given_Two_Seq(data,numseq1,numseq2,
				   dist_new,mod,&lk_new,&dlk,&d2lk);
		  n_iter_mov--;
		  if(!n_iter_mov) return dist_best;
		}while(lk_new < lk_best);
	      dlk = .0;
	    }
	}
      
      step = -dlk/d2lk;

      if(d2lk > 0.0) 
	{
	  step = .0;
	  dist_new /= 1.5;
	}
      
      while(fabs(step) > dist_new) step /= 1.5;
      dist_new += step;

      
      n_iter++;
      if(n_iter > 50) break;
    }
  return dist_new;
}

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

int Count_Swap(arbre *tree)
{
  int i;
  
  tree->n_swap = 0;
  For(i,2*tree->n_otu-3)
    {
      if((!tree->t_edges[i]->left->tax) &&
	 (!tree->t_edges[i]->rght->tax))
	{
	  if(tree->t_edges[i]->diff_lk > -2.0)
	    {
	      tree->n_swap++;
	    }
	}
    }
  return tree->n_swap;
}

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

void Optimiz_Ext_Br(arbre *tree)
{
  int i;
  edge *b;
  double l_infa,l_max,l_infb,l_init;
  double lk_init;
  
  lk_init = tree->tot_loglk;
  
  For(i,2*tree->n_otu-3)
    {
      b = tree->t_edges[i];
      if((b->left->tax) || (b->rght->tax))
	{

	  l_init = b->l;
  
	  l_infa = 100.;
	  l_max  = b->l;
	  l_infb = -10.;
	  
	  Br_Len_Brent(l_infa,l_max,l_infb, // Note: kept call here (might have wanted side-effects)
                       1.e-5,
                       &(b->l),
                       b,tree,1000);

	  b->ql[0] = b->l;
	  b->best_conf = 1;
	  b->l = l_init;
	}
    }
  tree->tot_loglk = lk_init; 
}

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

void Optimiz_All_Free_Param(arbre *tree, int verbose)
{
  int  init_both_sides, init_derivatives;
  

  init_both_sides     = tree->both_sides;
  init_derivatives    = tree->mod->s_opt->opt_bl;
  tree->both_sides    = 0;
  tree->mod->s_opt->opt_bl = 0;

  
  if((tree->mod->whichmodel == 7) ||
     ((tree->mod->whichmodel == 8) && 
      (tree->mod->s_opt->opt_rr_param) && 
      (tree->mod->n_diff_rr_param > 1)))
    {
        int failed;
        
        failed = 0;

        if(verbose) 
            {
                (tree->mod->whichmodel == 7)?
                    (printf("\n. Optimisation of the GTR parameters...\n")):	
                    (printf("\n. Optimisation of the custom model parameters...\n"));
            }
        
        tree->mod->update_eigen = 1;
        BFGS(tree,tree->mod->rr_param_values,tree->mod->n_diff_rr_param,1.e-5,1.e-7,
             &Return_Abs_Lk,
             &Num_Derivative_Several_Param ,
             &Lnsrch_RR_Param,&failed);
        
        if(failed)
            {
                int i;
                
                printf("\n. Optimising one-by-one...\n");
                
                For(i,tree->mod->n_diff_rr_param) 
                    if(i != 5)
                        Optimize_Single_Param_Generic(tree,&(tree->mod->rr_param_values[i]),tree->mod->rr_param_values[i],1.E-20,1.E+10,1000);
            }
        
        tree->mod->update_eigen = 0;
      
    }

  if(tree->mod->s_opt->opt_kappa)
    {
      if(verbose) printf("\n. Optimisation of the ts/tv ratio...\n");fflush(stdout);   
      Optimize_Single_Param_Generic(tree,&(tree->mod->kappa),tree->mod->kappa,0.1,100,100);
      /* 	      printf("kappa = %f\n",tree->mod->kappa); */
      
    }
  
  if(tree->mod->s_opt->opt_lambda) 
    {
      Optimize_Single_Param_Generic(tree,&(tree->mod->lambda),tree->mod->lambda,0.001,100,50);
      /* 	      printf("lambda = %f\n",tree->mod->lambda); */
    }
  
  if(tree->mod->s_opt->opt_pinvar) 
    {
      if(verbose) printf("\n. Optimisation of the proportion of invariable sites...\n");fflush(stdout);    
      tree->mod->pinvar = 0.5;
      Optimize_Single_Param_Generic(tree,&(tree->mod->pinvar),tree->mod->pinvar,0.0001,0.9999,100);            
/*       printf("p-invar = %f\n",tree->mod->pinvar); */
    }

  if(tree->mod->s_opt->opt_alpha) 
    { 
      if(verbose) printf("\n. Optimisation of the gamma shape parameter...\n");fflush(stdout);
      Optimize_Single_Param_Generic(tree,&(tree->mod->alpha),tree->mod->alpha,0.01,100,100);
/*       printf("alpha = %f %f\n",tree->mod->alpha,Return_Lk(tree));       */
    }

  if(tree->mod->s_opt->opt_bfreq)
    {
        int failed,i;
        
        failed = 0;
        tree->mod->update_eigen = 1;
        if(verbose) printf("\n. Optimisation of nucleotide frequencies...\n");
        BFGS(tree,tree->mod->pi,4,1.e-5,1.e-7,&Return_Abs_Lk,&Num_Derivative_Several_Param,&Lnsrch_Nucleotide_Frequencies,&failed);

        if(failed)
            {
                For(i,5) 
                    Optimize_Single_Param_Generic(tree,&(tree->mod->pi[i]),tree->mod->pi[i],1.E-10,0.999999,1000);                
            }
        tree->mod->update_eigen = 0;
    }


  tree->both_sides    = init_both_sides;
  tree->mod->s_opt->opt_bl = init_derivatives;

}



#define ITMAX 200
#define EPS   3.0e-8
#define TOLX (4*EPS)
#define STPMX 100.0
static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

void BFGS(arbre *tree, double *p, int n, double gtol, double step_size,
	  double(*func)(), void (*dfunc)(), void (*lnsrch)(),int *failed)
{

  int check,i,its,j;
  double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test,fret;
  double *dg,*g,*hdg,**hessin,*pnew,*xi;
  
  hessin = (double **)mCalloc(n,sizeof(double *));
  For(i,n) hessin[i] = (double *)mCalloc(n,sizeof(double));
  dg   = (double *)mCalloc(n,sizeof(double ));
  g    = (double *)mCalloc(n,sizeof(double ));
  pnew = (double *)mCalloc(n,sizeof(double ));
  hdg  = (double *)mCalloc(n,sizeof(double ));
  xi   = (double *)mCalloc(n,sizeof(double ));
  

  fp=(*func)(tree);
  (*dfunc)(tree,p,n,step_size,func,g);


  for (i=0;i<n;i++) 
    {
      for (j=0;j<n;j++) hessin[i][j]=0.0;
      hessin[i][i]=1.0;
      xi[i] = -g[i];
      sum += p[i]*p[i];
    }

  stpmax=STPMX*MAX(sqrt(sum),(double)n);

  for(its=1;its<=ITMAX;its++) 
    {
      lnsrch(tree,n,p,fp,g,xi,pnew,&fret,stpmax,&check);

/*       printf("BFGS -> %f\n",tree->tot_loglk); */

      fp = fret;
      
      for (i=0;i<n;i++) 
	{
	  xi[i]=pnew[i]-p[i];
	  p[i]=pnew[i];
	}

      test=0.0;
      for (i=0;i<n;i++) 
	{
	  temp=fabs(xi[i])/MAX(fabs(p[i]),1.0);
	  if (temp > test) test=temp;
	}
      if (test < TOLX) 
	{
	  (*func)(tree);
	  For(i,n) Free(hessin[i]);
	  free(hessin);
	  free(xi);
	  free(pnew);
	  free(hdg);
	  free(g);
	  free(dg);   

	  if(its == 1) 
              {
                  printf("\n. WARNING : BFGS failed ! \n");
                  *failed = 1;
              }
	  return;
	}

      for (i=0;i<n;i++) dg[i]=g[i];

      (*dfunc)(tree,p,n,step_size,func,g);

      test=0.0;
      den=MAX(fret,1.0);
      for (i=0;i<n;i++) 
	{
	  temp=fabs(g[i])*MAX(fabs(p[i]),1.0)/den;
	  if (temp > test) test=temp;
	}
      if (test < gtol) 
	{
	  (*func)(tree);
	  For(i,n) Free(hessin[i]);
	  free(hessin);
	  free(xi);
	  free(pnew);
	  free(hdg);
	  free(g);
	  free(dg);   
	  return;
	}

    for (i=0;i<n;i++) dg[i]=g[i]-dg[i];

    for (i=0;i<n;i++) 
      {
	hdg[i]=0.0;
	for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
      }

    fac=fae=sumdg=sumxi=0.0;
    for (i=0;i<n;i++) 
      {
	fac += dg[i]*xi[i];
	fae += dg[i]*hdg[i];
	sumdg += SQR(dg[i]);
	sumxi += SQR(xi[i]);
      }
    
    if(fac*fac > EPS*sumdg*sumxi) 
      {
	fac=1.0/fac;
	fad=1.0/fae;
	for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
	for (i=0;i<n;i++) 
	  {
	    for (j=0;j<n;j++) 
	      {
		hessin[i][j] += fac*xi[i]*xi[j]
		  -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
	      }
	  }
      }
    for (i=0;i<n;i++) 
      {
	xi[i]=0.0;
	for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
      }
    }
  Exit("\n. Too many iterations in BFGS...\n");
  For(i,n) Free(hessin[i]);
  free(hessin);
  free(xi);
  free(pnew);
  free(hdg);
  free(g);
  free(dg);   
}

#undef ITMAX
#undef EPS
#undef TOLX
#undef STPMX

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


#define ALF 1.0e-4
#define TOLX 1.0e-7

void Lnsrch_RR_Param(arbre *tree, int n, double *xold, double fold, 
		     double *g, double *p, double *x,
		     double *f, double stpmax, int *check)
{
  int i;
  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
  double *local_xold;

  alam = alam2 = f2 = fold2 = tmplam = .0;

  local_xold = (double *)mCalloc(n,sizeof(double));
  For(i,n) local_xold[i] = xold[i];


  *check=0;
  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
  sum=sqrt(sum);
  if(sum > stpmax)
    for(i=0;i<n;i++) p[i] *= stpmax/sum;
  for(slope=0.0,i=0;i<n;i++)
    slope += g[i]*p[i];
  test=0.0;
  for(i=0;i<n;i++) 
    {
      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
      if (temp > test) test=temp;
    }
  alamin=TOLX/test;
  alam=1.0;
  for (;;) 
    {
      for(i=0;i<n;i++) 
	{
	  x[i]=local_xold[i]+alam*p[i];
	}

      /**/
      for(i=0;i<n;i++)
	{
	  tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
	  if(tree->mod->rr_param_values[i] < 0.0) break;
	}
      /**/

      if(i==n) 
	{
	  *f=Return_Abs_Lk(tree);
/* 	  printf("loglk = %f\n",*f); */
	}
      else     *f=1.+fold+ALF*alam*slope;
      if (alam < alamin)
	{
	  for (i=0;i<n;i++) 
	    {
	      x[i]=local_xold[i];
	      if(x[i] < .0) break;
	    }
	  /**/      
	  for(i=0;i<n;i++)
	    {
	      tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
	      if(tree->mod->rr_param_values[i] < 0.0) 
		tree->mod->rr_param_values[i] = 0.0;
	    }
	  /**/

	  *check=1;
	  For(i,n) xold[i] = local_xold[i];
	  Free(local_xold);
	  return;
	} 
      else if (*f <= fold+ALF*alam*slope) 
	{
	  For(i,n) xold[i] = local_xold[i];
	  Free(local_xold); 
	  return;
	}
      else 
	{
	  if (alam == 1.0)
	    tmplam = -slope/(2.0*(*f-fold-slope));
	  else 
	    {
	      rhs1 = *f-fold-alam*slope;
	      rhs2=f2-fold2-alam2*slope;
	      a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
	      b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
	      if (a == 0.0) tmplam = -slope/(2.0*b);
	      else 
		{
		  disc=b*b-3.0*a*slope;
		  if (disc<0.0) tmplam = 0.5*alam;
		  else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
		  else tmplam = -slope/(b+sqrt(disc));
		}
	      if (tmplam>0.5*alam) tmplam=0.5*alam;
	    }
	}
      alam2=alam;
      f2 = *f;
      fold2=fold;
      alam=MAX(tmplam,0.1*alam);
    }
  Free(local_xold);
}

#undef ALF
#undef TOLX
#undef NRANSI

/*********************************************************/
#define ALF 1.0e-4
#define TOLX 1.0e-7

void Lnsrch_Nucleotide_Frequencies(arbre *tree, int n, double *xold, double fold, double *g, double *p, double *x,
				   double *f, double stpmax, int *check)
{
  int i;
  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
  double *local_xold;

  alam = alam2 = f2 = fold2 = tmplam = .0;

  local_xold = (double *)mCalloc(n,sizeof(double));
  For(i,n) local_xold[i] = xold[i];


  *check=0;
  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
  sum=sqrt(sum);
  if(sum > stpmax)
    for(i=0;i<n;i++) p[i] *= stpmax/sum;
  for(slope=0.0,i=0;i<n;i++)
    slope += g[i]*p[i];
  test=0.0;
  for(i=0;i<n;i++) 
    {
      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
      if (temp > test) test=temp;
    }
  alamin=TOLX/test;
  alam=1.0;
  for (;;) 
    {
      for(i=0;i<n;i++) x[i]=fabs(local_xold[i]+alam*p[i]);
      /**/      
      for(i=0;i<n;i++) 
	{
	  tree->mod->pi[i]=fabs(local_xold[i]+alam*p[i]);
/* 	  if( */
/* 	     (tree->mod->pi[i] < 0.001) || */
/* 	     (tree->mod->pi[i] > 0.999) */
/* 	     ) */
/* 	    break; */
	}
      /**/
      if(i==n) 
	{
	  *f=Return_Abs_Lk(tree);
	}
      else     *f=1.+fold+ALF*alam*slope;
      if (alam < alamin)
	{
	  for (i=0;i<n;i++) x[i]=local_xold[i];
	  for (i=0;i<n;i++) tree->mod->pi[i]=local_xold[i];
	  *check=1;
	  For(i,n) xold[i] = local_xold[i];
	  Free(local_xold);
	  return;
	} 
      else if (*f <= fold+ALF*alam*slope) 
	{
	  For(i,n) xold[i] = local_xold[i];
	  Free(local_xold); 
	  return;
	}
      else 
	{
	  if (alam == 1.0)
	    tmplam = -slope/(2.0*(*f-fold-slope));
	  else 
	    {
	      rhs1 = *f-fold-alam*slope;
	      rhs2=f2-fold2-alam2*slope;
	      a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
	      b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
	      if (a == 0.0) tmplam = -slope/(2.0*b);
	      else 
		{
		  disc=b*b-3.0*a*slope;
		  if (disc<0.0) 
		    {
		      disc=b*b-3.0*a*slope;
		      if (disc<0.0) tmplam = 0.5*alam;
		      else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
		      else tmplam = -slope/(b+sqrt(disc));
		    }
		  else tmplam=(-b+sqrt(disc))/(3.0*a);
		}
	      if (tmplam>0.5*alam) tmplam=0.5*alam;
	    }
	}
      alam2=alam;
      f2 = *f;
      fold2=fold;
      alam=MAX(tmplam,0.1*alam);
    }
  Free(local_xold);
}

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

/* void Optimize_Global_Rate(arbre *tree) */
/* { */
/*     printf("\n. Global rate (%f->)",tree->tot_loglk); */
/*     Optimize_Single_Param_Generic(tree,&(tree->tbl),tree->tbl,BL_MIN,1.E+4,100); */
/*     printf("%f)\n",tree->tot_loglk); */
/* } */


#undef ALF
#undef TOLX
#undef NRANSI