source: branches/profile/GDE/PHYML/optimiz.c

Last change on this file was 4073, checked in by westram, 13 years ago
  • phyml 2.4.5
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.3 KB
Line 
1/*
2
3PHYML :  a program that  computes maximum likelihood  phylogenies from
4DNA or AA homologous sequences
5
6Copyright (C) Stephane Guindon. Oct 2003 onward
7
8All parts of  the source except where indicated  are distributed under
9the GNU public licence.  See http://www.opensource.org for details.
10
11*/
12
13#include "utilities.h"
14#include "optimiz.h"
15#include "lk.h"
16#include "free.h"
17
18/* double UNLIKELY; */
19/* double ROUND_MAX; */
20/* double MIN_DIFF_LK; */
21
22/*********************************************************/
23
24void Optimize_Single_Param_Generic(arbre *tree, double *param, 
25                                   double start, 
26                                   double lim_inf, double lim_sup,
27                                   int n_max_iter)
28{
29  double ax,bx,cx;
30/*   double fa,fb,fc; */
31  double lk_init;
32
33  lk_init = tree->tot_loglk;
34
35  tree->mod->s_opt->opt_bl = 
36  tree->both_sides         = 0;
37
38  ax =  lim_inf;
39  if((*param < lim_inf) ||
40     (*param > lim_sup)) bx = (lim_sup-lim_inf)/2.;
41  else bx = start;
42  cx = lim_sup;
43 
44/*   Generic_Brak(param, */
45/*             &ax,&bx,&cx, */
46/*             &fa,&fb,&fc, */
47/*             lim_inf, lim_sup, */
48/*             tree); */
49 
50  Generic_Brent(param,
51                ax,bx,cx,1.e-10,
52                param,tree,n_max_iter);
53
54  if(tree->tot_loglk < lk_init-MIN_DIFF_LK) 
55    {
56      printf("\n %.10f < %.10f --> diff=%.10f\n",
57             tree->tot_loglk,lk_init,
58             tree->tot_loglk-lk_init);
59      Exit("\n. Optimisation failed !\n");
60    }
61}
62
63/*********************************************************/
64
65int Generic_Brak(double *param,
66                 double *ax, double *bx, double *cx, 
67                 double *fa, double *fb, double *fc,
68                 double lim_inf, double lim_sup,
69                 arbre *tree)
70{
71   double ulim,u,r,q,fu,dum;
72
73   u = 0.0;
74   *param = *ax;
75
76   if(*param > lim_sup) *param = lim_sup;
77   if(*param < lim_inf) *param = lim_inf;
78   *fa=-Return_Lk(tree);
79   *param = *bx;
80   if(*param > lim_sup) *param = lim_sup;
81   if(*param < lim_inf) *param = lim_inf;
82   *fb=-Return_Lk(tree);
83   if (*fb > *fa) {
84      SHFT(dum,*ax,*bx,dum)
85      SHFT(dum,*fb,*fa,dum)
86   }
87   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
88   *param = fabs(*cx);
89   if(*param > lim_sup) *param = lim_sup;
90   if(*param < lim_inf) *param = lim_inf;
91   *fc=-Return_Lk(tree); 
92   while (*fb > *fc) 
93     {
94       
95       if(*ax > lim_sup) *ax = lim_sup;
96       if(*ax < lim_inf) *ax = lim_inf;
97       if(*bx > lim_sup) *bx = lim_sup;
98       if(*bx < lim_inf) *bx = lim_inf;
99       if(*cx > lim_sup) *cx = lim_sup;
100       if(*cx < lim_inf) *cx = lim_inf;
101       if(u   > lim_sup) u   = lim_sup;
102       if(u   < lim_inf) u   = lim_inf;
103
104       r=(*bx-*ax)*(*fb-*fc);
105       q=(*bx-*cx)*(*fb-*fa);
106       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
107               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
108       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
109       
110       if ((*bx-u)*(u-*cx) > lim_inf) 
111         {
112           *param = fabs(u);
113           if(*param > lim_sup) {*param = u = lim_sup;}
114           if(*param < lim_inf) {*param = u = lim_inf;}
115           fu=-Return_Lk(tree);
116           if (fu < *fc) 
117             {
118               *ax=(*bx);
119               *bx=u;
120               *fa=(*fb);
121               *fb=fu;
122               (*ax)=fabs(*ax);
123               (*bx)=fabs(*bx);
124               (*cx)=fabs(*cx);
125               return(0);
126             } 
127           else if (fu > *fb) 
128             {
129               *cx=u;
130               *fc=fu; 
131               (*ax)=fabs(*ax);
132               (*bx)=fabs(*bx);
133               (*cx)=fabs(*cx);
134               return(0);
135             }
136           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
137           *param = fabs(u);
138           if(*param > lim_sup) {*param = u = lim_sup;}
139           if(*param < lim_inf) {*param = u = lim_inf;}
140           fu=-Return_Lk(tree);
141         } 
142       else if ((*cx-u)*(u-ulim) > lim_inf) 
143         {
144           *param = fabs(u);
145           if(*param > lim_sup) {*param = u = lim_sup;}
146           if(*param < lim_inf) {*param = u = lim_inf;}
147           fu=-Return_Lk(tree);
148           if (fu < *fc) 
149             {
150               SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
151               *param = fabs(u); 
152               SHFT(*fb,*fc,fu,-Return_Lk(tree))
153             }
154         } 
155       else if ((u-ulim)*(ulim-*cx) >= lim_inf) 
156         {
157           u=ulim;
158           *param = fabs(u);
159           if(*param > lim_sup) {*param = u = lim_sup;}
160           if(*param < lim_inf) {*param = u = lim_inf;}
161           fu=-Return_Lk(tree);
162         } 
163       else 
164         {
165           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
166           *param = fabs(u);
167           if(*param > lim_sup) {*param = u = lim_sup;}
168           if(*param < lim_inf) {*param = u = lim_inf;}
169           fu=-Return_Lk(tree);
170         }
171       SHFT(*ax,*bx,*cx,u)
172       SHFT(*fa,*fb,*fc,fu)
173
174
175     }
176   (*ax)=fabs(*ax);
177   (*bx)=fabs(*bx);
178   (*cx)=fabs(*cx);
179   return(0);
180}
181
182/*********************************************************/
183
184double Generic_Brent(double *param, 
185                     double ax, double bx, double cx, double tol, 
186                     double *xmin, arbre *tree, int n_iter_max)
187{
188  int iter;
189  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
190  double e=0.0;
191  double init_loglk, max_loglk;
192  double bestx;
193
194  d=0.0;
195  a=((ax < cx) ? ax : cx);
196  b=((ax > cx) ? ax : cx);
197  x=w=v=bx;
198  *param=bx;
199  Lk(tree,tree->data);
200  fw=fv=fx=-tree->tot_loglk;
201  init_loglk = tree->tot_loglk;
202  max_loglk = UNLIKELY;
203  bestx = bx;
204
205  for(iter=1;iter<=BRENT_ITMAX;iter++) 
206    {
207      xm=0.5*(a+b);
208      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
209      if(fabs(x-xm) <= (tol2-0.5*(b-a))) 
210        {
211          if(tree->tot_loglk < init_loglk - MIN_DIFF_LK)
212              {
213                  printf("\n. WARNING : Brent failed\n");
214                  *param = bestx;
215                  Lk(tree,tree->data);
216              }
217          *xmin=x;
218          return -fx;
219        }
220     
221      if(fabs(e) > tol1) 
222        {
223          r=(x-w)*(fx-fv);
224          q=(x-v)*(fx-fw);
225          p=(x-v)*q-(x-w)*r;
226          q=2.0*(q-r);
227          if(q > 0.0) p = -p;
228          q=fabs(q);
229          etemp=e;
230          e=d;
231          if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
232              {
233                  d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
234/*                   printf("Golden section step\n"); */
235              }
236          else
237              {
238                  d=p/q;
239                  u=x+d;
240                  if (u-a < tol2 || b-u < tol2)
241                      d=SIGN(tol1,xm-x);
242/*                   printf("Parabolic step\n"); */
243              }
244        }
245      else
246          {
247              d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
248/*               printf("Golden section step (default)\n"); */
249          }
250
251      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
252      *param=u;
253      Lk(tree,tree->data);
254      fu=-tree->tot_loglk;
255     
256      if(tree->tot_loglk > max_loglk)
257          {
258              max_loglk = tree->tot_loglk;
259              bestx = u;
260          }
261
262/*       printf("param=%f loglk=%f\n",*param,tree->tot_loglk); */
263
264      if(fu <= fx) 
265        {
266          if(iter > n_iter_max) 
267            {
268              if(tree->tot_loglk < init_loglk - MIN_DIFF_LK)
269                  printf("\n. WARNING : Brent failed\n");
270                 
271              return tree->tot_loglk;
272            }
273          if(u >= x) a=x; else b=x;
274          SHFT(v,w,x,u)
275            SHFT(fv,fw,fx,fu)
276            } 
277      else
278        {
279          if (u < x) a=u; else b=u;
280          if (fu <= fw || w == x) 
281            {
282              v=w;
283              w=u;
284              fv=fw;
285            fw=fu;
286            } 
287          else if (fu <= fv || v == x || v == w) {
288            v=u;
289            fv=fu;
290          }
291        }
292    }
293  Exit("\n. Too many iterations in BRENT !");
294  return(-1);
295  /* Not Reached ??  *xmin=x;   */
296  /* Not Reached ??  return fx; */
297}
298
299/*********************************************************/
300
301double RRparam_GTR_Golden(double ax, double bx, double cx, double tol, 
302                          double *xmin, arbre *tree, allseq *alldata, double *param, int n_iter_max)
303{
304   double f1,f2,x0,x1,x2,x3;
305   int n_iter;
306
307
308   x0=ax;
309   x3=cx;
310   if (fabs(cx-bx) > fabs(bx-ax)) 
311     {
312       x1=bx;
313       x2=bx+GOLDEN_C*(cx-bx);
314     } 
315   else 
316     {
317       x2=bx;
318       x1=bx-GOLDEN_C*(bx-ax);
319     }
320   (*param)=x1;
321
322   Lk(tree,alldata);
323   f1=-tree->tot_loglk;
324   (*param)=x2;
325
326   Lk(tree,alldata);
327   f2=-tree->tot_loglk;
328
329   n_iter = 0;
330   while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) 
331     {
332
333       if (f2 < f1) 
334         {
335           SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
336           (*param)=x2;
337           Lk(tree,alldata);
338           SHFT2(f1,f2,-tree->tot_loglk)
339         } 
340       else 
341         {
342           SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
343           (*param)=x1;
344           Lk(tree,alldata);
345           SHFT2(f2,f1,-tree->tot_loglk)
346         }
347       
348       if(n_iter++ > n_iter_max) break;
349       
350/*        printf("p=%E %f\n",(*param),tree->tot_loglk); */
351     }
352   if (f1 < f2) 
353    {
354       *xmin=x1;
355       return f1;
356     } 
357   else 
358     {
359       *xmin=x2;
360       return f2;
361     }
362}
363
364/*********************************************************/
365
366double Br_Len_Golden(double ax, double bx, double cx, double tol, 
367                     double *xmin, edge *b_fcus, arbre *tree)
368{
369   double f1,f2,x0,x1,x2,x3;
370
371   x0=ax;
372   x3=cx;
373   if (fabs(cx-bx) > fabs(bx-ax)) 
374     {
375       x1=bx;
376       x2=bx+GOLDEN_C*(cx-bx);
377     } 
378   else 
379     {
380       x2=bx;
381       x1=bx-GOLDEN_C*(bx-ax);
382     }
383   
384   b_fcus->l=x1;
385   f1 = -Lk_At_Given_Edge(tree,b_fcus);
386   b_fcus->l=x2;
387   f2 = -Lk_At_Given_Edge(tree,b_fcus);
388   while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) 
389     {
390       if (f2 < f1) 
391         {
392           SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
393           b_fcus->l=x2;
394           SHFT2(f1,f2,-Lk_At_Given_Edge(tree,b_fcus))
395         } 
396       else 
397         {
398           SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
399           b_fcus->l=x1;
400           SHFT2(f2,f1,-Lk_At_Given_Edge(tree,b_fcus))
401         }
402     }
403   if (f1 < f2) 
404     {
405       *xmin=fabs(x1);
406       return -f1;
407     } 
408   else 
409     {
410       *xmin=fabs(x2);
411       return -f2;
412     }
413}
414
415/*********************************************************/
416
417int Br_Len_Brak(double *ax, double *bx, double *cx, 
418                double *fa, double *fb, double *fc, 
419                edge *b_fcus, arbre *tree)
420{
421   double ulim,u,r,q,fu,dum;
422
423   b_fcus->l = *ax;
424   *fa=-Lk_At_Given_Edge(tree,b_fcus);
425   b_fcus->l = *bx;
426   *fb=-Lk_At_Given_Edge(tree,b_fcus);
427   if (*fb > *fa) {
428      SHFT(dum,*ax,*bx,dum)
429      SHFT(dum,*fb,*fa,dum)
430   }
431   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
432   b_fcus->l = fabs(*cx);
433   *fc=-Lk_At_Given_Edge(tree,b_fcus);
434   while (*fb > *fc) 
435     {
436
437       r=(*bx-*ax)*(*fb-*fc);
438       q=(*bx-*cx)*(*fb-*fa);
439       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
440               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
441       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
442       
443       if ((*bx-u)*(u-*cx) > 0.0) 
444         {
445           b_fcus->l = fabs(u);
446           fu=-Lk_At_Given_Edge(tree,b_fcus);
447           if (fu < *fc) 
448             {
449               *ax=(*bx);
450               *bx=u;
451               *fa=(*fb);
452               *fb=fu;
453               (*ax)=fabs(*ax);
454               (*bx)=fabs(*bx);
455               (*cx)=fabs(*cx);
456               return(0);
457             } 
458           else if (fu > *fb) 
459             {
460               *cx=u;
461               *fc=fu; 
462               (*ax)=fabs(*ax);
463               (*bx)=fabs(*bx);
464               (*cx)=fabs(*cx);
465               return(0);
466             }
467           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
468           b_fcus->l = fabs(u);
469           fu=-Lk_At_Given_Edge(tree,b_fcus);
470         } 
471       else if ((*cx-u)*(u-ulim) > 0.0) 
472         {
473           b_fcus->l = fabs(u);
474           fu=-Lk_At_Given_Edge(tree,b_fcus);
475           if (fu < *fc) 
476             {
477               SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
478               b_fcus->l = fabs(u); 
479               SHFT(*fb,*fc,fu,-Lk_At_Given_Edge(tree,b_fcus))
480             }
481         } 
482       else if ((u-ulim)*(ulim-*cx) >= 0.0) 
483         {
484           u=ulim;
485           b_fcus->l = fabs(u);
486           fu=-Lk_At_Given_Edge(tree,b_fcus);
487         } 
488       else 
489         {
490           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
491           b_fcus->l = fabs(u);
492           fu=-Lk_At_Given_Edge(tree,b_fcus);
493         }
494       SHFT(*ax,*bx,*cx,u)
495       SHFT(*fa,*fb,*fc,fu)
496      }
497   (*ax)=fabs(*ax);
498   (*bx)=fabs(*bx);
499   (*cx)=fabs(*cx);
500   return(0);
501}
502
503/*********************************************************/
504
505double Br_Len_Brent(double ax, double bx, double cx, double tol,
506                    double *xmin, edge *b_fcus, arbre *tree, int n_iter_max)
507{
508  int iter;
509  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
510  double e=0.0;
511 
512  d=0.0;
513  a=((ax < cx) ? ax : cx);
514  b=((ax > cx) ? ax : cx);
515  x=w=v=bx;
516  b_fcus->l = fabs(bx);
517  fw=fv=fx=-Lk_At_Given_Edge(tree,b_fcus);
518 
519  for(iter=1;iter<=BRENT_ITMAX;iter++)
520    {
521      xm=0.5*(a+b);
522      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
523      if(fabs(x-xm) <= (tol2-0.5*(b-a)))
524        {
525          *xmin=x;
526          Lk_At_Given_Edge(tree,b_fcus);
527          return -fx;
528        }
529     
530      if(fabs(e) > tol1)
531        {
532          r=(x-w)*(fx-fv);
533          q=(x-v)*(fx-fw);
534          p=(x-v)*q-(x-w)*r;
535          q=2.0*(q-r);
536          if(q > 0.0) p = -p;
537          q=fabs(q);
538          etemp=e;
539          e=d;
540          if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
541            d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
542          else{
543            d=p/q;
544            u=x+d;
545            if (u-a < tol2 || b-u < tol2)
546              d=SIGN(tol1,xm-x);
547          }
548        }
549      else
550        {
551          d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
552        }
553      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
554      if(u<BL_MIN) u = BL_MIN;
555      b_fcus->l=fabs(u);
556      fu=-Lk_At_Given_Edge(tree,b_fcus);
557/*       printf("edge %d l=%f lnL=%f\n",b_fcus->num,b_fcus->l,fu); */
558      if(fu <= fx)
559        {
560          if(iter > n_iter_max) 
561            {
562              printf("\n. WARNING : too many iterations in Brent\n");
563              b_fcus->l = fabs(bx);
564              Lk_At_Given_Edge(tree,b_fcus);
565              return tree->tot_loglk;
566            }
567
568          if(u >= x) a=x; else b=x;
569          SHFT(v,w,x,u)
570            SHFT(fv,fw,fx,fu)
571            }
572      else
573        {
574          if (u < x) a=u; else b=u;
575          if (fu <= fw || w == x)
576            {
577              v=w;
578              w=u;
579              fv=fw;
580            fw=fu;
581            }
582          else if (fu <= fv || v == x || v == w) {
583            v=u;
584            fv=fu;
585          }
586        }
587    }
588  printf("Too many iterations in BRENT");
589  return(-1);
590  /* Not Reached ??  *xmin=x;   */
591  /* Not Reached ??  return fx; */
592}
593
594/*********************************************************/
595
596int Dist_Seq_Brak(double *ax, double *bx, double *cx, 
597                  double *fa, double *fb, double *fc, 
598                  allseq *data, int numseq1, int numseq2, 
599                  model *mod)
600{
601   double ulim,u,r,q,fu,dum;
602   double dist;
603   double lk,dlk,d2lk;
604
605   dist = *ax;
606   *fa=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
607   dist = *bx;
608   *fb=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
609   if (*fb > *fa) {
610      SHFT(dum,*ax,*bx,dum)
611      SHFT(dum,*fb,*fa,dum)
612   }
613   *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
614   dist = fabs(*cx);
615   *fc=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
616   while (*fb > *fc) 
617     {
618       r=(*bx-*ax)*(*fb-*fc);
619       q=(*bx-*cx)*(*fb-*fa);
620       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
621               (2.0*SIGN(MAX(fabs(q-r),MNBRAK_TINY),q-r));
622       ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
623       
624       if ((*bx-u)*(u-*cx) > 0.0) 
625         {
626           dist = fabs(u);
627           fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
628           if (fu < *fc) 
629             {
630               *ax=(*bx);
631               *bx=u;
632               *fa=(*fb);
633               *fb=fu;
634               return(0);
635             } 
636           else if (fu > *fb) 
637             {
638               *cx=u;
639               *fc=fu;
640               return(0);
641             }
642           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
643           dist = fabs(u);
644           fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
645         } 
646       else if ((*cx-u)*(u-ulim) > 0.0) 
647         {
648           dist = fabs(u);
649           fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
650           if (fu < *fc) 
651             {
652               SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
653               dist = fabs(u); 
654               SHFT(*fb,*fc,fu,-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk))
655             }
656         } 
657       else if ((u-ulim)*(ulim-*cx) >= 0.0) 
658         {
659           u=ulim;
660           dist = fabs(u);
661           fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
662         } 
663       else 
664         {
665           u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
666           dist = fabs(u);
667           fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
668         }
669       SHFT(*ax,*bx,*cx,u)
670       SHFT(*fa,*fb,*fc,fu)
671      }
672   return(0);
673}
674
675/*********************************************************/
676
677double Dist_Seq_Brent(double ax, double bx, double cx, double tol, 
678                      double *xmin, allseq *data, 
679                      int numseq1, int numseq2, model *mod)
680{
681  int iter;
682  double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
683  double e=0.0;
684  double dist;
685  double lk,dlk,d2lk;
686 
687  d=0.0;
688  a=((ax < cx) ? ax : cx);
689  b=((ax > cx) ? ax : cx);
690  x=w=v=bx;
691  dist = fabs(bx);
692  fw=fv=fx=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
693  for(iter=1;iter<=BRENT_ITMAX;iter++) 
694    {
695      xm=0.5*(a+b);
696      tol2=2.0*(tol1=tol*fabs(x)+BRENT_ZEPS);
697      if(fabs(x-xm) <= (tol2-0.5*(b-a))) 
698        {
699          *xmin=x;
700          return -fx;
701        }
702     
703      if(fabs(e) > tol1) 
704        {
705          r=(x-w)*(fx-fv);
706          q=(x-v)*(fx-fw);
707          p=(x-v)*q-(x-w)*r;
708          q=2.0*(q-r);
709          if(q > 0.0) p = -p;
710          q=fabs(q);
711          etemp=e;
712          e=d;
713          if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
714            d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
715          else{
716            d=p/q;
717            u=x+d;
718            if (u-a < tol2 || b-u < tol2)
719              d=SIGN(tol1,xm-x);
720          }
721        } else
722          {
723            d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
724          }
725      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
726      dist=fabs(u);
727      fu=-Lk_Given_Two_Seq(data,numseq1,numseq2,dist,mod,&lk,&dlk,&d2lk);
728      if(fu <= fx) {
729        if(u >= x) a=x; else b=x;
730        SHFT(v,w,x,u)
731          SHFT(fv,fw,fx,fu)
732          } 
733      else
734        {
735          if (u < x) a=u; else b=u;
736          if (fu <= fw || w == x) 
737            {
738              v=w;
739              w=u;
740              fv=fw;
741              fw=fu;
742            } 
743          else if (fu <= fv || v == x || v == w) {
744            v=u;
745            fv=fu;
746          }
747        }
748    }
749  printf("\n . BRENT method failed, trying Newton-Raphson");
750  return(+1.0);
751  /* Not Reached ??  *xmin=x;   */
752  /* Not Reached ??  return fx; */
753}
754
755/*********************************************************/
756
757double Optimize_Dist(model *mod, double init, allseq *twoseqs)
758{
759  double d_infa,d_max,d_infb;
760  double lk_infa, lk_max, lk_infb, lk;
761
762  d_infa = 100.*BL_MIN;
763  d_max  = init;
764  d_infb = 3.*init;
765  if(init <= BL_MIN) {d_infa = -BL_START; d_max = .0; d_infb = BL_START;}
766  lk_infa = lk_max = lk_infb = .0;
767
768  Dist_Seq_Brak(&d_infa, &d_max, &d_infb,
769                &lk_infa,&lk_max,&lk_infb,
770                twoseqs,0,1,mod);
771
772  lk = (double)Dist_Seq_Brent(d_infa,d_max,d_infb,
773                              1.e-5,&d_max,twoseqs,0,1,mod);
774  if(lk > .0) return -1.0;
775  else        return d_max;
776
777}
778
779/*********************************************************/
780
781void Round_Optimize(arbre *tree, allseq *data)
782{
783  int n_round,each;
784  double lk_old, lk_new, tol;
785  node *root;
786
787  lk_new = tree->tot_loglk;
788  lk_old = UNLIKELY;
789  n_round = 0;
790  each = 1;
791  tol = 1.e-2;
792  root = tree->noeud[0];
793
794
795  tree->mod->s_opt->opt_bl = 0;
796  tree->both_sides         = 1;
797  Lk(tree,data);
798
799
800  while(n_round < ROUND_MAX)
801    {
802
803        (!((n_round+2)%2))?(root=tree->noeud[0]):(root=tree->noeud[tree->n_otu-1]);
804        Optimize_Br_Len_Serie(root,
805                              root->v[0],
806                              root->b[0],
807                              tree,
808                              data,
809                              5);
810       
811
812        tree->mod->s_opt->opt_bl = 0;
813        tree->both_sides         = 1;
814        Lk(tree,data);
815       
816        if(!each)
817            {
818                each = 1;
819                if(tree->mod->s_opt->print) printf("\n");
820                Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
821                tree->mod->s_opt->opt_bl = 0;
822                tree->both_sides         = 1;
823                Lk(tree,data);
824            }
825       
826        lk_new = tree->tot_loglk;
827       
828        if(tree->mod->s_opt->print)
829            {
830                if(lk_old < UNLIKELY+1)
831                    printf("\n. Log(lk) :               * -> %15.6f ",lk_new);
832                else
833                    printf("\n. Log(lk) : %15.6f -> %15.6f ",lk_old,lk_new);
834                fflush(NULL);
835            }
836       
837        if(lk_new < lk_old - MIN_DIFF_LK) Exit("\n. Optimisation failed ! (Round_Optimize)\n");
838        if(fabs(lk_new - lk_old) < MIN_DIFF_LK)  break;
839        else lk_old  = lk_new;
840        n_round++;
841        each--;
842    }
843 
844  if(tree->mod->s_opt->print) printf("\n");
845  Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
846}
847
848/*********************************************************/
849
850void Optimize_Br_Len_Serie(node *a, node *d, edge *b_fcus, 
851                           arbre *tree,allseq *alldata, int n_passes)
852{
853  int i;
854  double l_infa,l_max,l_infb;
855  double lk_init;
856 
857  lk_init = tree->tot_loglk;
858
859  l_infa = 10.*b_fcus->l;
860  l_max  = b_fcus->l;
861  l_infb = BL_MIN;
862 
863
864  Br_Len_Brent(l_infa,l_max,l_infb,
865               1.e-5,
866               &(b_fcus->l),
867               b_fcus,tree,1000);
868
869  /* Golden method is generally slower than Brent */
870  /*   Br_Len_Golden(l_infa,l_max,l_infb, */
871  /*            1.e-5, */
872  /*            &(b_fcus->l), */
873  /*            b_fcus,tree); */
874/*   printf("Edge %d -> %20f\n",b_fcus->num,tree->tot_loglk); */
875  if(tree->tot_loglk < lk_init - MIN_DIFF_LK)
876    {
877      printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
878      printf("%f -- %f \n",lk_init,tree->tot_loglk);
879      Exit("\n. Err. in Optimize_Br_Len_Serie\n");
880    }
881
882/*   printf("Edge %3d -> %f %f\n", */
883/*       b_fcus->num, */
884/*       tree->tot_loglk, */
885/*       b_fcus->l); fflush(NULL); */
886   
887 
888  if(d->tax) return;
889  else For(i,3) if(d->v[i] != a)
890    {
891      Update_P_Lk(tree,d->b[i],d);
892      Optimize_Br_Len_Serie(d,d->v[i],d->b[i],tree,alldata,n_passes);
893    }
894  For(i,3) if((d->v[i] == a) && !(d->v[i]->tax)) Update_P_Lk(tree,d->b[i],d);
895}
896
897/*********************************************************/
898
899double Br_Len_NR(arbre *tree,allseq *alldata, edge *b_fcus, int n_passes)
900{
901  int n_iter, n_iter_mov,edge_num;
902  double step_edge,lk_new,lk_old;
903  double l_branch_best;
904  int l_r,r_l;
905  double mult;
906  int old_opt;
907
908  old_opt =              tree->mod->s_opt->opt_bl;
909  tree->mod->s_opt->opt_bl = tree->both_sides = 1;
910  edge_num =                     b_fcus->num;
911  l_branch_best =                  b_fcus->l;
912  step_edge =                             .0;
913
914  l_r = b_fcus->l_r;
915  r_l = b_fcus->r_l;
916
917  lk_new = lk_old = UNLIKELY; 
918
919  n_iter = 0;
920  for(;;)
921    {     
922      lk_new = Lk_At_Given_Edge(tree,b_fcus);
923
924      if(fabs(lk_new - lk_old) < MIN_DIFF_LK)
925        break;
926      else
927        {
928          if((lk_new > lk_old+MIN_DIFF_LK) || (!n_iter))
929            {
930              lk_old = lk_new;
931              l_branch_best = b_fcus->l;
932            }
933
934          else if((lk_new < lk_old-MIN_DIFF_LK) && (lk_old != 0.0)) 
935            {
936              tree->mod->s_opt->opt_bl = tree->both_sides = 0;
937              n_iter_mov = 0;
938              while(lk_new < lk_old-MIN_DIFF_LK)
939                {
940                  step_edge *= .5;
941                  b_fcus->l = l_branch_best + step_edge;
942                  if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN;
943                  lk_new = Lk_At_Given_Edge(tree,b_fcus);
944                  n_iter_mov++;
945                  if(n_iter_mov > 10) 
946                    {
947                      b_fcus->l = l_branch_best;
948                      lk_new = Lk_At_Given_Edge(tree,b_fcus);
949                      tree->mod->s_opt->opt_bl = tree->both_sides = old_opt;
950                      return lk_new;
951                    }
952                }
953              tree->mod->s_opt->opt_bl = tree->both_sides = 1;
954              lk_new = Lk_At_Given_Edge(tree,b_fcus);
955              tree->tot_dloglk[edge_num] = .0;
956            }
957        }
958           
959      mult = 0;
960      step_edge = -tree->tot_dloglk[edge_num] / tree->tot_d2loglk[edge_num];
961
962      if(tree->tot_d2loglk[edge_num] > .0)
963        {
964          b_fcus->l /= 1.2;
965          step_edge = 0;
966        }
967     
968      while(fabs(step_edge) > b_fcus->l) step_edge /= 1.2;
969      b_fcus->l += step_edge;
970
971      if(b_fcus->l < BL_MIN) 
972        {
973          b_fcus->l = BL_MIN;
974          break;
975        }
976
977      n_iter++;
978      if(n_iter > n_passes) break;
979
980    }
981  tree->mod->s_opt->opt_bl = tree->both_sides = old_opt;
982  return lk_new;
983}
984
985/*********************************************************/
986
987void Moving_Backward(arbre *tree, allseq *alldata, 
988                     double **step, double lk_base, double *l_branch_init)
989{
990  double lk_new;
991  int i,n_iter_max;
992
993
994  tree->mod->s_opt->opt_bl = tree->both_sides = 0;
995  n_iter_max = 10;
996  lk_new = lk_base-MIN_DIFF_LK-1.;
997 
998
999  For(i,2*tree->n_otu-3) 
1000    {
1001      tree->t_edges[i]->l = l_branch_init[i];
1002      tree->tot_dloglk[i] = .0;
1003      tree->tot_d2loglk[i] = -.1;
1004    }
1005
1006  while(lk_new < lk_base-MIN_DIFF_LK)
1007    {
1008      if(!n_iter_max) 
1009        {
1010          For(i,2*tree->n_otu-3) tree->t_edges[i]->l = l_branch_init[i];
1011          tree->mod->s_opt->opt_bl = tree->both_sides = 1;
1012          return;
1013        }
1014
1015      For(i,2*tree->n_otu-3)
1016        {
1017          (*step)[i] *= .5;
1018          tree->t_edges[i]->l += (*step)[i];
1019          if(tree->t_edges[i]->l < BL_MIN) tree->t_edges[i]->l = BL_MIN;
1020        }
1021
1022      Lk(tree,alldata);
1023      lk_new = tree->tot_loglk;
1024      n_iter_max--;
1025    }
1026  tree->mod->s_opt->opt_bl = tree->both_sides = 1;
1027}
1028
1029
1030/*********************************************************/
1031
1032double Optimize_One_Dist(allseq *data, int numseq1, int numseq2, double init_dist, model *mod)
1033{
1034
1035  int n_iter, n_iter_mov;
1036  double lk_new,lk_old,lk_best;
1037  double step;
1038  double dlk,d2lk;
1039  double dist_new,dist_best;
1040
1041  step = .0;
1042  dlk  = d2lk = .0;
1043  lk_old = lk_new = lk_best = .0;
1044  dist_new = dist_best = init_dist;
1045 
1046  n_iter = 0;
1047  for(;;)
1048    {
1049     
1050      if(dist_new < BL_MIN) dist_new = BL_MIN;
1051
1052      Lk_Given_Two_Seq(data,numseq1,numseq2,dist_new,mod,&lk_new,&dlk,&d2lk);
1053
1054      if(((fabs(lk_new - lk_old) < MIN_DIFF_LK) && (lk_new >= lk_old)))
1055        break;
1056      else
1057        {
1058          if(((fabs(lk_new - lk_old) > MIN_DIFF_LK) 
1059              && (lk_new > lk_old)) 
1060              || (!n_iter))
1061            {
1062              dist_best = dist_new;
1063              lk_best   = lk_new;
1064              lk_old    = lk_new;
1065            }
1066          else
1067            {
1068              n_iter_mov = 20;
1069              do
1070                {
1071                  dist_new = dist_best;
1072                  step /= 2.;
1073                  dist_new += step;
1074                  Lk_Given_Two_Seq(data,numseq1,numseq2,
1075                                   dist_new,mod,&lk_new,&dlk,&d2lk);
1076                  n_iter_mov--;
1077                  if(!n_iter_mov) return dist_best;
1078                }while(lk_new < lk_best);
1079              dlk = .0;
1080            }
1081        }
1082     
1083      step = -dlk/d2lk;
1084
1085      if(d2lk > 0.0) 
1086        {
1087          step = .0;
1088          dist_new /= 1.5;
1089        }
1090     
1091      while(fabs(step) > dist_new) step /= 1.5;
1092      dist_new += step;
1093
1094     
1095      n_iter++;
1096      if(n_iter > 50) break;
1097    }
1098  return dist_new;
1099}
1100
1101/*********************************************************/
1102
1103void Print_Lk_Progress(arbre *tree, double lk_new, double lk_old, int n_iter)
1104{
1105  if(!n_iter)
1106    printf("\n. Log(lk) :               * -> %15.6f ",lk_new);
1107  else
1108    printf("\n. Log(lk) : %15.6f -> %15.6f ",lk_old,lk_new);
1109  fflush(stdout);
1110}
1111
1112/*********************************************************/
1113
1114int Count_Swap(arbre *tree)
1115{
1116  int i;
1117 
1118  tree->n_swap = 0;
1119  For(i,2*tree->n_otu-3)
1120    {
1121      if((!tree->t_edges[i]->left->tax) &&
1122         (!tree->t_edges[i]->rght->tax))
1123        {
1124          if(tree->t_edges[i]->diff_lk > -2.0)
1125            {
1126              tree->n_swap++;
1127            }
1128        }
1129    }
1130  return tree->n_swap;
1131}
1132
1133/*********************************************************/
1134
1135void Optimiz_Ext_Br(arbre *tree)
1136{
1137  int i;
1138  edge *b;
1139  double l_infa,l_max,l_infb,l_init;
1140  double lk, lk_init;
1141 
1142  lk_init = tree->tot_loglk;
1143 
1144  For(i,2*tree->n_otu-3)
1145    {
1146      b = tree->t_edges[i];
1147      if((b->left->tax) || (b->rght->tax))
1148        {
1149
1150          l_init = b->l;
1151 
1152          l_infa = 100.;
1153          l_max  = b->l;
1154          l_infb = -10.;
1155         
1156          lk = Br_Len_Brent(l_infa,l_max,l_infb,
1157                            1.e-5,
1158                            &(b->l),
1159                            b,tree,1000);
1160          b->ql[0] = b->l;
1161          b->best_conf = 1;
1162          b->l = l_init;
1163        }
1164    }
1165  tree->tot_loglk = lk_init; 
1166}
1167
1168/*********************************************************/
1169
1170void Optimiz_All_Free_Param(arbre *tree, int verbose)
1171{
1172  int  init_both_sides, init_derivatives;
1173 
1174
1175  init_both_sides     = tree->both_sides;
1176  init_derivatives    = tree->mod->s_opt->opt_bl;
1177  tree->both_sides    = 0;
1178  tree->mod->s_opt->opt_bl = 0;
1179
1180 
1181  if((tree->mod->whichmodel == 7) ||
1182     ((tree->mod->whichmodel == 8) && 
1183      (tree->mod->s_opt->opt_rr_param) && 
1184      (tree->mod->n_diff_rr_param > 1)))
1185    {
1186        int failed;
1187       
1188        failed = 0;
1189
1190        if(verbose) 
1191            {
1192                (tree->mod->whichmodel == 7)?
1193                    (printf("\n. Optimisation of the GTR parameters...\n")):   
1194                    (printf("\n. Optimisation of the custom model parameters...\n"));
1195            }
1196       
1197        tree->mod->update_eigen = 1;
1198        BFGS(tree,tree->mod->rr_param_values,tree->mod->n_diff_rr_param,1.e-5,1.e-7,
1199             &Return_Abs_Lk,
1200             &Num_Derivative_Several_Param ,
1201             &Lnsrch_RR_Param,&failed);
1202       
1203        if(failed)
1204            {
1205                int i;
1206               
1207                printf("\n. Optimising one-by-one...\n");
1208               
1209                For(i,tree->mod->n_diff_rr_param) 
1210                    if(i != 5)
1211                        Optimize_Single_Param_Generic(tree,&(tree->mod->rr_param_values[i]),tree->mod->rr_param_values[i],1.E-20,1.E+10,1000);
1212            }
1213       
1214        tree->mod->update_eigen = 0;
1215     
1216    }
1217
1218  if(tree->mod->s_opt->opt_kappa)
1219    {
1220      if(verbose) printf("\n. Optimisation of the ts/tv ratio...\n");fflush(stdout);   
1221      Optimize_Single_Param_Generic(tree,&(tree->mod->kappa),tree->mod->kappa,0.1,100,100);
1222      /*              printf("kappa = %f\n",tree->mod->kappa); */
1223     
1224    }
1225 
1226  if(tree->mod->s_opt->opt_lambda) 
1227    {
1228      Optimize_Single_Param_Generic(tree,&(tree->mod->lambda),tree->mod->lambda,0.001,100,50);
1229      /*              printf("lambda = %f\n",tree->mod->lambda); */
1230    }
1231 
1232  if(tree->mod->s_opt->opt_pinvar) 
1233    {
1234      if(verbose) printf("\n. Optimisation of the proportion of invariable sites...\n");fflush(stdout);   
1235      tree->mod->pinvar = 0.5;
1236      Optimize_Single_Param_Generic(tree,&(tree->mod->pinvar),tree->mod->pinvar,0.0001,0.9999,100);           
1237/*       printf("p-invar = %f\n",tree->mod->pinvar); */
1238    }
1239
1240  if(tree->mod->s_opt->opt_alpha) 
1241    { 
1242      if(verbose) printf("\n. Optimisation of the gamma shape parameter...\n");fflush(stdout);
1243      Optimize_Single_Param_Generic(tree,&(tree->mod->alpha),tree->mod->alpha,0.01,100,100);
1244/*       printf("alpha = %f %f\n",tree->mod->alpha,Return_Lk(tree));       */
1245    }
1246
1247  if(tree->mod->s_opt->opt_bfreq)
1248    {
1249        int failed,i;
1250       
1251        failed = 0;
1252        tree->mod->update_eigen = 1;
1253        if(verbose) printf("\n. Optimisation of nucleotide frequencies...\n");
1254        BFGS(tree,tree->mod->pi,4,1.e-5,1.e-7,&Return_Abs_Lk,&Num_Derivative_Several_Param,&Lnsrch_Nucleotide_Frequencies,&failed);
1255
1256        if(failed)
1257            {
1258                For(i,5) 
1259                    Optimize_Single_Param_Generic(tree,&(tree->mod->pi[i]),tree->mod->pi[i],1.E-10,0.999999,1000);               
1260            }
1261        tree->mod->update_eigen = 0;
1262    }
1263
1264
1265  tree->both_sides    = init_both_sides;
1266  tree->mod->s_opt->opt_bl = init_derivatives;
1267
1268}
1269
1270
1271
1272#define ITMAX 200
1273#define EPS   3.0e-8
1274#define TOLX (4*EPS)
1275#define STPMX 100.0
1276static double sqrarg;
1277#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
1278
1279void BFGS(arbre *tree, double *p, int n, double gtol, double step_size,
1280          double(*func)(), void (*dfunc)(), void (*lnsrch)(),int *failed)
1281{
1282
1283  int check,i,its,j;
1284  double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test,fret;
1285  double *dg,*g,*hdg,**hessin,*pnew,*xi;
1286 
1287  hessin = (double **)mCalloc(n,sizeof(double *));
1288  For(i,n) hessin[i] = (double *)mCalloc(n,sizeof(double));
1289  dg   = (double *)mCalloc(n,sizeof(double ));
1290  g    = (double *)mCalloc(n,sizeof(double ));
1291  pnew = (double *)mCalloc(n,sizeof(double ));
1292  hdg  = (double *)mCalloc(n,sizeof(double ));
1293  xi   = (double *)mCalloc(n,sizeof(double ));
1294 
1295
1296  fp=(*func)(tree);
1297  (*dfunc)(tree,p,n,step_size,func,g);
1298
1299
1300  for (i=0;i<n;i++) 
1301    {
1302      for (j=0;j<n;j++) hessin[i][j]=0.0;
1303      hessin[i][i]=1.0;
1304      xi[i] = -g[i];
1305      sum += p[i]*p[i];
1306    }
1307
1308  stpmax=STPMX*MAX(sqrt(sum),(double)n);
1309
1310  for(its=1;its<=ITMAX;its++) 
1311    {
1312      lnsrch(tree,n,p,fp,g,xi,pnew,&fret,stpmax,&check);
1313
1314/*       printf("BFGS -> %f\n",tree->tot_loglk); */
1315
1316      fp = fret;
1317     
1318      for (i=0;i<n;i++) 
1319        {
1320          xi[i]=pnew[i]-p[i];
1321          p[i]=pnew[i];
1322        }
1323
1324      test=0.0;
1325      for (i=0;i<n;i++) 
1326        {
1327          temp=fabs(xi[i])/MAX(fabs(p[i]),1.0);
1328          if (temp > test) test=temp;
1329        }
1330      if (test < TOLX) 
1331        {
1332          (*func)(tree);
1333          For(i,n) Free(hessin[i]);
1334          free(hessin);
1335          free(xi);
1336          free(pnew);
1337          free(hdg);
1338          free(g);
1339          free(dg);   
1340
1341          if(its == 1) 
1342              {
1343                  printf("\n. WARNING : BFGS failed ! \n");
1344                  *failed = 1;
1345              }
1346          return;
1347        }
1348
1349      for (i=0;i<n;i++) dg[i]=g[i];
1350
1351      (*dfunc)(tree,p,n,step_size,func,g);
1352
1353      test=0.0;
1354      den=MAX(fret,1.0);
1355      for (i=0;i<n;i++) 
1356        {
1357          temp=fabs(g[i])*MAX(fabs(p[i]),1.0)/den;
1358          if (temp > test) test=temp;
1359        }
1360      if (test < gtol) 
1361        {
1362          (*func)(tree);
1363          For(i,n) Free(hessin[i]);
1364          free(hessin);
1365          free(xi);
1366          free(pnew);
1367          free(hdg);
1368          free(g);
1369          free(dg);   
1370          return;
1371        }
1372
1373    for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
1374
1375    for (i=0;i<n;i++) 
1376      {
1377        hdg[i]=0.0;
1378        for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
1379      }
1380
1381    fac=fae=sumdg=sumxi=0.0;
1382    for (i=0;i<n;i++) 
1383      {
1384        fac += dg[i]*xi[i];
1385        fae += dg[i]*hdg[i];
1386        sumdg += SQR(dg[i]);
1387        sumxi += SQR(xi[i]);
1388      }
1389   
1390    if(fac*fac > EPS*sumdg*sumxi) 
1391      {
1392        fac=1.0/fac;
1393        fad=1.0/fae;
1394        for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
1395        for (i=0;i<n;i++) 
1396          {
1397            for (j=0;j<n;j++) 
1398              {
1399                hessin[i][j] += fac*xi[i]*xi[j]
1400                  -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
1401              }
1402          }
1403      }
1404    for (i=0;i<n;i++) 
1405      {
1406        xi[i]=0.0;
1407        for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
1408      }
1409    }
1410  Exit("\n. Too many iterations in BFGS...\n");
1411  For(i,n) Free(hessin[i]);
1412  free(hessin);
1413  free(xi);
1414  free(pnew);
1415  free(hdg);
1416  free(g);
1417  free(dg);   
1418}
1419
1420#undef ITMAX
1421#undef EPS
1422#undef TOLX
1423#undef STPMX
1424
1425/*********************************************************/
1426
1427
1428#define ALF 1.0e-4
1429#define TOLX 1.0e-7
1430
1431void Lnsrch_RR_Param(arbre *tree, int n, double *xold, double fold, 
1432                     double *g, double *p, double *x,
1433                     double *f, double stpmax, int *check)
1434{
1435  int i;
1436  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1437  double *local_xold;
1438
1439  alam = alam2 = f2 = fold2 = tmplam = .0;
1440
1441  local_xold = (double *)mCalloc(n,sizeof(double));
1442  For(i,n) local_xold[i] = xold[i];
1443
1444
1445  *check=0;
1446  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1447  sum=sqrt(sum);
1448  if(sum > stpmax)
1449    for(i=0;i<n;i++) p[i] *= stpmax/sum;
1450  for(slope=0.0,i=0;i<n;i++)
1451    slope += g[i]*p[i];
1452  test=0.0;
1453  for(i=0;i<n;i++) 
1454    {
1455      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
1456      if (temp > test) test=temp;
1457    }
1458  alamin=TOLX/test;
1459  alam=1.0;
1460  for (;;) 
1461    {
1462      for(i=0;i<n;i++) 
1463        {
1464          x[i]=local_xold[i]+alam*p[i];
1465        }
1466
1467      /**/
1468      for(i=0;i<n;i++)
1469        {
1470          tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
1471          if(tree->mod->rr_param_values[i] < 0.0) break;
1472        }
1473      /**/
1474
1475      if(i==n) 
1476        {
1477          *f=Return_Abs_Lk(tree);
1478/*        printf("loglk = %f\n",*f); */
1479        }
1480      else     *f=1.+fold+ALF*alam*slope;
1481      if (alam < alamin)
1482        {
1483          for (i=0;i<n;i++) 
1484            {
1485              x[i]=local_xold[i];
1486              if(x[i] < .0) break;
1487            }
1488          /**/     
1489          for(i=0;i<n;i++)
1490            {
1491              tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
1492              if(tree->mod->rr_param_values[i] < 0.0) 
1493                tree->mod->rr_param_values[i] = 0.0;
1494            }
1495          /**/
1496
1497          *check=1;
1498          For(i,n) xold[i] = local_xold[i];
1499          Free(local_xold);
1500          return;
1501        } 
1502      else if (*f <= fold+ALF*alam*slope) 
1503        {
1504          For(i,n) xold[i] = local_xold[i];
1505          Free(local_xold); 
1506          return;
1507        }
1508      else 
1509        {
1510          if (alam == 1.0)
1511            tmplam = -slope/(2.0*(*f-fold-slope));
1512          else 
1513            {
1514              rhs1 = *f-fold-alam*slope;
1515              rhs2=f2-fold2-alam2*slope;
1516              a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
1517              b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
1518              if (a == 0.0) tmplam = -slope/(2.0*b);
1519              else 
1520                {
1521                  disc=b*b-3.0*a*slope;
1522                  if (disc<0.0) tmplam = 0.5*alam;
1523                  else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
1524                  else tmplam = -slope/(b+sqrt(disc));
1525                }
1526              if (tmplam>0.5*alam) tmplam=0.5*alam;
1527            }
1528        }
1529      alam2=alam;
1530      f2 = *f;
1531      fold2=fold;
1532      alam=MAX(tmplam,0.1*alam);
1533    }
1534  Free(local_xold);
1535}
1536
1537#undef ALF
1538#undef TOLX
1539#undef NRANSI
1540
1541/*********************************************************/
1542#define ALF 1.0e-4
1543#define TOLX 1.0e-7
1544
1545void Lnsrch_Nucleotide_Frequencies(arbre *tree, int n, double *xold, double fold, double *g, double *p, double *x,
1546                                   double *f, double stpmax, int *check)
1547{
1548  int i;
1549  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1550  double *local_xold;
1551
1552  alam = alam2 = f2 = fold2 = tmplam = .0;
1553
1554  local_xold = (double *)mCalloc(n,sizeof(double));
1555  For(i,n) local_xold[i] = xold[i];
1556
1557
1558  *check=0;
1559  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1560  sum=sqrt(sum);
1561  if(sum > stpmax)
1562    for(i=0;i<n;i++) p[i] *= stpmax/sum;
1563  for(slope=0.0,i=0;i<n;i++)
1564    slope += g[i]*p[i];
1565  test=0.0;
1566  for(i=0;i<n;i++) 
1567    {
1568      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
1569      if (temp > test) test=temp;
1570    }
1571  alamin=TOLX/test;
1572  alam=1.0;
1573  for (;;) 
1574    {
1575      for(i=0;i<n;i++) x[i]=fabs(local_xold[i]+alam*p[i]);
1576      /**/     
1577      for(i=0;i<n;i++) 
1578        {
1579          tree->mod->pi[i]=fabs(local_xold[i]+alam*p[i]);
1580/*        if( */
1581/*           (tree->mod->pi[i] < 0.001) || */
1582/*           (tree->mod->pi[i] > 0.999) */
1583/*           ) */
1584/*          break; */
1585        }
1586      /**/
1587      if(i==n) 
1588        {
1589          *f=Return_Abs_Lk(tree);
1590        }
1591      else     *f=1.+fold+ALF*alam*slope;
1592      if (alam < alamin)
1593        {
1594          for (i=0;i<n;i++) x[i]=local_xold[i];
1595          for (i=0;i<n;i++) tree->mod->pi[i]=local_xold[i];
1596          *check=1;
1597          For(i,n) xold[i] = local_xold[i];
1598          Free(local_xold);
1599          return;
1600        } 
1601      else if (*f <= fold+ALF*alam*slope) 
1602        {
1603          For(i,n) xold[i] = local_xold[i];
1604          Free(local_xold); 
1605          return;
1606        }
1607      else 
1608        {
1609          if (alam == 1.0)
1610            tmplam = -slope/(2.0*(*f-fold-slope));
1611          else 
1612            {
1613              rhs1 = *f-fold-alam*slope;
1614              rhs2=f2-fold2-alam2*slope;
1615              a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
1616              b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
1617              if (a == 0.0) tmplam = -slope/(2.0*b);
1618              else 
1619                {
1620                  disc=b*b-3.0*a*slope;
1621                  if (disc<0.0) 
1622                    {
1623                      disc=b*b-3.0*a*slope;
1624                      if (disc<0.0) tmplam = 0.5*alam;
1625                      else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
1626                      else tmplam = -slope/(b+sqrt(disc));
1627                    }
1628                  else tmplam=(-b+sqrt(disc))/(3.0*a);
1629                }
1630              if (tmplam>0.5*alam) tmplam=0.5*alam;
1631            }
1632        }
1633      alam2=alam;
1634      f2 = *f;
1635      fold2=fold;
1636      alam=MAX(tmplam,0.1*alam);
1637    }
1638  Free(local_xold);
1639}
1640
1641/*********************************************************/
1642
1643/* void Optimize_Global_Rate(arbre *tree) */
1644/* { */
1645/*     printf("\n. Global rate (%f->)",tree->tot_loglk); */
1646/*     Optimize_Single_Param_Generic(tree,&(tree->tbl),tree->tbl,BL_MIN,1.E+4,100); */
1647/*     printf("%f)\n",tree->tot_loglk); */
1648/* } */
1649
1650
1651#undef ALF
1652#undef TOLX
1653#undef NRANSI
1654
1655
Note: See TracBrowser for help on using the repository browser.