source: trunk/GDE/PHYML/optimiz.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.0 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;
785  node *root;
786
787  lk_new = tree->tot_loglk;
788  lk_old = UNLIKELY;
789  n_round = 0;
790  each = 1;
791  root = tree->noeud[0];
792
793
794  tree->mod->s_opt->opt_bl = 0;
795  tree->both_sides         = 1;
796  Lk(tree,data);
797
798
799  while(n_round < ROUND_MAX)
800    {
801
802        (!((n_round+2)%2))?(root=tree->noeud[0]):(root=tree->noeud[tree->n_otu-1]);
803        Optimize_Br_Len_Serie(root,
804                              root->v[0],
805                              root->b[0],
806                              tree,
807                              data,
808                              5);
809       
810
811        tree->mod->s_opt->opt_bl = 0;
812        tree->both_sides         = 1;
813        Lk(tree,data);
814       
815        if(!each)
816            {
817                each = 1;
818                if(tree->mod->s_opt->print) printf("\n");
819                Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
820                tree->mod->s_opt->opt_bl = 0;
821                tree->both_sides         = 1;
822                Lk(tree,data);
823            }
824       
825        lk_new = tree->tot_loglk;
826       
827        if(tree->mod->s_opt->print)
828            {
829                if(lk_old < UNLIKELY+1)
830                    printf("\n. Log(lk) :               * -> %15.6f ",lk_new);
831                else
832                    printf("\n. Log(lk) : %15.6f -> %15.6f ",lk_old,lk_new);
833                fflush(NULL);
834            }
835       
836        if(lk_new < lk_old - MIN_DIFF_LK) Exit("\n. Optimisation failed ! (Round_Optimize)\n");
837        if(fabs(lk_new - lk_old) < MIN_DIFF_LK)  break;
838        else lk_old  = lk_new;
839        n_round++;
840        each--;
841    }
842 
843  if(tree->mod->s_opt->print) printf("\n");
844  Optimiz_All_Free_Param(tree,tree->mod->s_opt->print);
845}
846
847/*********************************************************/
848
849void Optimize_Br_Len_Serie(node *a, node *d, edge *b_fcus, 
850                           arbre *tree,allseq *alldata, int n_passes)
851{
852  int i;
853  double l_infa,l_max,l_infb;
854  double lk_init;
855 
856  lk_init = tree->tot_loglk;
857
858  l_infa = 10.*b_fcus->l;
859  l_max  = b_fcus->l;
860  l_infb = BL_MIN;
861 
862
863  Br_Len_Brent(l_infa,l_max,l_infb,
864               1.e-5,
865               &(b_fcus->l),
866               b_fcus,tree,1000);
867
868  /* Golden method is generally slower than Brent */
869  /*   Br_Len_Golden(l_infa,l_max,l_infb, */
870  /*            1.e-5, */
871  /*            &(b_fcus->l), */
872  /*            b_fcus,tree); */
873/*   printf("Edge %d -> %20f\n",b_fcus->num,tree->tot_loglk); */
874  if(tree->tot_loglk < lk_init - MIN_DIFF_LK)
875    {
876      printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
877      printf("%f -- %f \n",lk_init,tree->tot_loglk);
878      Exit("\n. Err. in Optimize_Br_Len_Serie\n");
879    }
880
881/*   printf("Edge %3d -> %f %f\n", */
882/*       b_fcus->num, */
883/*       tree->tot_loglk, */
884/*       b_fcus->l); fflush(NULL); */
885   
886 
887  if(d->tax) return;
888  else For(i,3) if(d->v[i] != a)
889    {
890      Update_P_Lk(tree,d->b[i],d);
891      Optimize_Br_Len_Serie(d,d->v[i],d->b[i],tree,alldata,n_passes);
892    }
893  For(i,3) if((d->v[i] == a) && !(d->v[i]->tax)) Update_P_Lk(tree,d->b[i],d);
894}
895
896/*********************************************************/
897
898void Moving_Backward(arbre *tree, allseq *alldata, 
899                     double **step, double lk_base, double *l_branch_init)
900{
901  double lk_new;
902  int i,n_iter_max;
903
904
905  tree->mod->s_opt->opt_bl = tree->both_sides = 0;
906  n_iter_max = 10;
907  lk_new = lk_base-MIN_DIFF_LK-1.;
908 
909
910  For(i,2*tree->n_otu-3) 
911    {
912      tree->t_edges[i]->l = l_branch_init[i];
913      tree->tot_dloglk[i] = .0;
914      tree->tot_d2loglk[i] = -.1;
915    }
916
917  while(lk_new < lk_base-MIN_DIFF_LK)
918    {
919      if(!n_iter_max) 
920        {
921          For(i,2*tree->n_otu-3) tree->t_edges[i]->l = l_branch_init[i];
922          tree->mod->s_opt->opt_bl = tree->both_sides = 1;
923          return;
924        }
925
926      For(i,2*tree->n_otu-3)
927        {
928          (*step)[i] *= .5;
929          tree->t_edges[i]->l += (*step)[i];
930          if(tree->t_edges[i]->l < BL_MIN) tree->t_edges[i]->l = BL_MIN;
931        }
932
933      Lk(tree,alldata);
934      lk_new = tree->tot_loglk;
935      n_iter_max--;
936    }
937  tree->mod->s_opt->opt_bl = tree->both_sides = 1;
938}
939
940
941/*********************************************************/
942
943double Optimize_One_Dist(allseq *data, int numseq1, int numseq2, double init_dist, model *mod)
944{
945
946  int n_iter, n_iter_mov;
947  double lk_new,lk_old,lk_best;
948  double step;
949  double dlk,d2lk;
950  double dist_new,dist_best;
951
952  step = .0;
953  dlk  = d2lk = .0;
954  lk_old = lk_new = lk_best = .0;
955  dist_new = dist_best = init_dist;
956 
957  n_iter = 0;
958  for(;;)
959    {
960     
961      if(dist_new < BL_MIN) dist_new = BL_MIN;
962
963      Lk_Given_Two_Seq(data,numseq1,numseq2,dist_new,mod,&lk_new,&dlk,&d2lk);
964
965      if(((fabs(lk_new - lk_old) < MIN_DIFF_LK) && (lk_new >= lk_old)))
966        break;
967      else
968        {
969          if(((fabs(lk_new - lk_old) > MIN_DIFF_LK) 
970              && (lk_new > lk_old)) 
971              || (!n_iter))
972            {
973              dist_best = dist_new;
974              lk_best   = lk_new;
975              lk_old    = lk_new;
976            }
977          else
978            {
979              n_iter_mov = 20;
980              do
981                {
982                  dist_new = dist_best;
983                  step /= 2.;
984                  dist_new += step;
985                  Lk_Given_Two_Seq(data,numseq1,numseq2,
986                                   dist_new,mod,&lk_new,&dlk,&d2lk);
987                  n_iter_mov--;
988                  if(!n_iter_mov) return dist_best;
989                }while(lk_new < lk_best);
990              dlk = .0;
991            }
992        }
993     
994      step = -dlk/d2lk;
995
996      if(d2lk > 0.0) 
997        {
998          step = .0;
999          dist_new /= 1.5;
1000        }
1001     
1002      while(fabs(step) > dist_new) step /= 1.5;
1003      dist_new += step;
1004
1005     
1006      n_iter++;
1007      if(n_iter > 50) break;
1008    }
1009  return dist_new;
1010}
1011
1012/*********************************************************/
1013
1014int Count_Swap(arbre *tree)
1015{
1016  int i;
1017 
1018  tree->n_swap = 0;
1019  For(i,2*tree->n_otu-3)
1020    {
1021      if((!tree->t_edges[i]->left->tax) &&
1022         (!tree->t_edges[i]->rght->tax))
1023        {
1024          if(tree->t_edges[i]->diff_lk > -2.0)
1025            {
1026              tree->n_swap++;
1027            }
1028        }
1029    }
1030  return tree->n_swap;
1031}
1032
1033/*********************************************************/
1034
1035void Optimiz_Ext_Br(arbre *tree)
1036{
1037  int i;
1038  edge *b;
1039  double l_infa,l_max,l_infb,l_init;
1040  double lk_init;
1041 
1042  lk_init = tree->tot_loglk;
1043 
1044  For(i,2*tree->n_otu-3)
1045    {
1046      b = tree->t_edges[i];
1047      if((b->left->tax) || (b->rght->tax))
1048        {
1049
1050          l_init = b->l;
1051 
1052          l_infa = 100.;
1053          l_max  = b->l;
1054          l_infb = -10.;
1055         
1056          Br_Len_Brent(l_infa,l_max,l_infb, // Note: kept call here (might have wanted side-effects)
1057                       1.e-5,
1058                       &(b->l),
1059                       b,tree,1000);
1060
1061          b->ql[0] = b->l;
1062          b->best_conf = 1;
1063          b->l = l_init;
1064        }
1065    }
1066  tree->tot_loglk = lk_init; 
1067}
1068
1069/*********************************************************/
1070
1071void Optimiz_All_Free_Param(arbre *tree, int verbose)
1072{
1073  int  init_both_sides, init_derivatives;
1074 
1075
1076  init_both_sides     = tree->both_sides;
1077  init_derivatives    = tree->mod->s_opt->opt_bl;
1078  tree->both_sides    = 0;
1079  tree->mod->s_opt->opt_bl = 0;
1080
1081 
1082  if((tree->mod->whichmodel == 7) ||
1083     ((tree->mod->whichmodel == 8) && 
1084      (tree->mod->s_opt->opt_rr_param) && 
1085      (tree->mod->n_diff_rr_param > 1)))
1086    {
1087        int failed;
1088       
1089        failed = 0;
1090
1091        if(verbose) 
1092            {
1093                (tree->mod->whichmodel == 7)?
1094                    (printf("\n. Optimisation of the GTR parameters...\n")):   
1095                    (printf("\n. Optimisation of the custom model parameters...\n"));
1096            }
1097       
1098        tree->mod->update_eigen = 1;
1099        BFGS(tree,tree->mod->rr_param_values,tree->mod->n_diff_rr_param,1.e-5,1.e-7,
1100             &Return_Abs_Lk,
1101             &Num_Derivative_Several_Param ,
1102             &Lnsrch_RR_Param,&failed);
1103       
1104        if(failed)
1105            {
1106                int i;
1107               
1108                printf("\n. Optimising one-by-one...\n");
1109               
1110                For(i,tree->mod->n_diff_rr_param) 
1111                    if(i != 5)
1112                        Optimize_Single_Param_Generic(tree,&(tree->mod->rr_param_values[i]),tree->mod->rr_param_values[i],1.E-20,1.E+10,1000);
1113            }
1114       
1115        tree->mod->update_eigen = 0;
1116     
1117    }
1118
1119  if(tree->mod->s_opt->opt_kappa)
1120    {
1121      if(verbose) printf("\n. Optimisation of the ts/tv ratio...\n");fflush(stdout);   
1122      Optimize_Single_Param_Generic(tree,&(tree->mod->kappa),tree->mod->kappa,0.1,100,100);
1123      /*              printf("kappa = %f\n",tree->mod->kappa); */
1124     
1125    }
1126 
1127  if(tree->mod->s_opt->opt_lambda) 
1128    {
1129      Optimize_Single_Param_Generic(tree,&(tree->mod->lambda),tree->mod->lambda,0.001,100,50);
1130      /*              printf("lambda = %f\n",tree->mod->lambda); */
1131    }
1132 
1133  if(tree->mod->s_opt->opt_pinvar) 
1134    {
1135      if(verbose) printf("\n. Optimisation of the proportion of invariable sites...\n");fflush(stdout);   
1136      tree->mod->pinvar = 0.5;
1137      Optimize_Single_Param_Generic(tree,&(tree->mod->pinvar),tree->mod->pinvar,0.0001,0.9999,100);           
1138/*       printf("p-invar = %f\n",tree->mod->pinvar); */
1139    }
1140
1141  if(tree->mod->s_opt->opt_alpha) 
1142    { 
1143      if(verbose) printf("\n. Optimisation of the gamma shape parameter...\n");fflush(stdout);
1144      Optimize_Single_Param_Generic(tree,&(tree->mod->alpha),tree->mod->alpha,0.01,100,100);
1145/*       printf("alpha = %f %f\n",tree->mod->alpha,Return_Lk(tree));       */
1146    }
1147
1148  if(tree->mod->s_opt->opt_bfreq)
1149    {
1150        int failed,i;
1151       
1152        failed = 0;
1153        tree->mod->update_eigen = 1;
1154        if(verbose) printf("\n. Optimisation of nucleotide frequencies...\n");
1155        BFGS(tree,tree->mod->pi,4,1.e-5,1.e-7,&Return_Abs_Lk,&Num_Derivative_Several_Param,&Lnsrch_Nucleotide_Frequencies,&failed);
1156
1157        if(failed)
1158            {
1159                For(i,5) 
1160                    Optimize_Single_Param_Generic(tree,&(tree->mod->pi[i]),tree->mod->pi[i],1.E-10,0.999999,1000);               
1161            }
1162        tree->mod->update_eigen = 0;
1163    }
1164
1165
1166  tree->both_sides    = init_both_sides;
1167  tree->mod->s_opt->opt_bl = init_derivatives;
1168
1169}
1170
1171
1172
1173#define ITMAX 200
1174#define EPS   3.0e-8
1175#define TOLX (4*EPS)
1176#define STPMX 100.0
1177static double sqrarg;
1178#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
1179
1180void BFGS(arbre *tree, double *p, int n, double gtol, double step_size,
1181          double(*func)(), void (*dfunc)(), void (*lnsrch)(),int *failed)
1182{
1183
1184  int check,i,its,j;
1185  double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test,fret;
1186  double *dg,*g,*hdg,**hessin,*pnew,*xi;
1187 
1188  hessin = (double **)mCalloc(n,sizeof(double *));
1189  For(i,n) hessin[i] = (double *)mCalloc(n,sizeof(double));
1190  dg   = (double *)mCalloc(n,sizeof(double ));
1191  g    = (double *)mCalloc(n,sizeof(double ));
1192  pnew = (double *)mCalloc(n,sizeof(double ));
1193  hdg  = (double *)mCalloc(n,sizeof(double ));
1194  xi   = (double *)mCalloc(n,sizeof(double ));
1195 
1196
1197  fp=(*func)(tree);
1198  (*dfunc)(tree,p,n,step_size,func,g);
1199
1200
1201  for (i=0;i<n;i++) 
1202    {
1203      for (j=0;j<n;j++) hessin[i][j]=0.0;
1204      hessin[i][i]=1.0;
1205      xi[i] = -g[i];
1206      sum += p[i]*p[i];
1207    }
1208
1209  stpmax=STPMX*MAX(sqrt(sum),(double)n);
1210
1211  for(its=1;its<=ITMAX;its++) 
1212    {
1213      lnsrch(tree,n,p,fp,g,xi,pnew,&fret,stpmax,&check);
1214
1215/*       printf("BFGS -> %f\n",tree->tot_loglk); */
1216
1217      fp = fret;
1218     
1219      for (i=0;i<n;i++) 
1220        {
1221          xi[i]=pnew[i]-p[i];
1222          p[i]=pnew[i];
1223        }
1224
1225      test=0.0;
1226      for (i=0;i<n;i++) 
1227        {
1228          temp=fabs(xi[i])/MAX(fabs(p[i]),1.0);
1229          if (temp > test) test=temp;
1230        }
1231      if (test < TOLX) 
1232        {
1233          (*func)(tree);
1234          For(i,n) Free(hessin[i]);
1235          free(hessin);
1236          free(xi);
1237          free(pnew);
1238          free(hdg);
1239          free(g);
1240          free(dg);   
1241
1242          if(its == 1) 
1243              {
1244                  printf("\n. WARNING : BFGS failed ! \n");
1245                  *failed = 1;
1246              }
1247          return;
1248        }
1249
1250      for (i=0;i<n;i++) dg[i]=g[i];
1251
1252      (*dfunc)(tree,p,n,step_size,func,g);
1253
1254      test=0.0;
1255      den=MAX(fret,1.0);
1256      for (i=0;i<n;i++) 
1257        {
1258          temp=fabs(g[i])*MAX(fabs(p[i]),1.0)/den;
1259          if (temp > test) test=temp;
1260        }
1261      if (test < gtol) 
1262        {
1263          (*func)(tree);
1264          For(i,n) Free(hessin[i]);
1265          free(hessin);
1266          free(xi);
1267          free(pnew);
1268          free(hdg);
1269          free(g);
1270          free(dg);   
1271          return;
1272        }
1273
1274    for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
1275
1276    for (i=0;i<n;i++) 
1277      {
1278        hdg[i]=0.0;
1279        for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
1280      }
1281
1282    fac=fae=sumdg=sumxi=0.0;
1283    for (i=0;i<n;i++) 
1284      {
1285        fac += dg[i]*xi[i];
1286        fae += dg[i]*hdg[i];
1287        sumdg += SQR(dg[i]);
1288        sumxi += SQR(xi[i]);
1289      }
1290   
1291    if(fac*fac > EPS*sumdg*sumxi) 
1292      {
1293        fac=1.0/fac;
1294        fad=1.0/fae;
1295        for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
1296        for (i=0;i<n;i++) 
1297          {
1298            for (j=0;j<n;j++) 
1299              {
1300                hessin[i][j] += fac*xi[i]*xi[j]
1301                  -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
1302              }
1303          }
1304      }
1305    for (i=0;i<n;i++) 
1306      {
1307        xi[i]=0.0;
1308        for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
1309      }
1310    }
1311  Exit("\n. Too many iterations in BFGS...\n");
1312  For(i,n) Free(hessin[i]);
1313  free(hessin);
1314  free(xi);
1315  free(pnew);
1316  free(hdg);
1317  free(g);
1318  free(dg);   
1319}
1320
1321#undef ITMAX
1322#undef EPS
1323#undef TOLX
1324#undef STPMX
1325
1326/*********************************************************/
1327
1328
1329#define ALF 1.0e-4
1330#define TOLX 1.0e-7
1331
1332void Lnsrch_RR_Param(arbre *tree, int n, double *xold, double fold, 
1333                     double *g, double *p, double *x,
1334                     double *f, double stpmax, int *check)
1335{
1336  int i;
1337  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1338  double *local_xold;
1339
1340  alam = alam2 = f2 = fold2 = tmplam = .0;
1341
1342  local_xold = (double *)mCalloc(n,sizeof(double));
1343  For(i,n) local_xold[i] = xold[i];
1344
1345
1346  *check=0;
1347  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1348  sum=sqrt(sum);
1349  if(sum > stpmax)
1350    for(i=0;i<n;i++) p[i] *= stpmax/sum;
1351  for(slope=0.0,i=0;i<n;i++)
1352    slope += g[i]*p[i];
1353  test=0.0;
1354  for(i=0;i<n;i++) 
1355    {
1356      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
1357      if (temp > test) test=temp;
1358    }
1359  alamin=TOLX/test;
1360  alam=1.0;
1361  for (;;) 
1362    {
1363      for(i=0;i<n;i++) 
1364        {
1365          x[i]=local_xold[i]+alam*p[i];
1366        }
1367
1368      /**/
1369      for(i=0;i<n;i++)
1370        {
1371          tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
1372          if(tree->mod->rr_param_values[i] < 0.0) break;
1373        }
1374      /**/
1375
1376      if(i==n) 
1377        {
1378          *f=Return_Abs_Lk(tree);
1379/*        printf("loglk = %f\n",*f); */
1380        }
1381      else     *f=1.+fold+ALF*alam*slope;
1382      if (alam < alamin)
1383        {
1384          for (i=0;i<n;i++) 
1385            {
1386              x[i]=local_xold[i];
1387              if(x[i] < .0) break;
1388            }
1389          /**/     
1390          for(i=0;i<n;i++)
1391            {
1392              tree->mod->rr_param_values[i]=local_xold[i]+alam*p[i];
1393              if(tree->mod->rr_param_values[i] < 0.0) 
1394                tree->mod->rr_param_values[i] = 0.0;
1395            }
1396          /**/
1397
1398          *check=1;
1399          For(i,n) xold[i] = local_xold[i];
1400          Free(local_xold);
1401          return;
1402        } 
1403      else if (*f <= fold+ALF*alam*slope) 
1404        {
1405          For(i,n) xold[i] = local_xold[i];
1406          Free(local_xold); 
1407          return;
1408        }
1409      else 
1410        {
1411          if (alam == 1.0)
1412            tmplam = -slope/(2.0*(*f-fold-slope));
1413          else 
1414            {
1415              rhs1 = *f-fold-alam*slope;
1416              rhs2=f2-fold2-alam2*slope;
1417              a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
1418              b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
1419              if (a == 0.0) tmplam = -slope/(2.0*b);
1420              else 
1421                {
1422                  disc=b*b-3.0*a*slope;
1423                  if (disc<0.0) tmplam = 0.5*alam;
1424                  else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
1425                  else tmplam = -slope/(b+sqrt(disc));
1426                }
1427              if (tmplam>0.5*alam) tmplam=0.5*alam;
1428            }
1429        }
1430      alam2=alam;
1431      f2 = *f;
1432      fold2=fold;
1433      alam=MAX(tmplam,0.1*alam);
1434    }
1435  Free(local_xold);
1436}
1437
1438#undef ALF
1439#undef TOLX
1440#undef NRANSI
1441
1442/*********************************************************/
1443#define ALF 1.0e-4
1444#define TOLX 1.0e-7
1445
1446void Lnsrch_Nucleotide_Frequencies(arbre *tree, int n, double *xold, double fold, double *g, double *p, double *x,
1447                                   double *f, double stpmax, int *check)
1448{
1449  int i;
1450  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1451  double *local_xold;
1452
1453  alam = alam2 = f2 = fold2 = tmplam = .0;
1454
1455  local_xold = (double *)mCalloc(n,sizeof(double));
1456  For(i,n) local_xold[i] = xold[i];
1457
1458
1459  *check=0;
1460  for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1461  sum=sqrt(sum);
1462  if(sum > stpmax)
1463    for(i=0;i<n;i++) p[i] *= stpmax/sum;
1464  for(slope=0.0,i=0;i<n;i++)
1465    slope += g[i]*p[i];
1466  test=0.0;
1467  for(i=0;i<n;i++) 
1468    {
1469      temp=fabs(p[i])/MAX(fabs(local_xold[i]),1.0);
1470      if (temp > test) test=temp;
1471    }
1472  alamin=TOLX/test;
1473  alam=1.0;
1474  for (;;) 
1475    {
1476      for(i=0;i<n;i++) x[i]=fabs(local_xold[i]+alam*p[i]);
1477      /**/     
1478      for(i=0;i<n;i++) 
1479        {
1480          tree->mod->pi[i]=fabs(local_xold[i]+alam*p[i]);
1481/*        if( */
1482/*           (tree->mod->pi[i] < 0.001) || */
1483/*           (tree->mod->pi[i] > 0.999) */
1484/*           ) */
1485/*          break; */
1486        }
1487      /**/
1488      if(i==n) 
1489        {
1490          *f=Return_Abs_Lk(tree);
1491        }
1492      else     *f=1.+fold+ALF*alam*slope;
1493      if (alam < alamin)
1494        {
1495          for (i=0;i<n;i++) x[i]=local_xold[i];
1496          for (i=0;i<n;i++) tree->mod->pi[i]=local_xold[i];
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) 
1523                    {
1524                      disc=b*b-3.0*a*slope;
1525                      if (disc<0.0) tmplam = 0.5*alam;
1526                      else if(b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
1527                      else tmplam = -slope/(b+sqrt(disc));
1528                    }
1529                  else tmplam=(-b+sqrt(disc))/(3.0*a);
1530                }
1531              if (tmplam>0.5*alam) tmplam=0.5*alam;
1532            }
1533        }
1534      alam2=alam;
1535      f2 = *f;
1536      fold2=fold;
1537      alam=MAX(tmplam,0.1*alam);
1538    }
1539  Free(local_xold);
1540}
1541
1542/*********************************************************/
1543
1544/* void Optimize_Global_Rate(arbre *tree) */
1545/* { */
1546/*     printf("\n. Global rate (%f->)",tree->tot_loglk); */
1547/*     Optimize_Single_Param_Generic(tree,&(tree->tbl),tree->tbl,BL_MIN,1.E+4,100); */
1548/*     printf("%f)\n",tree->tot_loglk); */
1549/* } */
1550
1551
1552#undef ALF
1553#undef TOLX
1554#undef NRANSI
1555
1556
Note: See TracBrowser for help on using the repository browser.