1 /*
2 PhyML : a program that computes maximum likelihood phylogenies from
3 DNA or AA homologous sequences
4 Copyright (C) Stephane Guindon. Oct 2003 onward
5 All parts of the source except where indicated are distributed under
6 the GNU public licence. See http://www.opensource.org for details.
7 */
8
9 #include "optimiz.h"
10
11 static phydbl Br_Len_Spline(phydbl *l, t_edge *b, int n_iter_max, phydbl tol, t_tree *tree);
12
13 //////////////////////////////////////////////////////////////
14 //////////////////////////////////////////////////////////////
15
Optimize_Single_Param_Generic(t_tree * tree,phydbl * param,phydbl lim_inf,phydbl lim_sup,phydbl tol,int n_max_iter,int quickdirty)16 void Optimize_Single_Param_Generic(t_tree *tree, phydbl *param, phydbl lim_inf, phydbl lim_sup, phydbl tol, int n_max_iter, int quickdirty)
17 {
18 phydbl lk_init;
19
20 lk_init = tree->c_lnL;
21
22 Generic_Brent_Lk(param,
23 lim_inf,
24 lim_sup,
25 tol,
26 n_max_iter,
27 quickdirty,
28 Wrap_Lk,
29 NULL,
30 tree,
31 NULL,
32 NO);
33
34
35 if(tree->c_lnL < lk_init - tree->mod->s_opt->min_diff_lk_global)
36 {
37 PhyML_Fprintf(stderr,"\n. %.10f < %.10f --> diff=%.10f param value = %f\n",tree->c_lnL,lk_init,tree->c_lnL-lk_init,*param);
38 assert(FALSE);
39 }
40 }
41
42 //////////////////////////////////////////////////////////////
43 //////////////////////////////////////////////////////////////
44
Generic_Brak(phydbl * param,phydbl * ax,phydbl * bx,phydbl * cx,phydbl * fa,phydbl * fb,phydbl * fc,phydbl lim_inf,phydbl lim_sup,t_tree * tree)45 int Generic_Brak(phydbl *param,
46 phydbl *ax, phydbl *bx, phydbl *cx,
47 phydbl *fa, phydbl *fb, phydbl *fc,
48 phydbl lim_inf, phydbl lim_sup,
49 t_tree *tree)
50 {
51 phydbl ulim,u,r,q,fu,dum;
52
53 u = 0.0;
54 *param = *ax;
55
56 if(*param > lim_sup) *param = lim_sup;
57 if(*param < lim_inf) *param = lim_inf;
58 *fa=-Lk(NULL,tree);
59 *param = *bx;
60 if(*param > lim_sup) *param = lim_sup;
61 if(*param < lim_inf) *param = lim_inf;
62 *fb=-Lk(NULL,tree);
63 if (*fb > *fa) {
64 SHFT(dum,*ax,*bx,dum)
65 SHFT(dum,*fb,*fa,dum)
66 }
67 *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
68 *param = FABS(*cx);
69 if(*param > lim_sup) *param = lim_sup;
70 if(*param < lim_inf) *param = lim_inf;
71 *fc=-Lk(NULL,tree);
72 while (*fb > *fc)
73 {
74
75 if(*ax > lim_sup) *ax = lim_sup;
76 if(*ax < lim_inf) *ax = lim_inf;
77 if(*bx > lim_sup) *bx = lim_sup;
78 if(*bx < lim_inf) *bx = lim_inf;
79 if(*cx > lim_sup) *cx = lim_sup;
80 if(*cx < lim_inf) *cx = lim_inf;
81 if(u > lim_sup) u = lim_sup;
82 if(u < lim_inf) u = lim_inf;
83
84 r=(*bx-*ax)*(*fb-*fc);
85 q=(*bx-*cx)*(*fb-*fa);
86 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
87 (2.0*SIGN(MAX(FABS(q-r),MNBRAK_TINY),q-r));
88 ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
89
90 if ((*bx-u)*(u-*cx) > lim_inf)
91 {
92 *param = FABS(u);
93 if(*param > lim_sup) {*param = u = lim_sup;}
94 if(*param < lim_inf) {*param = u = lim_inf;}
95 fu=-Lk(NULL,tree);
96 if (fu < *fc)
97 {
98 *ax=(*bx);
99 *bx=u;
100 *fa=(*fb);
101 *fb=fu;
102 (*ax)=FABS(*ax);
103 (*bx)=FABS(*bx);
104 (*cx)=FABS(*cx);
105 return(0);
106 }
107 else if (fu > *fb)
108 {
109 *cx=u;
110 *fc=fu;
111 (*ax)=FABS(*ax);
112 (*bx)=FABS(*bx);
113 (*cx)=FABS(*cx);
114 return(0);
115 }
116 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
117 *param = FABS(u);
118 if(*param > lim_sup) {*param = u = lim_sup;}
119 if(*param < lim_inf) {*param = u = lim_inf;}
120 fu=-Lk(NULL,tree);
121 }
122 else if ((*cx-u)*(u-ulim) > lim_inf)
123 {
124 *param = FABS(u);
125 if(*param > lim_sup) {*param = u = lim_sup;}
126 if(*param < lim_inf) {*param = u = lim_inf;}
127 fu=-Lk(NULL,tree);
128 if (fu < *fc)
129 {
130 SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
131 *param = FABS(u);
132 SHFT(*fb,*fc,fu,-Lk(NULL,tree))
133 }
134 }
135 else if ((u-ulim)*(ulim-*cx) >= lim_inf)
136 {
137 u=ulim;
138 *param = FABS(u);
139 if(*param > lim_sup) {*param = u = lim_sup;}
140 if(*param < lim_inf) {*param = u = lim_inf;}
141 fu=-Lk(NULL,tree);
142 }
143 else
144 {
145 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
146 *param = FABS(u);
147 if(*param > lim_sup) {*param = u = lim_sup;}
148 if(*param < lim_inf) {*param = u = lim_inf;}
149 fu=-Lk(NULL,tree);
150 }
151 SHFT(*ax,*bx,*cx,u)
152 SHFT(*fa,*fb,*fc,fu)
153
154
155 }
156 (*ax)=FABS(*ax);
157 (*bx)=FABS(*bx);
158 (*cx)=FABS(*cx);
159 return(0);
160 }
161
162 //////////////////////////////////////////////////////////////
163 //////////////////////////////////////////////////////////////
164
Generic_Brak_Lk(phydbl * param,phydbl * ax,phydbl * bx,phydbl * cx,phydbl * fa,phydbl * fb,phydbl * fc,phydbl min,phydbl max,phydbl (* obj_func)(t_edge *,t_tree *,supert_tree *),t_edge * branch,t_tree * tree,supert_tree * stree)165 int Generic_Brak_Lk(phydbl *param,
166 phydbl *ax, phydbl *bx, phydbl *cx,
167 phydbl *fa, phydbl *fb, phydbl *fc,
168 phydbl min, phydbl max,
169 phydbl (*obj_func)(t_edge *,t_tree *,supert_tree *),
170 t_edge *branch, t_tree *tree, supert_tree *stree)
171 {
172 phydbl ulim,u,r,q,fu,dum;
173
174 *param = *ax;
175 if(*param < min) *param = min;
176 if(*param > max) *param = max;
177 *fa = -(*obj_func)(branch,tree,stree);
178
179 *param = *bx;
180 if(*param < min) *param = min;
181 if(*param > max) *param = max;
182 *fb = -(*obj_func)(branch,tree,stree);
183
184 if (*fb > *fa)
185 {
186 SHFT(dum,*ax,*bx,dum)
187 SHFT(dum,*fb,*fa,dum)
188 }
189
190 *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
191 *param = *cx;
192 if(*param < min) *param = min;
193 if(*param > max) *param = max;
194 *fc = -(*obj_func)(branch,tree,stree);
195
196 /* printf("\nx la: %f lb: %f lc: %f fa: %f fb: %f fc: %f",*ax,*bx,*cx,*fa,*fb,*fc); */
197 while (*fb > *fc)
198 {
199 /* printf("\nx la: %f lb: %f lc: %f fa: %f fb: %f fc: %f",*ax,*bx,*cx,*fa,*fb,*fc); */
200 r=(*bx-*ax)*(*fb-*fc);
201 q=(*bx-*cx)*(*fb-*fa);
202 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
203 (2.0*SIGN(MAX(FABS(q-r),MNBRAK_TINY),q-r));
204 ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
205
206 if((*bx-u)*(u-*cx) > 0.0)
207 {
208 *param = u;
209 if(*param < min) *param = min;
210 if(*param > max) *param = max;
211 fu = -(*obj_func)(branch,tree,stree);
212
213 if (fu < *fc)
214 {
215 *ax=(*bx);
216 *bx=u;
217 *fa=(*fb);
218 *fb=fu;
219 return(0);
220 }
221 else if (fu > *fb)
222 {
223 *cx=u;
224 *fc=fu;
225 return(0);
226 }
227 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
228 *param = u;
229 if(*param < min) *param = min;
230 if(*param > max) *param = max;
231 fu = -(*obj_func)(branch,tree,stree);
232 }
233 else if ((*cx-u)*(u-ulim) > 0.0)
234 {
235 *param = u;
236 if(*param < min) *param = min;
237 if(*param > max) *param = max;
238 fu = -(*obj_func)(branch,tree,stree);
239
240 if (fu < *fc)
241 {
242 SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
243 *param = u;
244 if(*param < min) *param = min;
245 if(*param > max) *param = max;
246 SHFT(*fb,*fc,fu,-(*obj_func)(branch,tree,stree))
247 }
248 }
249 else if ((u-ulim)*(ulim-*cx) >= 0.0)
250 {
251 u=ulim;
252 *param = u;
253 if(*param < min) *param = min;
254 if(*param > max) *param = max;
255 fu = -(*obj_func)(branch,tree,stree);
256 }
257 else
258 {
259 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
260 *param = u;
261 if(*param < min) *param = min;
262 if(*param > max) *param = max;
263 fu = -(*obj_func)(branch,tree,stree);
264 }
265 SHFT(*ax,*bx,*cx,u)
266 SHFT(*fa,*fb,*fc,fu)
267 }
268
269 return(0);
270 }
271
272 //////////////////////////////////////////////////////////////
273 //////////////////////////////////////////////////////////////
274
Generic_Brent(phydbl * param,phydbl ax,phydbl cx,phydbl tol,int n_iter_max,phydbl (* obj_func)(t_tree *),t_tree * tree)275 phydbl Generic_Brent(phydbl *param, phydbl ax, phydbl cx, phydbl tol,
276 int n_iter_max,
277 phydbl (*obj_func)(t_tree *),
278 t_tree *tree)
279 {
280 int iter;
281 phydbl a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
282 phydbl e=0.0;
283 phydbl old_score;
284 phydbl bx = *param;
285 int n_opt_step;
286
287 n_opt_step = 0;
288 d=0.0;
289 a=((ax < cx) ? ax : cx);
290 b=((ax > cx) ? ax : cx);
291 x=w=v=bx;
292 (*param) = bx;
293 fw=fv=fx=fu=(*obj_func)(tree);
294 old_score = fw;
295
296 /* PhyML_Printf("\n. %p init_score=%f fu=%f ax=%f cx=%f param=%f",tree,init_score,fu,ax,cx,*param); */
297
298 for(iter=1;iter<=BRENT_IT_MAX;iter++)
299 {
300 xm=0.5*(a+b);
301 tol2=2.0*(tol1=tol*x+BRENT_ZEPS);
302
303 if((fabs(fu-old_score) < tol && iter > 1) || (iter > n_iter_max - 1))
304 {
305 (*param) = x;
306 fu = (*obj_func)(tree);
307 /* printf("\n. return %f [%f] %d",*param,fu,iter); */
308 return fu;
309 }
310
311 if(fabs(e) > tol1)
312 {
313 r=(x-w)*(fx-fv);
314 q=(x-v)*(fx-fw);
315 p=(x-v)*q-(x-w)*r;
316 q=2.0*(q-r);
317 if(q > 0.0) p = -p;
318 q=fabs(q);
319 etemp=e;
320 e=d;
321 if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
322 {
323 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
324 /* PhyML_Printf("\n. Golden section step"); */
325 }
326 else
327 {
328 d=p/q;
329 u=x+d;
330 if (u-a < tol2 || b-u < tol2) d=SIGN(tol1,xm-x);
331 /* PhyML_Printf("\n. Parabolic step [e=%f]",e); */
332 }
333 }
334 else
335 {
336 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
337 /* PhyML_Printf("\n. Golden section step (default) [e=%f tol1=%f a=%f b=%f d=%f x=%f]",e,tol1,a,b,d,x); */
338 }
339
340 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
341 (*param) = u;
342 n_opt_step++;
343 old_score = fu;
344 fu = (*obj_func)(tree);
345 /* PhyML_Printf("\n. iter=%d/%d param=%f lnL=%f u: %f x: %f d: %f logt: %d",iter,BRENT_IT_MAX,*param,fu,u,x,d,logt); */
346
347 if(fu <= fx)
348 {
349 if(u >= x) a=x; else b=x;
350 SHFT(v,w,x,u)
351 SHFT(fv,fw,fx,fu)
352 }
353 else
354 {
355 if (u < x) a=u; else b=u;
356 if (fu < fw || fabs(w-x) < SMALL)
357 {
358 v=w;
359 w=u;
360 fv=fw;
361 fw=fu;
362 }
363 /* else if (fu <= fv || v == x || v == w) */
364 else if (fu < fv || fabs(v-x) < SMALL || fabs(v-w) < SMALL)
365 {
366 v=u;
367 fv=fu;
368 }
369 }
370 }
371
372 PhyML_Printf("\n. Too many iterations in Generic_Brent !");
373 assert(FALSE);
374 return((*obj_func)(tree));
375 /* Not Reached ?? *param=x; */
376 /* Not Reached ?? return fx; */
377 }
378
379 //////////////////////////////////////////////////////////////
380 //////////////////////////////////////////////////////////////
381
382
RRparam_GTR_Golden(phydbl ax,phydbl bx,phydbl cx,phydbl tol,phydbl * xmin,t_tree * tree,calign * cdata,phydbl * param,int n_iter_max)383 phydbl RRparam_GTR_Golden(phydbl ax, phydbl bx, phydbl cx, phydbl tol,
384 phydbl *xmin, t_tree *tree, calign *cdata, phydbl *param, int n_iter_max)
385 {
386 phydbl f1,f2,x0,x1,x2,x3;
387 int n_iter;
388
389
390 x0=ax;
391 x3=cx;
392 if (FABS(cx-bx) > FABS(bx-ax))
393 {
394 x1=bx;
395 x2=bx+GOLDEN_C*(cx-bx);
396 }
397 else
398 {
399 x2=bx;
400 x1=bx-GOLDEN_C*(bx-ax);
401 }
402 (*param)=x1;
403
404 Lk(NULL,tree);
405 f1=-tree->c_lnL;
406 (*param)=x2;
407
408 Lk(NULL,tree);
409 f2=-tree->c_lnL;
410
411 n_iter = 0;
412 while (FABS(x3-x0) > tol*(FABS(x1)+FABS(x2)))
413 {
414
415 if (f2 < f1)
416 {
417 SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
418 (*param)=x2;
419 Lk(NULL,tree);
420 SHFT2(f1,f2,-tree->c_lnL)
421 }
422 else
423 {
424 SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
425 (*param)=x1;
426 Lk(NULL,tree);
427 SHFT2(f2,f1,-tree->c_lnL)
428 }
429
430 if(n_iter++ > n_iter_max) break;
431
432 /* PhyML_Printf("p=%E %f\n",(*param),tree->c_lnL); */
433 }
434 if (f1 < f2)
435 {
436 *xmin=x1;
437 return f1;
438 }
439 else
440 {
441 *xmin=x2;
442 return f2;
443 }
444 }
445
446 //////////////////////////////////////////////////////////////
447 //////////////////////////////////////////////////////////////
448
449
Br_Len_Golden(phydbl ax,phydbl bx,phydbl cx,phydbl tol,phydbl * xmin,t_edge * b_fcus,t_tree * tree)450 phydbl Br_Len_Golden(phydbl ax, phydbl bx, phydbl cx, phydbl tol,
451 phydbl *xmin, t_edge *b_fcus, t_tree *tree)
452 {
453 phydbl f1,f2,x0,x1,x2,x3;
454
455 x0=ax;
456 x3=cx;
457 if(FABS(cx-bx) > FABS(bx-ax))
458 {
459 x1=bx;
460 x2=bx+GOLDEN_C*(cx-bx);
461 }
462 else
463 {
464 x2=bx;
465 x1=bx-GOLDEN_C*(bx-ax);
466 }
467
468 b_fcus->l->v = x1;
469 f1 = -Lk(b_fcus,tree);
470 b_fcus->l->v = x2;
471 f2 = -Lk(b_fcus,tree);
472
473 while (FABS(x3-x0) > tol*(FABS(x1)+FABS(x2)))
474 {
475 if (f2 < f1)
476 {
477 SHFT3(x0,x1,x2,GOLDEN_R*x1+GOLDEN_C*x3)
478 b_fcus->l->v = x2;
479 SHFT2(f1,f2,-Lk(b_fcus,tree))
480 }
481 else
482 {
483 SHFT3(x3,x2,x1,GOLDEN_R*x2+GOLDEN_C*x0)
484 b_fcus->l->v = x1;
485 SHFT2(f2,f1,-Lk(b_fcus,tree))
486 }
487 }
488
489 if (f1 < f2)
490 {
491 *xmin = x1;
492 return -f1;
493 }
494 else
495 {
496 *xmin = x2;
497 return -f2;
498 }
499 }
500
501 //////////////////////////////////////////////////////////////
502 //////////////////////////////////////////////////////////////
503
504
Br_Len_Brak(phydbl * ax,phydbl * bx,phydbl * cx,phydbl * fa,phydbl * fb,phydbl * fc,t_edge * b_fcus,t_tree * tree)505 int Br_Len_Brak(phydbl *ax, phydbl *bx, phydbl *cx,
506 phydbl *fa, phydbl *fb, phydbl *fc,
507 t_edge *b_fcus, t_tree *tree)
508 {
509 phydbl ulim,u,r,q,fu,dum;
510
511 b_fcus->l->v = *ax;
512 *fa=-Lk(b_fcus,tree);
513 b_fcus->l->v = *bx;
514 *fb=-Lk(b_fcus,tree);
515 if (*fb > *fa) {
516 SHFT(dum,*ax,*bx,dum)
517 SHFT(dum,*fb,*fa,dum)
518 }
519 *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
520 b_fcus->l->v = *cx;
521 *fc=-Lk(b_fcus,tree);
522 while (*fb > *fc + tree->mod->s_opt->min_diff_lk_local)
523 {
524 PhyML_Printf("fb=%f fc=%f\n",*fb,*fc);
525 r=(*bx-*ax)*(*fb-*fc);
526 q=(*bx-*cx)*(*fb-*fa);
527 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
528 (2.0*SIGN(MAX(FABS(q-r),MNBRAK_TINY),q-r));
529 ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
530
531 if ((*bx-u)*(u-*cx) > 0.0)
532 {
533 b_fcus->l->v = u;
534 fu=-Lk(b_fcus,tree);
535 if (fu < *fc)
536 {
537 *ax=(*bx);
538 *bx=u;
539 *fa=(*fb);
540 *fb=fu;
541 /* (*ax)=FABS(*ax); */
542 /* (*bx)=FABS(*bx); */
543 /* (*cx)=FABS(*cx); */
544 return(0);
545 }
546 else if (fu > *fb)
547 {
548 *cx=u;
549 *fc=fu;
550 /* (*ax)=FABS(*ax); */
551 /* (*bx)=FABS(*bx); */
552 /* (*cx)=FABS(*cx); */
553 return(0);
554 }
555 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
556 b_fcus->l->v = u;
557 fu=-Lk(b_fcus,tree);
558 }
559 else if ((*cx-u)*(u-ulim) > 0.0)
560 {
561 b_fcus->l->v = FABS(u);
562 fu=-Lk(b_fcus,tree);
563 if (fu < *fc)
564 {
565 SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
566 b_fcus->l->v = u;
567 SHFT(*fb,*fc,fu,-Lk(b_fcus,tree))
568 }
569 }
570 else if ((u-ulim)*(ulim-*cx) >= 0.0)
571 {
572 u=ulim;
573 b_fcus->l->v = u;
574 fu=-Lk(b_fcus,tree);
575 }
576 else
577 {
578 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
579 b_fcus->l->v = u;
580 fu=-Lk(b_fcus,tree);
581 }
582 SHFT(*ax,*bx,*cx,u)
583 SHFT(*fa,*fb,*fc,fu)
584 }
585 (*ax)=FABS(*ax);
586 (*bx)=FABS(*bx);
587 (*cx)=FABS(*cx);
588 return(0);
589 }
590
591 //////////////////////////////////////////////////////////////
592 //////////////////////////////////////////////////////////////
593
Fast_Br_Len(t_edge * b,t_tree * tree,int approx)594 phydbl Fast_Br_Len(t_edge *b, t_tree *tree, int approx)
595 {
596 phydbl init_min_diff_lk_local = tree->mod->s_opt->min_diff_lk_local;
597
598 tree->mod->s_opt->min_diff_lk_local = 1.E-1;
599
600 if(tree->is_mixt_tree) MIXT_Br_Len_Opt(b,tree);
601 else Br_Len_Opt(&(b->l->v),b,tree);
602
603 tree->mod->s_opt->min_diff_lk_local = init_min_diff_lk_local;
604
605 return tree->c_lnL;
606 }
607
608 //////////////////////////////////////////////////////////////
609 //////////////////////////////////////////////////////////////
610
Br_Len_Opt(phydbl * l,t_edge * b,t_tree * tree)611 phydbl Br_Len_Opt(phydbl *l, t_edge *b, t_tree *tree)
612 {
613 phydbl lk_begin, lk_end;
614
615
616 if(tree->is_mixt_tree == YES && tree->ignore_mixt_info == NO)
617 {
618 MIXT_Br_Len_Opt(b,tree);
619 return tree->c_lnL;
620 }
621
622 if(b->l->onoff == OFF) return tree->c_lnL;
623
624 lk_begin = UNLIKELY;
625 lk_end = UNLIKELY;
626
627 Set_Update_Eigen_Lr(YES,tree);
628 Set_Use_Eigen_Lr(NO,tree);
629
630 lk_begin = Lk(b,tree); /* We can't assume that the log-lk value is up-to-date */
631
632 Set_Update_Eigen_Lr(NO,tree);
633 Set_Use_Eigen_Lr(YES,tree);
634
635 Br_Len_Spline(l,b,tree->mod->s_opt->brent_it_max,tree->mod->s_opt->min_diff_lk_local,tree);
636
637 Update_PMat_At_Given_Edge(b,tree);
638
639 Set_Update_Eigen_Lr(NO,tree);
640 Set_Use_Eigen_Lr(NO,tree);
641
642 /* lk_begin = Lk(b,tree); */
643 /* tree->n_tot_bl_opt += Generic_Brent_Lk(l, */
644 /* tree->mod->l_min, */
645 /* tree->mod->l_max, */
646 /* tree->mod->s_opt->min_diff_lk_local, */
647 /* tree->mod->s_opt->brent_it_max, */
648 /* tree->mod->s_opt->quickdirty, */
649 /* Wrap_Lk_At_Given_Edge, */
650 /* b,tree,NULL,NO); */
651
652 /* printf("\n. b->num: %4d l=%12G lnL: %12G",b->num,b->l->v,tree->c_lnL); */
653
654
655 /* lk_end = Lk(b,tree); /\* We can't assume that the log-lk value is up-to-date *\/ */
656 lk_end = tree->c_lnL;
657
658 if(lk_end < lk_begin - tree->mod->s_opt->min_diff_lk_local)
659 {
660 PhyML_Fprintf(stderr,"\n. lk_beg = %f lk_end = %f",lk_begin, lk_end);
661 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d",__FILE__,__LINE__);
662 Exit("\n");
663 }
664
665 return tree->c_lnL;
666 }
667
668 //////////////////////////////////////////////////////////////
669 //////////////////////////////////////////////////////////////
670
Round_Optimize(t_tree * tree,int n_round_max)671 void Round_Optimize(t_tree *tree, int n_round_max)
672 {
673 int n_round,each,freq;
674 phydbl lk_old, lk_new;
675
676 lk_new = tree->c_lnL;
677 lk_old = UNLIKELY;
678 n_round = 0;
679 each = 0;
680 freq = 1;
681
682 while(n_round < n_round_max)
683 {
684 if(tree->mod->s_opt->opt_bl || tree->mod->s_opt->constrained_br_len) Optimize_Br_Len_Serie(n_round_max,tree);
685
686
687 if((tree->mod->s_opt->opt_bl || tree->mod->s_opt->constrained_br_len) &&
688 (tree->verbose > VL2) &&
689 (tree->io->quiet == NO)) Print_Lk(tree,"[Branch lengths ]");
690
691 if(!each)
692 {
693 each = freq;
694 Optimiz_All_Free_Param(tree,(tree->io->quiet)?(0):(tree->verbose > VL2));
695 }
696
697 lk_new = tree->c_lnL;
698
699 if(lk_new < lk_old - tree->mod->s_opt->min_diff_lk_local)
700 {
701 PhyML_Fprintf(stderr,"\n. lk_new = %f lk_old = %f diff = %f",lk_new,lk_old,lk_new-lk_old);
702 assert(FALSE);
703 }
704
705 if((fabs(lk_new - lk_old) < tree->mod->s_opt->min_diff_lk_local) && (each == freq)) break;
706 /* if(fabs(lk_new - lk_old) < tree->mod->s_opt->min_diff_lk_local) break; */
707 lk_old = lk_new;
708
709 n_round++;
710 each--;
711 }
712 }
713
714 //////////////////////////////////////////////////////////////
715 //////////////////////////////////////////////////////////////
716
Optimize_Br_Len_Serie(int n_max_iter,t_tree * tree)717 void Optimize_Br_Len_Serie(int n_max_iter, t_tree *tree)
718 {
719 phydbl lk_init,lk_end;
720 int iter;
721
722
723 /* if(tree->n_tot_bl_opt == 0) */
724 /* { */
725 /* for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl(1.E-6,tree->a_edges[i]->l); */
726 /* for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Min_Thresh(tree->mod->l_min,tree->a_edges[i]->l); */
727 /* for(int i=0;i<2*tree->n_otu-3;++i) Set_Scalar_Dbl_Max_Thresh(tree->mod->l_max,tree->a_edges[i]->l); */
728 /* } */
729
730
731
732 Set_Both_Sides(NO,tree);
733 Lk(NULL,tree);
734
735 lk_init = tree->c_lnL;
736
737 if(tree->mod->gamma_mgf_bl == YES)
738 {
739 Generic_Brent_Lk(&(tree->mod->l_var_sigma),
740 tree->mod->l_var_min,
741 tree->mod->l_var_max,
742 tree->mod->s_opt->min_diff_lk_local,
743 tree->mod->s_opt->brent_it_max,
744 tree->mod->s_opt->quickdirty,
745 Wrap_Lk,NULL,tree,NULL,NO);
746
747 if(tree->c_lnL < lk_init - tree->mod->s_opt->min_diff_lk_local)
748 {
749 PhyML_Printf("\n. %f -- %f",lk_init,tree->c_lnL);
750 PhyML_Printf("\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
751 }
752
753 if((tree->io->quiet)?(0):(tree->verbose > VL2))
754 {
755 Print_Lk(tree,"[Branch len. var. ]");
756 PhyML_Printf("[%10f]",tree->mod->l_var_sigma);
757 }
758 }
759
760 iter = 0;
761 do
762 {
763 lk_init = tree->c_lnL;
764 if(tree->n_root && tree->ignore_root == NO)
765 {
766 Update_Partial_Lk(tree,tree->n_root->b[1],tree->n_root);
767 Optimize_Br_Len_Serie_Post(tree->n_root,tree->n_root->v[1],tree->n_root->b[1],tree);
768 Update_Partial_Lk(tree,tree->n_root->b[2],tree->n_root);
769 Optimize_Br_Len_Serie_Post(tree->n_root,tree->n_root->v[2],tree->n_root->b[2],tree);
770 }
771 else if(tree->n_root && tree->ignore_root == YES)
772 {
773 Optimize_Br_Len_Serie_Post(tree->e_root->rght,
774 tree->e_root->left,
775 tree->e_root,tree);
776
777 Optimize_Br_Len_Serie_Post(tree->e_root->left,
778 tree->e_root->rght,
779 tree->e_root,tree);
780 }
781 else
782 {
783 Optimize_Br_Len_Serie_Post(tree->a_nodes[tree->tip_root],
784 tree->a_nodes[tree->tip_root]->v[0],
785 tree->a_nodes[tree->tip_root]->b[0],
786 tree);
787 }
788
789 lk_end = tree->c_lnL;
790
791 if(lk_end < lk_init - tree->mod->s_opt->min_diff_lk_local)
792 {
793 PhyML_Fprintf(stderr,"\n. lk_init: %f lk_end: %f",lk_init,lk_end);
794 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
795 }
796
797 iter++;
798
799
800 /* PhyML_Printf("\n. lnL: %f %f %f %d",tree->c_lnL,Rgamma((phydbl)(iter+1),(phydbl)(1./(iter+1))),(phydbl)(iter+1)*(phydbl)(1./(iter+1))*(phydbl)(1./(iter+1)),iter); */
801
802 }
803 /* while(lk_end - lk_init > tree->mod->s_opt->min_diff_lk_local && iter < n_max_iter); */
804 while(iter < 1);
805 }
806
807 /*////////////////////////////////////////////////////////////
808 ////////////////////////////////////////////////////////////*/
809
Optimize_Br_Len_Multiplier(t_tree * mixt_tree,int verbose)810 void Optimize_Br_Len_Multiplier(t_tree *mixt_tree, int verbose)
811 {
812 phydbl lk_init;
813 t_tree *tree;
814
815
816 tree = mixt_tree;
817
818 do
819 {
820 if(tree->mod->s_opt->opt_br_len_mult == YES)
821 {
822 lk_init = Get_Lk(tree);
823 /* Generic_Brent_Lk(&(tree->mod->br_len_mult_unscaled->v), */
824 Generic_Brent_Lk(&(tree->mod->br_len_mult->v),
825 1.E-2,1.E+1,
826 tree->mod->s_opt->min_diff_lk_local,
827 tree->mod->s_opt->brent_it_max,
828 tree->mod->s_opt->quickdirty,
829 Wrap_Lk,NULL,mixt_tree,NULL,NO);
830
831 if(Get_Lk(tree) < lk_init - tree->mod->s_opt->min_diff_lk_local)
832 {
833 PhyML_Fprintf(stderr,"\n. %f -- %f",lk_init,tree->c_lnL);
834 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
835 }
836 }
837 tree = tree->next_mixt;
838 }
839 while(tree);
840
841 tree = mixt_tree;
842 do
843 {
844 if(verbose && tree->mod->s_opt->opt_br_len_mult == YES)
845 {
846 Print_Lk(tree,"[Tree scale ]");
847 PhyML_Printf("[%10f]",tree->mod->br_len_mult->v);
848 }
849
850 tree = tree->next_mixt;
851 }
852 while(tree);
853 }
854
855 /*////////////////////////////////////////////////////////////
856 ////////////////////////////////////////////////////////////*/
857
Optimize_Br_Len_Serie_Post(t_node * a,t_node * d,t_edge * b_fcus,t_tree * tree)858 void Optimize_Br_Len_Serie_Post(t_node *a, t_node *d, t_edge *b_fcus, t_tree *tree)
859 {
860 int i;
861 phydbl lk_init;
862
863
864 lk_init = tree->c_lnL;
865
866 if(tree->mod->s_opt->constrained_br_len == YES)
867 {
868 Generic_Brent_Lk(&(tree->mod->br_len_mult->v),
869 1.E-2,1.E+1,
870 tree->mod->s_opt->min_diff_lk_local,
871 tree->mod->s_opt->brent_it_max,
872 tree->mod->s_opt->quickdirty,
873 Wrap_Lk,NULL,tree,NULL,NO);
874
875 if(tree->c_lnL < lk_init - tree->mod->s_opt->min_diff_lk_local)
876 {
877 PhyML_Fprintf(stderr,"\n. %f -- %f",lk_init,tree->c_lnL);
878 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
879 }
880
881 return;
882 }
883
884 if(tree->io->mod->s_opt->opt_bl == YES) Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
885
886 if(tree->c_lnL < lk_init - tree->mod->s_opt->min_diff_lk_local)
887 {
888 PhyML_Fprintf(stderr,"\n. %f -- %f",lk_init,tree->c_lnL);
889 PhyML_Fprintf(stderr,"\n. Edge: %d",b_fcus->num);
890 PhyML_Fprintf(stderr,"\n. is_mixt_tree: %d",tree->is_mixt_tree);
891 Generic_Exit(__FILE__,__LINE__,__FUNCTION__);
892 }
893
894 if(d->tax) return;
895
896 if(tree->n_root && tree->ignore_root == NO)
897 {
898 for(i=0;i<3;++i)
899 {
900 if(d->v[i] != a && d->b[i] != tree->e_root)
901 {
902 Update_Partial_Lk(tree,d->b[i],d);
903 Optimize_Br_Len_Serie_Post(d,d->v[i],d->b[i],tree);
904 }
905 }
906
907 for(i=0;i<3;++i)
908 if(d->v[i] == a || d->b[i] == tree->e_root)
909 {
910 Update_Partial_Lk(tree,d->b[i],d);
911 if(tree->io->mod->s_opt->opt_bl == YES) Br_Len_Opt(&(d->b[i]->l->v),d->b[i],tree);
912 }
913 }
914 else
915 {
916 // Ok if root exists but requires traversal to be initiated from a node != root
917 for(i=0;i<3;++i)
918 {
919 if(d->v[i] != a)
920 {
921 Update_Partial_Lk(tree,d->b[i],d);
922 Optimize_Br_Len_Serie_Post(d,d->v[i],d->b[i],tree);
923 }
924 }
925 Update_Partial_Lk(tree,b_fcus,d);
926 if(tree->io->mod->s_opt->opt_bl == YES) Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
927 }
928 }
929
930 //////////////////////////////////////////////////////////////
931 //////////////////////////////////////////////////////////////
932
Optimiz_Ext_Br(t_tree * tree)933 void Optimiz_Ext_Br(t_tree *tree)
934 {
935 int i;
936 t_edge *b;
937 phydbl lk_init;
938 scalar_dbl *l_init,*v_init;
939
940 lk_init = tree->c_lnL;
941
942 For(i,2*tree->n_otu-3)
943 {
944 b = tree->a_edges[i];
945 if((b->left->tax) || (b->rght->tax))
946 {
947
948 l_init = Duplicate_Scalar_Dbl(b->l);
949 v_init = Duplicate_Scalar_Dbl(b->l_var);
950
951 Br_Len_Opt(&(b->l->v),b,tree);
952
953 if(b->nni->best_l == NULL)
954 {
955 b->nni->best_l = Duplicate_Scalar_Dbl(b->l);
956 b->nni->best_v = Duplicate_Scalar_Dbl(b->l_var);
957 }
958 else
959 {
960 Copy_Scalar_Dbl(b->l,b->nni->best_l);
961 Copy_Scalar_Dbl(b->l_var,b->nni->best_v);
962 }
963
964 if(b->nni->l0 == NULL)
965 {
966 b->nni->l0 = Duplicate_Scalar_Dbl(b->l);
967 b->nni->v0 = Duplicate_Scalar_Dbl(b->l_var);
968 }
969 else
970 {
971 Copy_Scalar_Dbl(b->l,b->nni->l0);
972 Copy_Scalar_Dbl(b->l_var,b->nni->v0);
973 }
974
975 // Revert of original edge lengths
976 Copy_Scalar_Dbl(l_init,b->l);
977 Copy_Scalar_Dbl(v_init,b->l_var);
978
979 Free_Scalar_Dbl(l_init);
980 Free_Scalar_Dbl(v_init);
981
982 /* ori = b; */
983 /* do */
984 /* { */
985 /* b->nni->best_l->v = b->l->v; */
986 /* b->nni->l0->v = b->l->v; */
987 /* b->nni->best_conf = 0; */
988 /* b->l->v = l_init; */
989 /* b = b->next; */
990 /* } */
991 /* while(b); */
992 /* b = ori; */
993 }
994 }
995 tree->c_lnL = lk_init;
996 }
997
998 //////////////////////////////////////////////////////////////
999 //////////////////////////////////////////////////////////////
1000
Optimiz_All_Free_Param(t_tree * tree,int verbose)1001 void Optimiz_All_Free_Param(t_tree *tree, int verbose)
1002 {
1003 int init_both_sides;
1004
1005 if(!tree) return;
1006
1007 if(tree->mixt_tree && tree->mod->ras->invar == YES) return;
1008
1009 init_both_sides = tree->both_sides;
1010
1011
1012 Set_Both_Sides(NO,tree);
1013 Lk(NULL,tree);
1014
1015 Optimize_RR_Params(tree,verbose);
1016 Optimize_TsTv(tree,verbose);
1017 Optimize_Lambda(tree,verbose);
1018 Optimiz_Alpha_And_Pinv(tree,verbose);
1019 Optimize_Pinv(tree,verbose);
1020 Optimize_Alpha(tree,verbose);
1021 Optimize_State_Freqs(tree,verbose);
1022 Optimize_Rmat_Weights(tree,verbose);
1023 Optimize_Efrq_Weights(tree,verbose);
1024 Optimize_Free_Rate(tree,verbose);
1025 Optimize_Br_Len_Multiplier(tree,verbose);
1026
1027 if(tree->io->print_json_trace == YES) JSON_Tree_Io(tree,tree->io->fp_out_json_trace);
1028
1029 if(tree->mod->use_m4mod)
1030 {
1031 int i;
1032
1033 if(tree->mod->s_opt->opt_cov_delta)
1034 {
1035 Set_Update_Eigen(YES,tree->mod);
1036
1037 Generic_Brent_Lk(&(tree->mod->m4mod->delta),
1038 0.01,10.,
1039 tree->mod->s_opt->min_diff_lk_local,
1040 tree->mod->s_opt->brent_it_max,
1041 tree->mod->s_opt->quickdirty,
1042 Wrap_Lk,NULL,tree,NULL,NO);
1043
1044 if(verbose)
1045 {
1046 Print_Lk(tree,"[Switching param. ]");
1047 PhyML_Printf("[%10f]",tree->mod->m4mod->delta);
1048 }
1049
1050 Set_Update_Eigen(NO,tree->mod);
1051
1052 }
1053
1054
1055 if(tree->mod->s_opt->opt_cov_free_rates)
1056 {
1057 int rcat;
1058
1059 Set_Update_Eigen(YES,tree->mod);
1060
1061 for(rcat=0;rcat<tree->mod->m4mod->n_h;rcat++)
1062 {
1063 /* Optimize_Single_Param_Generic(tree,&(tree->mod->m4mod->multipl_unscaled[rcat]), */
1064 /* .01,10., */
1065 /* tree->mod->s_opt->min_diff_lk_local, */
1066 /* tree->mod->s_opt->brent_it_max, */
1067 /* tree->mod->s_opt->quickdirty); */
1068
1069 Generic_Brent_Lk(&(tree->mod->m4mod->multipl_unscaled[rcat]),
1070 0.1,100.,
1071 tree->mod->s_opt->min_diff_lk_local,
1072 tree->mod->s_opt->brent_it_max,
1073 tree->mod->s_opt->quickdirty,
1074 Wrap_Lk,NULL,tree,NULL,NO);
1075
1076 if(verbose)
1077 {
1078 Print_Lk(tree,"[Rel. subst. rate ]");
1079 PhyML_Printf("[%10f]",tree->mod->m4mod->multipl[rcat]);
1080 }
1081 }
1082
1083 for(rcat=0;rcat<tree->mod->m4mod->n_h;rcat++)
1084 {
1085
1086 /* Optimize_Single_Param_Generic(tree,&(tree->mod->m4mod->h_fq_unscaled[rcat]), */
1087 /* .01,100., */
1088 /* tree->mod->s_opt->min_diff_lk_local, */
1089 /* tree->mod->s_opt->brent_it_max, */
1090 /* tree->mod->s_opt->quickdirty); */
1091
1092 Generic_Brent_Lk(&(tree->mod->m4mod->h_fq_unscaled[rcat]),
1093 0.1,100.,
1094 tree->mod->s_opt->min_diff_lk_local,
1095 tree->mod->s_opt->brent_it_max,
1096 tree->mod->s_opt->quickdirty,
1097 Wrap_Lk,NULL,tree,NULL,NO);
1098
1099
1100 if(verbose)
1101 {
1102 Print_Lk(tree,"[Subst. class freq ]");
1103 PhyML_Printf("[%10f]",tree->mod->m4mod->h_fq[rcat]);
1104 }
1105 }
1106
1107 Set_Update_Eigen(NO,tree->mod);
1108
1109 }
1110
1111 if(tree->mod->s_opt->opt_cov_alpha)
1112 {
1113
1114 Set_Update_Eigen(YES,tree->mod);
1115
1116 /* Optimize_Single_Param_Generic(tree,&(tree->mod->m4mod->ras->alpha), */
1117 /* .01,10., */
1118 /* tree->mod->s_opt->min_diff_lk_local, */
1119 /* tree->mod->s_opt->brent_it_max, */
1120 /* tree->mod->s_opt->quickdirty); */
1121
1122 Generic_Brent_Lk(&(tree->mod->m4mod->alpha),
1123 0.01,10.,
1124 tree->mod->s_opt->min_diff_lk_local,
1125 tree->mod->s_opt->brent_it_max,
1126 tree->mod->s_opt->quickdirty,
1127 Wrap_Lk,NULL,tree,NULL,NO);
1128
1129
1130 if(verbose)
1131 {
1132 Print_Lk(tree,"[Alpha (covarion) ]");
1133 PhyML_Printf("[%10f]",tree->mod->m4mod->alpha);
1134 }
1135
1136 Set_Update_Eigen(NO,tree->mod);
1137
1138 }
1139
1140
1141 /* Substitutions between nucleotides are considered to follow a
1142 GTR model */
1143
1144 if(tree->mod->io->datatype == NT)
1145 {
1146 if(tree->mod->whichmodel == GTR || tree->mod->whichmodel == CUSTOM)
1147 {
1148 Set_Update_Eigen(YES,tree->mod);
1149
1150 int *permut = Permutate(tree->mod->r_mat->n_diff_rr);
1151
1152 for(i=0;i<5;i++)
1153 if(permut[i] != 5)
1154 {
1155 Generic_Brent_Lk(&(tree->mod->m4mod->o_rr[permut[i]]),
1156 1.E-4,1.E+4,
1157 tree->mod->s_opt->min_diff_lk_local,
1158 tree->mod->s_opt->brent_it_max,
1159 tree->mod->s_opt->quickdirty,
1160 Wrap_Lk,NULL,tree,NULL,NO);
1161
1162 }
1163
1164 Free(permut);
1165
1166 if(verbose) Print_Lk(tree,"[GTR parameters ]");
1167
1168 Set_Update_Eigen(NO,tree->mod);
1169 }
1170 }
1171 }
1172
1173 Set_Both_Sides(init_both_sides,tree);
1174
1175 if(tree->both_sides == YES) Lk(NULL,tree); /* Needed to update all partial likelihoods */
1176
1177
1178 /* if(tree->next) Optimiz_All_Free_Param(tree->next,verbose); */
1179 /* else Optimiz_All_Free_Param(tree->next,verbose); */
1180
1181 /* if(tree->nextree) Optimiz_All_Free_Param(tree->nextree,verbose); */
1182 }
1183
1184
1185 #define ITMAX 200
1186 #define EPS 3.0e-8
1187 #define TOLX (4*EPS)
1188 #define STPMX 100.0
1189 static phydbl sqrarg;
1190 #define SQR(a) ((sqrarg=(a)) < SMALL ? 0.0 : sqrarg*sqrarg)
1191
BFGS(t_tree * tree,phydbl * p,int n,phydbl gtol,phydbl difff,phydbl step_size,int logt,int is_positive,phydbl (* func)(t_tree * tree),int (* dfunc)(t_tree * tree,phydbl * param,int n_param,phydbl stepsize,int logt,phydbl (* func)(t_tree * tree),phydbl * derivatives,int is_positive),int (* lnsrch)(t_tree * tree,int n,phydbl * xold,phydbl fold,phydbl * g,phydbl * p,phydbl * x,phydbl * f,phydbl stpmax,int * check,int logt,int is_positive),int * failed)1192 void BFGS(t_tree *tree,
1193 phydbl *p,
1194 int n,
1195 phydbl gtol,
1196 phydbl difff,
1197 phydbl step_size,
1198 int logt,
1199 int is_positive,
1200 phydbl(*func)(t_tree *tree),
1201 int(*dfunc)(t_tree *tree,phydbl *param,int n_param,phydbl stepsize,int logt, phydbl(*func)(t_tree *tree),phydbl *derivatives, int is_positive),
1202 int(*lnsrch)(t_tree *tree, int n, phydbl *xold, phydbl fold, phydbl *g, phydbl *p, phydbl *x,phydbl *f, phydbl stpmax, int *check, int logt, int is_positive),
1203 int *failed)
1204 {
1205
1206 int check,i,its,j;
1207 phydbl den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test,fret;
1208 phydbl *dg,*g,*hdg,**hessin,*pnew,*xi;
1209 phydbl fp_old;
1210 phydbl *init;
1211
1212 hessin = (phydbl **)mCalloc(n,sizeof(phydbl *));
1213 for(i=0;i<n;i++) hessin[i] = (phydbl *)mCalloc(n,sizeof(phydbl));
1214 dg = (phydbl *)mCalloc(n,sizeof(phydbl ));
1215 g = (phydbl *)mCalloc(n,sizeof(phydbl ));
1216 pnew = (phydbl *)mCalloc(n,sizeof(phydbl ));
1217 hdg = (phydbl *)mCalloc(n,sizeof(phydbl ));
1218 xi = (phydbl *)mCalloc(n,sizeof(phydbl ));
1219 init = (phydbl *)mCalloc(n,sizeof(phydbl ));
1220
1221
1222 for(i=0;i<n;i++) init[i] = p[i];
1223
1224
1225 /*! p is log transformed */
1226 if(logt == YES) for(i=0;i<n;i++) p[i] = exp(MIN(1.E+2,p[i]));
1227 fp=(*func)(tree);
1228 if(logt == YES) for(i=0;i<n;i++) p[i] = log(p[i]);
1229
1230 /* PhyML_Printf("\n. ENTER BFGS WITH: %f\n",fp); */
1231
1232 fp_old = fp;
1233
1234 (*dfunc)(tree,p,n,step_size,logt,func,g,is_positive);
1235
1236 /* PhyML_Printf("\n. BFGS step_size: %f",step_size); */
1237
1238 for (i=0;i<n;i++)
1239 {
1240 for (j=0;j<n;j++) hessin[i][j]=0.0;
1241 hessin[i][i]=1.0;
1242 xi[i] = -g[i];
1243 sum += p[i]*p[i];
1244 /* PhyML_Printf("\n. BFGS x[%d]: %f p[i]: %f",i,xi[i],p[i]); */
1245 }
1246
1247 stpmax=STPMX*MAX(SQRT(sum),(phydbl)n);
1248 for(its=1;its<=ITMAX;its++)
1249 {
1250
1251 lnsrch(tree,n,p,fp,g,xi,pnew,&fret,stpmax,&check,logt,is_positive);
1252
1253 /* PhyML_Printf("\n. BFGS: %f stpmax: %f fret: %g\n",tree->c_lnL,stpmax,fret); */
1254
1255 fp_old = fp;
1256 fp = fret;
1257
1258 for (i=0;i<n;i++)
1259 {
1260 xi[i]=pnew[i]-p[i];
1261 p[i]=pnew[i];
1262 }
1263
1264 test=0.0;
1265 for (i=0;i<n;i++)
1266 {
1267 temp=xi[i]/MAX(p[i],1.0);
1268 /* printf("\n. x[i]=%f p[i]=%f",xi[i],p[i]); */
1269 if (temp > test) test=temp;
1270 }
1271
1272 if (test < TOLX || (fabs(fp-fp_old) < difff && its > 1))
1273 {
1274 if(fp > fp_old)
1275 {
1276 for(i=0;i<n;i++) p[i] = init[i];
1277 *failed = YES;
1278 }
1279
1280 if(logt == YES) for(i=0;i<n;i++) p[i] = exp(MIN(1.E+2,p[i]));
1281 (*func)(tree);
1282 if(logt == YES) for(i=0;i<n;i++) p[i] = log(p[i]);
1283
1284 for(i=0;i<n;i++) Free(hessin[i]);
1285 free(hessin);
1286 free(xi);
1287 free(pnew);
1288 free(hdg);
1289 free(g);
1290 free(dg);
1291 free(init);
1292 return;
1293 }
1294
1295 for (i=0;i<n;i++) dg[i]=g[i];
1296
1297 (*dfunc)(tree,p,n,step_size,logt,func,g,is_positive);
1298
1299 test=0.0;
1300 den=MAX(fret,1.0);
1301 for (i=0;i<n;i++)
1302 {
1303 temp=g[i]*MAX(p[i],1.0)/den;
1304 if (temp > test) test=temp;
1305 }
1306 if (test < gtol)
1307 {
1308 *failed = NO;
1309 if(logt == YES) for(i=0;i<n;i++) p[i] = exp(MIN(1.E+2,p[i]));
1310 (*func)(tree);
1311 if(logt == YES) for(i=0;i<n;i++) p[i] = log(p[i]);
1312
1313 for(i=0;i<n;i++) Free(hessin[i]);
1314 free(hessin);
1315 free(xi);
1316 free(pnew);
1317 free(hdg);
1318 free(g);
1319 free(dg);
1320 free(init);
1321 return;
1322 }
1323
1324 for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
1325
1326 for (i=0;i<n;i++)
1327 {
1328 hdg[i]=0.0;
1329 for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
1330 }
1331
1332 fac=fae=sumdg=sumxi=0.0;
1333 for (i=0;i<n;i++)
1334 {
1335 fac += dg[i]*xi[i];
1336 fae += dg[i]*hdg[i];
1337 sumdg += SQR(dg[i]);
1338 sumxi += SQR(xi[i]);
1339 }
1340
1341 if(fac*fac > EPS*sumdg*sumxi)
1342 {
1343 fac=1.0/fac;
1344 fad=1.0/fae;
1345 for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
1346 for (i=0;i<n;i++)
1347 {
1348 for (j=0;j<n;j++)
1349 {
1350 hessin[i][j] += fac*xi[i]*xi[j]
1351 -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
1352 }
1353 }
1354 }
1355 for (i=0;i<n;i++)
1356 {
1357 xi[i]=0.0;
1358 for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
1359 }
1360 }
1361 /* PhyML_Printf("\n. Too many iterations in BFGS...\n"); */
1362 *failed = YES;
1363 for(i=0;i<n;i++) Free(hessin[i]);
1364 free(hessin);
1365 free(xi);
1366 free(pnew);
1367 free(hdg);
1368 free(g);
1369 free(dg);
1370 }
1371
1372 #undef ITMAX
1373 #undef EPS
1374 #undef TOLX
1375 #undef STPMX
1376
1377 //////////////////////////////////////////////////////////////
1378 //////////////////////////////////////////////////////////////
1379
1380
1381 #define ITMAX 2000
1382 #define EPS 3.0e-8
1383 #define TOLX (4*EPS)
1384 #define STPMX 100.0
1385 static phydbl sqrarg;
1386 #define SQR(a) ((sqrarg=(a)) < SMALL ? 0.0 : sqrarg*sqrarg)
1387
BFGS_Nonaligned(t_tree * tree,phydbl ** p,int n,phydbl gtol,phydbl difff,phydbl step_size,int logt,int is_positive,phydbl (* func)(t_tree * tree),int (* dfunc_nonaligned)(t_tree * tree,phydbl ** param,int n_param,phydbl stepsize,int logt,phydbl (* func)(t_tree * tree),phydbl * derivatives,int is_positive),int (* lnsrch_nonaligned)(t_tree * tree,int n,phydbl ** xold,phydbl fold,phydbl * g,phydbl * p,phydbl * x,phydbl * f,phydbl stpmax,int * check,int logt,int is_positive),int * failed)1388 void BFGS_Nonaligned(t_tree *tree,
1389 phydbl **p,
1390 int n,
1391 phydbl gtol,
1392 phydbl difff,
1393 phydbl step_size,
1394 int logt,
1395 int is_positive,
1396 phydbl(*func)(t_tree *tree),
1397 int(*dfunc_nonaligned)(t_tree *tree,phydbl **param,int n_param,phydbl stepsize,int logt,phydbl(*func)(t_tree *tree),phydbl *derivatives, int is_positive),
1398 int(*lnsrch_nonaligned)(t_tree *tree, int n, phydbl **xold, phydbl fold,phydbl *g, phydbl *p, phydbl *x,phydbl *f, phydbl stpmax, int *check, int logt, int is_positive),
1399 int *failed)
1400 {
1401
1402 int check,i,its,j;
1403 phydbl den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test,fret;
1404 phydbl *dg,*g,*hdg,**hessin,*pnew,*xi;
1405 phydbl fp_old;
1406 phydbl *init,*sign;
1407
1408 hessin = (phydbl **)mCalloc(n,sizeof(phydbl *));
1409 for(i=0;i<n;i++) hessin[i] = (phydbl *)mCalloc(n,sizeof(phydbl));
1410 dg = (phydbl *)mCalloc(n,sizeof(phydbl ));
1411 g = (phydbl *)mCalloc(n,sizeof(phydbl ));
1412 pnew = (phydbl *)mCalloc(n,sizeof(phydbl ));
1413 hdg = (phydbl *)mCalloc(n,sizeof(phydbl ));
1414 xi = (phydbl *)mCalloc(n,sizeof(phydbl ));
1415 init = (phydbl *)mCalloc(n,sizeof(phydbl ));
1416 sign = (phydbl *)mCalloc(n,sizeof(phydbl ));
1417
1418 for(i=0;i<n;i++) init[i] = (*(p[i]));
1419
1420 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = exp(MIN(1.E+2,*(p[i])));
1421 fp=(*func)(tree);
1422 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = log(*(p[i]));
1423
1424 fp_old = fp;
1425
1426 /* PhyML_Printf("\n- ENTER BFGS WITH: %f\n",fp); */
1427
1428 (*dfunc_nonaligned)(tree,p,n,step_size,logt,func,g,is_positive);
1429
1430 for (i=0;i<n;i++)
1431 {
1432 for (j=0;j<n;j++) hessin[i][j]=0.0;
1433 hessin[i][i]=1.0;
1434 xi[i] = -g[i];
1435 sum += (*(p[i]))*(*(p[i]));
1436 }
1437
1438
1439 stpmax=STPMX*MAX(SQRT(sum),(phydbl)n);
1440 /* stpmax = 0.01*MAX(SQRT(sum),(phydbl)n); */
1441 for(its=1;its<=ITMAX;its++)
1442 {
1443 lnsrch_nonaligned(tree,n,p,fp,g,xi,pnew,&fret,stpmax,&check,logt,is_positive);
1444
1445 /* PhyML_Printf("\n. BFGS -> %f\n",tree->c_lnL); */
1446
1447 fp_old = fp;
1448 fp = fret;
1449
1450 for (i=0;i<n;i++)
1451 {
1452 xi[i]=pnew[i]-(*(p[i]));
1453 (*(p[i]))=pnew[i];
1454 }
1455
1456 test=0.0;
1457 for (i=0;i<n;i++)
1458 {
1459 temp=xi[i]/MAX(*(p[i]),1.0);
1460 /* printf("\n. x[i]=%G p[i]=%f",xi[i],*(p[i])); */
1461 if (temp > test) test=temp;
1462 }
1463
1464 if (test < TOLX || (FABS(fp_old-fp) < difff && its > 1))
1465 {
1466
1467 if(fp > fp_old)
1468 {
1469 for(i=0;i<n;i++) (*(p[i])) = init[i];
1470 *failed = 1;
1471 }
1472
1473 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = exp(MIN(1.E+2,*(p[i])));
1474 for(i=0;i<n;i++) sign[i] = *(p[i]) > .0 ? 1. : -1.;
1475 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) = FABS(*(p[i]));
1476 (*func)(tree);
1477 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) *= sign[i];
1478 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = log(*(p[i]));
1479
1480 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) = FABS(*(p[i]));
1481
1482 for(i=0;i<n;i++) Free(hessin[i]);
1483 free(hessin);
1484 free(xi);
1485 free(pnew);
1486 free(hdg);
1487 free(g);
1488 free(dg);
1489 free(init);
1490 free(sign);
1491 return;
1492 }
1493
1494 for (i=0;i<n;i++) dg[i]=g[i];
1495
1496 (*dfunc_nonaligned)(tree,p,n,step_size,logt,func,g,is_positive);
1497
1498 test=0.0;
1499 den=MAX(fret,1.0);
1500 for (i=0;i<n;i++)
1501 {
1502 temp=g[i]*MAX(*(p[i]),1.0)/den;
1503 if (temp > test) test=temp;
1504 }
1505
1506 if (test < gtol)
1507 {
1508 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = exp(MIN(1.E+2,*(p[i])));
1509 for(i=0;i<n;i++) sign[i] = *(p[i]) > .0 ? 1. : -1.;
1510 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) = FABS(*(p[i]));
1511 (*func)(tree);
1512 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) *= sign[i];
1513 if(logt == YES) for(i=0;i<n;i++) (*(p[i])) = log(*(p[i]));
1514
1515 if(is_positive == YES) for(i=0;i<n;i++) *(p[i]) = FABS(*(p[i]));
1516
1517 for(i=0;i<n;i++) Free(hessin[i]);
1518 free(hessin);
1519 free(xi);
1520 free(pnew);
1521 free(hdg);
1522 free(g);
1523 free(dg);
1524 free(init);
1525 free(sign);
1526 return;
1527 }
1528
1529 for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
1530
1531 for (i=0;i<n;i++)
1532 {
1533 hdg[i]=0.0;
1534 for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
1535 }
1536
1537 fac=fae=sumdg=sumxi=0.0;
1538 for (i=0;i<n;i++)
1539 {
1540 fac += dg[i]*xi[i];
1541 fae += dg[i]*hdg[i];
1542 sumdg += SQR(dg[i]);
1543 sumxi += SQR(xi[i]);
1544 }
1545
1546 if(fac*fac > EPS*sumdg*sumxi)
1547 {
1548 fac=1.0/fac;
1549 fad=1.0/fae;
1550 for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
1551 for (i=0;i<n;i++)
1552 {
1553 for (j=0;j<n;j++)
1554 {
1555 hessin[i][j] += fac*xi[i]*xi[j]
1556 -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
1557 }
1558 }
1559 }
1560 for (i=0;i<n;i++)
1561 {
1562 xi[i]=0.0;
1563 for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
1564 }
1565 }
1566 PhyML_Printf("\n. Too many iterations in BFGS...\n");
1567 *failed = YES;
1568 for(i=0;i<n;i++) Free(hessin[i]);
1569 free(hessin);
1570 free(xi);
1571 free(pnew);
1572 free(hdg);
1573 free(g);
1574 free(dg);
1575 free(sign);
1576 }
1577
1578 #undef ITMAX
1579 #undef EPS
1580 #undef TOLX
1581 #undef STPMX
1582
1583 //////////////////////////////////////////////////////////////
1584 //////////////////////////////////////////////////////////////
1585
1586 #define ALF 1.0e-4
1587 #define TOLX 1.0e-7
1588
Lnsrch(t_tree * tree,int n,phydbl * xold,phydbl fold,phydbl * g,phydbl * p,phydbl * x,phydbl * f,phydbl stpmax,int * check,int logt,int is_positive)1589 int Lnsrch(t_tree *tree, int n, phydbl *xold, phydbl fold,
1590 phydbl *g, phydbl *p, phydbl *x,
1591 phydbl *f, phydbl stpmax, int *check, int logt, int is_positive)
1592 {
1593 int i;
1594 phydbl a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1595 phydbl *local_xold,*sign;
1596
1597 alam = alam2 = f2 = fold2 = tmplam = .0;
1598
1599 local_xold = (phydbl *)mCalloc(n,sizeof(phydbl));
1600 sign = (phydbl *)mCalloc(n,sizeof(phydbl));
1601
1602 for(i=0;i<n;i++) local_xold[i] = xold[i];
1603
1604 *check=0;
1605 for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1606 sum=SQRT(sum);
1607 /* PhyML_Printf("\n. lnsrch sum: %f",sum); */
1608 if(sum > stpmax) for(i=0;i<n;i++) p[i] *= stpmax/sum;
1609 /* for(i=0;i<n;i++) PhyML_Printf("\n. lnsrch p[i]: %f",p[i]); */
1610 slope=0.0;
1611 for(i=0;i<n;i++) slope += g[i]*p[i];
1612 /* PhyML_Printf("\n. lnsrch slope: %f",slope); */
1613 test=0.0;
1614 for(i=0;i<n;i++)
1615 {
1616 temp=p[i]/MAX(local_xold[i],1.0);
1617 if (temp > test) test=temp;
1618 }
1619 alamin=TOLX/test;
1620 alam=1.0;
1621 for (;;)
1622 {
1623 for(i=0;i<n;i++)
1624 {
1625 x[i]=local_xold[i]+alam*p[i];
1626 xold[i] = x[i];
1627 /* PhyML_Printf("\n. lnsrch x[i]: %f",x[i]); */
1628 }
1629
1630 /* PhyML_Printf("\n. lnsrch loop slope: %f alam: %f alam2: %f",slope,alam,alam2); */
1631
1632 if(i==n)
1633 {
1634 if(logt == YES) for(i=0;i<n;i++) xold[i] = exp(MIN(1.E+2,xold[i]));
1635 for(i=0;i<n;i++) sign[i] = xold[i] < .0 ? -1. : 1.;
1636 if(is_positive == YES) for(i=0;i<n;i++) xold[i] = FABS(xold[i]);
1637 /* for(i=0;i<n;i++) PhyML_Printf("\n. <<>> %f",xold[i]); */
1638 *f=Return_Abs_Lk(tree);
1639 if(is_positive == YES) for(i=0;i<n;i++) xold[i] *= sign[i];
1640 if(logt == YES) for(i=0;i<n;i++) xold[i] = log(xold[i]);
1641 }
1642 else *f=1.+fold+ALF*alam*slope;
1643 if (alam < alamin)
1644 {
1645 *check=1;
1646 for(i=0;i<n;i++) xold[i] = local_xold[i];
1647 if(is_positive == YES) for(i=0;i<n;i++) xold[i] = FABS(xold[i]);
1648 Free(local_xold);
1649 Free(sign);
1650 return 0;
1651 }
1652 else if (*f <= fold+ALF*alam*slope)
1653 {
1654 for(i=0;i<n;i++) xold[i] = local_xold[i];
1655 if(is_positive == YES) for(i=0;i<n;i++) xold[i] = FABS(xold[i]);
1656 Free(local_xold);
1657 Free(sign);
1658 return 0;
1659 }
1660 else
1661 {
1662 /* if (alam == 1.0) */
1663 if ((alam < 1.0+SMALL) && (alam > 1.0-SMALL))
1664 tmplam = -slope/(2.0*(*f-fold-slope));
1665 else
1666 {
1667 rhs1 = *f-fold-alam*slope;
1668 rhs2=f2-fold2-alam2*slope;
1669 a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
1670 b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
1671 if (a < SMALL && a > -SMALL) tmplam = -slope/(2.0*b);
1672 else
1673 {
1674 disc=b*b-3.0*a*slope;
1675 if (disc<0.0) tmplam = 0.5*alam;
1676 else if(b <= 0.0) tmplam=(-b+SQRT(disc))/(3.0*a);
1677 else tmplam = -slope/(b+SQRT(disc));
1678 }
1679 if (tmplam>0.5*alam) tmplam=0.5*alam;
1680 }
1681 }
1682 alam2=alam;
1683 f2 = *f;
1684 fold2=fold;
1685 alam=MAX(tmplam,0.1*alam);
1686 }
1687 Free(sign);
1688 Free(local_xold);
1689 return 1;
1690 }
1691
1692 #undef ALF
1693 #undef TOLX
1694 #undef NRANSI
1695
1696 //////////////////////////////////////////////////////////////
1697 //////////////////////////////////////////////////////////////
1698
1699 #define ALF 1.0e-4
1700 #define TOLX 1.0e-7
1701
Lnsrch_Nonaligned(t_tree * tree,int n,phydbl ** xold,phydbl fold,phydbl * g,phydbl * p,phydbl * x,phydbl * f,phydbl stpmax,int * check,int logt,int is_positive)1702 int Lnsrch_Nonaligned(t_tree *tree, int n, phydbl **xold, phydbl fold,
1703 phydbl *g, phydbl *p, phydbl *x,
1704 phydbl *f, phydbl stpmax, int *check, int logt, int is_positive)
1705 {
1706 int i;
1707 phydbl a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
1708 phydbl *local_xold,*sign;
1709
1710
1711 alam = alam2 = f2 = fold2 = tmplam = .0;
1712
1713 local_xold = (phydbl *)mCalloc(n,sizeof(phydbl));
1714 sign = (phydbl *)mCalloc(n,sizeof(phydbl));
1715
1716 for(i=0;i<n;i++) local_xold[i] = *(xold[i]);
1717
1718 *check=0;
1719 for(sum=0.0,i=0;i<n;i++) sum += p[i]*p[i];
1720 sum=SQRT(sum);
1721 if(sum > stpmax)
1722 for(i=0;i<n;i++) p[i] *= stpmax/sum;
1723 for(slope=0.0,i=0;i<n;i++)
1724 slope += g[i]*p[i];
1725 test=0.0;
1726 for(i=0;i<n;i++)
1727 {
1728 temp=p[i]/MAX(local_xold[i],1.0);
1729 if (temp > test) test=temp;
1730 }
1731 alamin=TOLX/test;
1732 alam=1.0;
1733 for (;;)
1734 {
1735 for(i=0;i<n;i++)
1736 {
1737 x[i]=local_xold[i]+alam*p[i];
1738 *(xold[i]) = x[i];
1739 }
1740
1741 if(i==n)
1742 {
1743 if(logt == YES) for(i=0;i<n;i++) *(xold[i]) = exp(MIN(1.E+2,*(xold[i])));
1744 for(i=0;i<n;i++) sign[i] = *(xold[i]) < .0 ? -1. : 1.;
1745 if(is_positive == YES) for(i=0;i<n;i++) *(xold[i]) = FABS(*(xold[i]));
1746 *f=Return_Abs_Lk(tree);
1747 if(is_positive == YES) for(i=0;i<n;i++) *(xold[i]) *= sign[i];
1748 if(logt == YES) for(i=0;i<n;i++) *(xold[i]) = log(*(xold[i]));
1749 }
1750 else *f=1.+fold+ALF*alam*slope;
1751
1752 if (alam < alamin)
1753 {
1754 *check=1;
1755 for(i=0;i<n;i++) *(xold[i]) = local_xold[i];
1756 if(is_positive == YES) for(i=0;i<n;i++) *(xold[i]) = FABS(*(xold[i]));
1757 Free(local_xold);
1758 Free(sign);
1759 return 0;
1760 }
1761 else if (*f <= fold+ALF*alam*slope)
1762 {
1763 for(i=0;i<n;i++) *(xold[i]) = local_xold[i];
1764 if(is_positive) for(i=0;i<n;i++) *(xold[i]) = FABS(*(xold[i]));
1765 Free(local_xold);
1766 Free(sign);
1767 return 0;
1768 }
1769 else
1770 {
1771 /* if (alam == 1.0) */
1772 if ((alam < 1.0+SMALL) && (alam > 1.0-SMALL))
1773 tmplam = -slope/(2.0*(*f-fold-slope));
1774 else
1775 {
1776 rhs1 = *f-fold-alam*slope;
1777 rhs2=f2-fold2-alam2*slope;
1778 a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
1779 b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
1780 if (a < SMALL && a > -SMALL) tmplam = -slope/(2.0*b);
1781 else
1782 {
1783 disc=b*b-3.0*a*slope;
1784 if (disc<0.0) tmplam = 0.5*alam;
1785 else if(b <= 0.0) tmplam=(-b+SQRT(disc))/(3.0*a);
1786 else tmplam = -slope/(b+SQRT(disc));
1787 }
1788 if (tmplam>0.5*alam) tmplam=0.5*alam;
1789 }
1790 }
1791 alam2=alam;
1792 f2 = *f;
1793 fold2=fold;
1794 alam=MAX(tmplam,0.1*alam);
1795 }
1796 Free(local_xold);
1797 Free(sign);
1798 return 1;
1799 }
1800
1801 #undef ALF
1802 #undef TOLX
1803 #undef NRANSI
1804
1805 //////////////////////////////////////////////////////////////
1806 //////////////////////////////////////////////////////////////
1807
Dist_F_Brak(phydbl * ax,phydbl * bx,phydbl * cx,phydbl * F,phydbl * param,t_mod * mod)1808 int Dist_F_Brak(phydbl *ax, phydbl *bx, phydbl *cx, phydbl *F, phydbl *param, t_mod *mod)
1809 {
1810 phydbl ulim,u,r,q,dum;
1811 phydbl fa, fb, fc, fu;
1812
1813 fa = -Lk_Dist(F,FABS(*ax),mod);
1814 fb = -Lk_Dist(F,FABS(*bx),mod);
1815
1816 if(fb > fa)
1817 {
1818 SHFT(dum,*ax,*bx,dum)
1819 SHFT(dum,fb,fa,dum)
1820 }
1821
1822 *cx=(*bx)+MNBRAK_GOLD*(*bx-*ax);
1823 fc = -Lk_Dist(F,FABS(*cx),mod);
1824
1825 while (fb > fc)
1826 {
1827 r=(*bx-*ax)*(fb-fc);
1828 q=(*bx-*cx)*(fb-fa);
1829 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
1830 (2.0*SIGN(MAX(FABS(q-r),MNBRAK_TINY),q-r));
1831 ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
1832
1833 if ((*bx-u)*(u-*cx) > 0.0)
1834 {
1835 fu = -Lk_Dist(F,FABS(u),mod);
1836 if (fu < fc)
1837 {
1838 *ax=(*bx);
1839 *bx=u;
1840 fa=fb;
1841 fb=fu;
1842 return(0);
1843 }
1844 else if (fu > fb)
1845 {
1846 *cx=u;
1847 fc=fu;
1848 return(0);
1849 }
1850 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
1851 fu = -Lk_Dist(F,FABS(u),mod);
1852 }
1853 else if ((*cx-u)*(u-ulim) > 0.0)
1854 {
1855 fu = -Lk_Dist(F,FABS(u),mod);
1856 if (fu < fc)
1857 {
1858 SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
1859 SHFT(fb,fc,fu,-Lk_Dist(F,FABS(u),mod))
1860 }
1861 }
1862 else if ((u-ulim)*(ulim-*cx) >= 0.0)
1863 {
1864 u = ulim;
1865 fu = -Lk_Dist(F,FABS(u),mod);
1866 }
1867 else
1868 {
1869 u =(*cx)+MNBRAK_GOLD*(*cx-*bx);
1870 fu = -Lk_Dist(F,FABS(u),mod);
1871 }
1872
1873 SHFT(*ax,*bx,*cx,u)
1874 SHFT(fa,fb,fc,fu)
1875 }
1876 return(0);
1877 }
1878
1879 //////////////////////////////////////////////////////////////
1880 //////////////////////////////////////////////////////////////
1881
Dist_F_Brent(phydbl ax,phydbl bx,phydbl cx,phydbl tol,int n_iter_max,phydbl * param,phydbl * F,t_mod * mod)1882 phydbl Dist_F_Brent(phydbl ax, phydbl bx, phydbl cx, phydbl tol, int n_iter_max,
1883 phydbl *param, phydbl *F, t_mod *mod)
1884 {
1885 int iter;
1886 phydbl a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
1887 phydbl e=0.0;
1888 phydbl old_lnL,init_lnL, curr_lnL;
1889
1890 d=0.0;
1891 a=((ax < cx) ? ax : cx);
1892 b=((ax > cx) ? ax : cx);
1893 x = w = v = bx;
1894 old_lnL = UNLIKELY;
1895 fw = fv = fx = -Lk_Dist(F,fabs(bx),mod);
1896 curr_lnL = init_lnL = -fw;
1897
1898 /* PhyML_Printf("\n. bx=%f f: %f %f %f %f fx: %f",bx,mod->e_frq->pi->v[0],mod->e_frq->pi->v[1],mod->e_frq->pi->v[2],mod->e_frq->pi->v[3],fx); */
1899 assert(isnan(fx) == FALSE);
1900 assert(isinf(fx) == FALSE);
1901
1902 for(iter=1;iter<=BRENT_IT_MAX;iter++)
1903 {
1904 xm=0.5*(a+b);
1905
1906 tol2=2.0*(tol1=tol*FABS(x)+BRENT_ZEPS);
1907
1908 if(
1909 ((FABS(curr_lnL-old_lnL) < mod->s_opt->min_diff_lk_local) &&
1910 (curr_lnL > init_lnL - mod->s_opt->min_diff_lk_local)) ||
1911 (iter > n_iter_max - 1)
1912 )
1913 {
1914 *param = x;
1915 curr_lnL = Lk_Dist(F,*param,mod);
1916 return -curr_lnL;
1917 }
1918
1919 if(FABS(e) > tol1)
1920 {
1921 r=(x-w)*(fx-fv);
1922 q=(x-v)*(fx-fw);
1923 p=(x-v)*q-(x-w)*r;
1924 q=2.0*(q-r);
1925 if(q > 0.0) p = -p;
1926 q=FABS(q);
1927 etemp=e;
1928 e=d;
1929 if(FABS(p) >= FABS(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
1930 {
1931 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
1932 /* PhyML_Printf("Golden section step\n"); */
1933 }
1934 else
1935 {
1936 d=p/q;
1937 u=x+d;
1938 if (u-a < tol2 || b-u < tol2)
1939 d=SIGN(tol1,xm-x);
1940 /* PhyML_Printf("Parabolic step\n"); */
1941 }
1942 }
1943 else
1944 {
1945 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
1946 /* PhyML_Printf("Golden section step (default)\n"); */
1947 }
1948
1949 u=(FABS(d) >= tol1 ? x+d : x+SIGN(tol1,d));
1950 (*param) = FABS(u);
1951 old_lnL = curr_lnL;
1952 fu = -Lk_Dist(F,FABS(u),mod);
1953 curr_lnL = -fu;
1954 /* PhyML_Printf("param=%f loglk=%f\n",*param,fu); */
1955
1956 /* if(fu <= fx) */
1957 if(fu < fx)
1958 {
1959 if(iter > n_iter_max) return -fu;
1960
1961 if(u >= x) a=x; else b=x;
1962 SHFT(v,w,x,u)
1963 SHFT(fv,fw,fx,fu)
1964 }
1965 else
1966 {
1967 if (u < x) a=u; else b=u;
1968 /* if (fu <= fw || w == x) */
1969 if (fu < fw || FABS(w-x) < SMALL)
1970 {
1971 v=w;
1972 w=u;
1973 fv=fw;
1974 fw=fu;
1975 }
1976 /* else if (fu <= fv || v == x || v == w) */
1977 else if (fu < fv || FABS(v-x) < SMALL || FABS(v-w) < SMALL)
1978 {
1979 v=u;
1980 fv=fu;
1981 }
1982 }
1983 }
1984
1985 Exit("\n. Too many iterations in Dist_F_Brent !\n");
1986 return(-1);
1987 }
1988
1989 //////////////////////////////////////////////////////////////
1990 //////////////////////////////////////////////////////////////
1991
Opt_Dist_F(phydbl * dist,phydbl * F,t_mod * mod)1992 void Opt_Dist_F(phydbl *dist, phydbl *F, t_mod *mod)
1993 {
1994 phydbl ax,bx,cx;
1995
1996 if(*dist < mod->l_min) *dist = mod->l_min;
1997
1998 ax = mod->l_min;
1999 bx = (*dist);
2000 cx = mod->l_max;
2001
2002 /* PhyML_Printf("\n. bx: %g",bx); */
2003
2004 /* Dist_F_Brak(&ax,&bx,&cx,F,dist,mod); */
2005 Dist_F_Brent(ax,bx,cx,1.E-10,1000,dist,F,mod);
2006 }
2007
2008 //////////////////////////////////////////////////////////////
2009 //////////////////////////////////////////////////////////////
2010
2011
Missing_Dist_Brak(phydbl * ax,phydbl * bx,phydbl * cx,int x,int y,matrix * mat)2012 int Missing_Dist_Brak(phydbl *ax, phydbl *bx, phydbl *cx, int x, int y, matrix *mat)
2013 {
2014 phydbl ulim,u,r,q,dum;
2015 phydbl fa, fb, fc, fu;
2016
2017 fa = Least_Square_Missing_Dist_XY(x,y,FABS(*ax),mat);
2018 fb = Least_Square_Missing_Dist_XY(x,y,FABS(*bx),mat);
2019
2020 if(fb > fa)
2021 {
2022 SHFT(dum,*ax,*bx,dum)
2023 SHFT(dum,fb,fa,dum)
2024 }
2025
2026 *cx=(*bx)+MNBRAK_GOLD*((*bx)-(*ax));
2027 fc = Least_Square_Missing_Dist_XY(x,y,FABS(*cx),mat);
2028
2029 while (fb > fc)
2030 {
2031 r=((*bx)-(*ax))*(fb-fc);
2032 q=((*bx)-(*cx))*(fb-fa);
2033 u=(*bx)-(((*bx)-(*cx))*q-((*bx)-(*ax))*r)/
2034 (2.0*SIGN(MAX(FABS(q-r),MNBRAK_TINY),q-r));
2035 ulim=(*bx)+MNBRAK_GLIMIT*(*cx-*bx);
2036
2037 if ((*bx-u)*(u-*cx) > 0.0)
2038 {
2039 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2040 if (fu < fc)
2041 {
2042 *ax=(*bx);
2043 *bx=u;
2044 fa=fb;
2045 fb=fu;
2046 return(0);
2047 }
2048 else if (fu > fb)
2049 {
2050 *cx=u;
2051 fc=fu;
2052 return(0);
2053 }
2054 u=(*cx)+MNBRAK_GOLD*(*cx-*bx);
2055 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2056 }
2057 else if ((*cx-u)*(u-ulim) > 0.0)
2058 {
2059 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2060 if (fu < fc)
2061 {
2062 SHFT(*bx,*cx,u,*cx+MNBRAK_GOLD*(*cx-*bx))
2063 SHFT(fb,fc,fu,Least_Square_Missing_Dist_XY(x,y,FABS(u),mat))
2064 }
2065 }
2066 else if ((u-ulim)*(ulim-*cx) >= 0.0)
2067 {
2068 u = ulim;
2069 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2070 }
2071 else
2072 {
2073 u =(*cx)+MNBRAK_GOLD*(*cx-*bx);
2074 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2075 }
2076
2077 SHFT(*ax,*bx,*cx,u)
2078 SHFT(fa,fb,fc,fu)
2079 }
2080 return(0);
2081 }
2082
2083 //////////////////////////////////////////////////////////////
2084 //////////////////////////////////////////////////////////////
2085
2086
Missing_Dist_Brent(phydbl ax,phydbl bx,phydbl cx,phydbl tol,int n_iter_max,int x,int y,matrix * mat)2087 phydbl Missing_Dist_Brent(phydbl ax, phydbl bx, phydbl cx, phydbl tol, int n_iter_max,
2088 int x, int y, matrix *mat)
2089 {
2090 int iter;
2091 phydbl a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,xx,xm;
2092 phydbl e=0.0;
2093
2094 d=0.0;
2095 a=((ax < cx) ? ax : cx);
2096 b=((ax > cx) ? ax : cx);
2097 xx=w=v=bx;
2098 fx=Least_Square_Missing_Dist_XY(x,y,FABS(bx),mat);
2099 fw=fv=-fx;
2100
2101 for(iter=1;iter<=BRENT_IT_MAX;iter++)
2102 {
2103 xm=0.5*(a+b);
2104 tol2=2.0*(tol1=tol*FABS(xx)+BRENT_ZEPS);
2105
2106 if(FABS(xx-xm) <= (tol2-0.5*(b-a)))
2107 {
2108 mat->dist[x][y] = xx;
2109 Least_Square_Missing_Dist_XY(x,y,mat->dist[x][y],mat);
2110 return -fx;
2111 }
2112
2113 if(FABS(e) > tol1)
2114 {
2115 r=(xx-w)*(fx-fv);
2116 q=(xx-v)*(fx-fw);
2117 p=(xx-v)*q-(xx-w)*r;
2118 q=2.0*(q-r);
2119 if(q > 0.0) p = -p;
2120 q=FABS(q);
2121 etemp=e;
2122 e=d;
2123 if(FABS(p) >= FABS(0.5*q*etemp) || p <= q*(a-xx) || p >= q*(b-xx))
2124 {
2125 d=BRENT_CGOLD*(e=(xx >= xm ? a-xx : b-xx));
2126 /* PhyML_Printf("Golden section step\n"); */
2127 }
2128 else
2129 {
2130 d=p/q;
2131 u=xx+d;
2132 if (u-a < tol2 || b-u < tol2)
2133 d=SIGN(tol1,xm-xx);
2134 /* PhyML_Printf("Parabolic step\n"); */
2135 }
2136 }
2137 else
2138 {
2139 d=BRENT_CGOLD*(e=(xx >= xm ? a-xx : b-xx));
2140 /* PhyML_Printf("Golden section step (default)\n"); */
2141 }
2142
2143 u=(FABS(d) >= tol1 ? xx+d : xx+SIGN(tol1,d));
2144 fu = Least_Square_Missing_Dist_XY(x,y,FABS(u),mat);
2145
2146 /* PhyML_Printf("param=%f loglk=%f\n",u,fu); */
2147
2148 /* if(fu <= fx) */
2149 if(fu < fx)
2150 {
2151 if(iter > n_iter_max) return -fu;
2152
2153 if(u >= xx) a=xx; else b=xx;
2154 SHFT(v,w,xx,u)
2155 SHFT(fv,fw,fx,fu)
2156 }
2157 else
2158 {
2159 if (u < xx) a=u; else b=u;
2160 /* if (fu <= fw || w == xx) */
2161 if (fu < fw || FABS(w-xx) < SMALL)
2162 {
2163 v=w;
2164 w=u;
2165 fv=fw;
2166 fw=fu;
2167 }
2168 /* else if (fu <= fv || v == xx || v == w) */
2169 else if (fu < fv || FABS(v-xx) < SMALL || FABS(v-w) < SMALL)
2170 {
2171 v=u;
2172 fv=fu;
2173 }
2174 }
2175 }
2176 Exit("\n. Too many iterations in Missing_Dist_Brent !");
2177 return(-1);
2178 }
2179
2180 //////////////////////////////////////////////////////////////
2181 //////////////////////////////////////////////////////////////
2182
2183
Opt_Missing_Dist(int x,int y,matrix * mat)2184 void Opt_Missing_Dist(int x, int y, matrix *mat)
2185 {
2186 phydbl ax,bx,cx;
2187
2188 ax = DIST_MAX;
2189 bx = DIST_MAX/4.;
2190
2191 Missing_Dist_Brak(&ax,&bx,&cx,x,y,mat);
2192 PhyML_Printf("ax=%f bx=%f cx=%f\n",FABS(ax),FABS(bx),FABS(cx));
2193 Missing_Dist_Brent(FABS(ax),FABS(bx),FABS(cx),1.E-5,100,x,y,mat);
2194 }
2195
2196 //////////////////////////////////////////////////////////////
2197 //////////////////////////////////////////////////////////////
2198
Optimiz_Alpha_And_Pinv(t_tree * mixt_tree,int verbose)2199 int Optimiz_Alpha_And_Pinv(t_tree *mixt_tree, int verbose)
2200 {
2201 scalar_dbl **alpha;
2202 int n_alpha;
2203 t_tree *tree;
2204 int i;
2205
2206 Set_Update_Eigen(NO,mixt_tree->mod);
2207
2208 alpha = NULL;
2209 n_alpha = 0;
2210 tree = mixt_tree;
2211
2212 do
2213 {
2214 if(tree->mod->s_opt->opt_alpha == YES && tree->mod->ras->n_catg > 1 && tree->mod->s_opt->opt_pinvar == YES)
2215 {
2216 for(i=0;i<n_alpha;i++) if(tree->mod->ras->alpha == alpha[i]) break;
2217
2218 if(i == n_alpha)
2219 {
2220 if(!alpha) alpha = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
2221 else alpha = (scalar_dbl **)mRealloc(alpha,n_alpha+1,sizeof(scalar_dbl *));
2222 alpha[n_alpha] = tree->mod->ras->alpha;
2223 n_alpha++;
2224
2225 if(tree->mod->s_opt->opt_alpha == YES &&
2226 tree->mod->ras->free_mixt_rates == NO)
2227 {
2228 if(tree->mod->ras->n_catg > 1)
2229 {
2230 Generic_Brent_Lk(&(tree->mod->ras->alpha->v),
2231 ALPHA_MIN,ALPHA_MAX,
2232 tree->mod->s_opt->min_diff_lk_local,
2233 tree->mod->s_opt->brent_it_max,
2234 tree->mod->s_opt->quickdirty,
2235 Wrap_Lk,NULL,mixt_tree,NULL,NO);
2236 }
2237 if(verbose == YES)
2238 {
2239 Print_Lk(mixt_tree,"[Alpha ]");
2240 PhyML_Printf("[%10f]",tree->mod->ras->alpha->v);
2241 }
2242 }
2243
2244 if(tree->mod->s_opt->opt_pinvar == YES &&
2245 tree->mod->ras->free_mixt_rates == NO)
2246 {
2247 tree->mod->s_opt->skip_tree_traversal = YES;
2248
2249 Optimize_Single_Param_Generic(mixt_tree,&(tree->mod->ras->pinvar->v),.0001,0.9999,
2250 tree->mod->s_opt->min_diff_lk_local,
2251 tree->mod->s_opt->brent_it_max,
2252 tree->mod->s_opt->quickdirty);
2253
2254 tree->mod->s_opt->skip_tree_traversal = NO;
2255
2256 if(verbose == YES)
2257 {
2258 Print_Lk(mixt_tree,"[P-inv ]");
2259 PhyML_Printf("[%10f]",tree->mod->ras->pinvar->v);
2260 }
2261 }
2262 }
2263 }
2264 tree = tree->next_mixt;
2265 }
2266 while(tree);
2267
2268 if(alpha) Free(alpha);
2269
2270 return 1;
2271 }
2272
2273 //////////////////////////////////////////////////////////////
2274 //////////////////////////////////////////////////////////////
2275
Br_Len_Spline(phydbl * l,t_edge * b,int n_iter_max,phydbl tol,t_tree * tree)2276 static phydbl Br_Len_Spline(phydbl *l, t_edge *b, int n_iter_max, phydbl tol, t_tree *tree)
2277 {
2278 short int converged;
2279 phydbl init_lnL,old_lnL;
2280 int iter;
2281 phydbl best_l, new_l, init_l, init_dl, best_lnL;
2282 phydbl u, v;
2283 phydbl fu, fv;
2284 phydbl dfu, dfv;
2285 phydbl mult;
2286 phydbl a_,b_,A_,B_,C_,D_,root1,root2;
2287 short int ok1, ok2;
2288 // Warning: make sure eigen_lr vectors are already up-to-date
2289
2290 Set_Use_Eigen_Lr(YES,tree);
2291
2292
2293 best_l = init_l = *l;
2294 best_lnL = old_lnL = init_lnL = tree->c_lnL;
2295 mult = 1.2;
2296 ok1 = ok2 = NO;
2297 a_ = b_ = A_ = B_ = D_ = root1 = root2 = -1.;
2298 u = v = fu = fv = dfu = dfv = -1.;
2299 new_l = -1.;
2300
2301 dLk(l,b,tree);
2302 init_dl = tree->c_dlnL;
2303
2304 if(*l > tree->mod->l_max) *l = 0.5;
2305 if(*l < tree->mod->l_min) *l = 0.001;
2306
2307 // Find value of l where first derivative is < 0;
2308 tree->c_dlnL = init_dl;
2309 while(tree->c_dlnL < 0.0)
2310 {
2311 *l /= mult;
2312 tree->n_tot_bl_opt++;
2313 if(*l < tree->mod->l_min)
2314 {
2315 *l = best_l;
2316 tree->c_lnL = best_lnL;
2317 return best_lnL;
2318 }
2319 dLk(l,b,tree);
2320 if(tree->c_lnL > best_lnL)
2321 {
2322 best_lnL = tree->c_lnL;
2323 best_l = *l;
2324 }
2325 }
2326 u = *l;
2327 fu = tree->c_lnL;
2328 dfu = tree->c_dlnL;
2329
2330
2331
2332 // Find good upper bound
2333 *l = init_l;
2334 tree->c_dlnL = init_dl;
2335 tree->c_lnL = init_lnL;
2336
2337 while(tree->c_dlnL > 0.0)
2338 {
2339 *l *= mult;
2340 tree->n_tot_bl_opt++;
2341 if(*l > tree->mod->l_max)
2342 {
2343 *l = best_l;
2344 tree->c_lnL = best_lnL;
2345 return best_lnL;
2346 }
2347 dLk(l,b,tree);
2348 if(tree->c_lnL > best_lnL)
2349 {
2350 best_lnL = tree->c_lnL;
2351 best_l = *l;
2352 }
2353 }
2354 v = *l;
2355 fv = tree->c_lnL;
2356 dfv = tree->c_dlnL;
2357
2358
2359
2360 /* PhyML_Printf("\n Begin NR loop (lnL: %12G dlnL: %12G) l: %12G num: %d",tree->c_lnL,tree->c_dlnL,*l,b->num); */
2361
2362 converged = NO;
2363 iter = 0;
2364 do
2365 {
2366 /* PhyML_Printf("\n. l=%12f lnL=%12f iter:%d u=%12f v=%12f root1=%12f root2=%12f dfu=%12f dfv=%12f fu=%12f fv=%12f diff=%12f tol=%12f init=%12f", */
2367 /* *l,tree->c_lnL,iter, */
2368 /* u,v, */
2369 /* root1,root2, */
2370 /* dfu,dfv, */
2371 /* fu,fv, */
2372 /* tree->c_lnL-old_lnL, */
2373 /* tol,init_l); */
2374
2375 /* PhyML_Printf("\n. l: %f dlnL: %f lnL: %f",*l,tree->c_dlnL,tree->c_lnL); */
2376 // Spline interpolation (https://en.wikipedia.org/wiki/Spline_interpolation)
2377 a_ = dfu*(v-u) - (fv-fu);
2378 b_ = -dfv*(v-u) + (fv-fu);
2379
2380 A_ = 3.*a_ - 3.*b_;
2381 B_ = -4.*a_ + 2.*b_;
2382 C_ = fv-fu+a_;
2383 D_ = sqrt(B_*B_-4.*A_*C_);
2384
2385 root1 = (-B_-D_)/(2.*A_);
2386 root2 = (-B_+D_)/(2.*A_);
2387
2388 root1 = root1*(v-u) + u;
2389 root2 = root2*(v-u) + u;
2390
2391 ok1 = NO;
2392 ok2 = NO;
2393 if(root1 > u && root1 < v) ok1 = YES;
2394 if(root2 > u && root2 < v) ok2 = YES;
2395
2396 if(Are_Equal(root1,u,1.E-5) == YES) ok1 = YES;
2397 if(Are_Equal(root2,u,1.E-5) == YES) ok2 = YES;
2398 if(Are_Equal(root1,v,1.E-5) == YES) ok1 = YES;
2399 if(Are_Equal(root2,v,1.E-5) == YES) ok2 = YES;
2400
2401
2402 if(ok1 == YES && ok2 == YES) new_l = root1 < root2 ? root1 : root2;
2403 else if(ok1 == YES) new_l = root1;
2404 else if(ok2 == YES) new_l = root2;
2405 else if(u/v > 1.1 || u/v < 0.9)
2406 {
2407 PhyML_Printf("\n. iter=%4d u=%12G fu=%12G dfu=%12G v=%12G fv=%12G dfv=%12G root1=%12G root2=%12G\n",iter,u,fu,dfu,v,fv,dfv,root1,root2);
2408 assert(FALSE);
2409 }
2410
2411
2412 *l = new_l;
2413 tree->n_tot_bl_opt++;
2414
2415 old_lnL = tree->c_lnL;
2416 dLk(l,b,tree);
2417 if(tree->c_lnL > best_lnL)
2418 {
2419 best_lnL = tree->c_lnL;
2420 best_l = *l;
2421 }
2422
2423 if(tree->c_dlnL > 0.0)
2424 {
2425 u = new_l;
2426 fu = tree->c_lnL;
2427 dfu = tree->c_dlnL;
2428 }
2429 else
2430 {
2431 v = new_l;
2432 fv = tree->c_lnL;
2433 dfv = tree->c_dlnL;
2434 }
2435
2436
2437
2438 if(u - v < DBL_MIN) converged = YES;
2439 if(fabs(tree->c_lnL-old_lnL) < tol) converged = YES;
2440 if(++iter == n_iter_max+20) converged = YES;
2441 if(iter >= n_iter_max) PhyML_Fprintf(stderr,"\n. Edge length optimization took too long... l=%G lnL=%G iter:%d u=%G v=%G root1=%G root2=%G dfu=%G dfv=%G fu=%G fv=%G diff=%G tol=%G",
2442 *l,tree->c_lnL,iter,
2443 u,v,
2444 root1,root2,
2445 dfu,dfv,
2446 fu,fv,
2447 tree->c_lnL-old_lnL,
2448 tol);
2449
2450 if(converged == NO)
2451 {
2452 if(!(u < v)) PhyML_Printf("\n. u=%g v=%g.\n",u,v);
2453 if(!(dfu > 0.0)) PhyML_Printf("\n. dfu=%g l=%g u=%g v=%g\n",dfu,*l,u,v);
2454 if(!(dfv < 0.0)) PhyML_Printf("\n. dfv=%g l=%g u=%g v=%g\n",dfv,*l);
2455
2456 assert(u < v);
2457 assert(dfu > 0.0);
2458 assert(dfv < 0.0);
2459 }
2460 }
2461 while(converged == NO);
2462
2463 /* PhyML_Printf("\n. l = %g l_max = %g diff=%g",*l,tree->mod->l_max,*l-tree->mod->l_max); */
2464
2465 assert(!(*l > tree->mod->l_max));
2466 assert(!(*l < tree->mod->l_min));
2467
2468 /* if(dfu > 1.) */
2469 /* { */
2470 /* PhyML_Printf("\n> [%4d] l=%f lnL=%f iter:%d u=%f v=%f root1=%f root2=%f dfu=%f dfv=%f fu=%f fv=%f diff=%f tol=%f init=%f", */
2471 /* tree->n_tot_bl_opt, */
2472 /* *l,tree->c_lnL,iter, */
2473 /* u,v, */
2474 /* root1,root2, */
2475 /* dfu,dfv, */
2476 /* fu,fv, */
2477 /* tree->c_lnL-old_lnL, */
2478 /* tol,init_l); */
2479 /* } */
2480 /* if(dfv < -1) */
2481 /* { */
2482 /* PhyML_Printf("\n< [%4d] l=%f lnL=%f iter:%d u=%f v=%f root1=%f root2=%f dfu=%f dfv=%f fu=%f fv=%f diff=%f tol=%f init=%f", */
2483 /* tree->n_tot_bl_opt, */
2484 /* *l,tree->c_lnL,iter, */
2485 /* u,v, */
2486 /* root1,root2, */
2487 /* dfu,dfv, */
2488 /* fu,fv, */
2489 /* tree->c_lnL-old_lnL, */
2490 /* tol,init_l); */
2491 /* } */
2492
2493 *l = best_l;
2494 tree->c_lnL = best_lnL;
2495
2496 if(iter == n_iter_max)
2497 {
2498 PhyML_Printf("\n. Too many iterations in edge length optimization routine (l=%G init=%G).\n",best_l,init_l);
2499 assert(FALSE);
2500 }
2501
2502 return tree->c_lnL;
2503 }
2504
2505 //////////////////////////////////////////////////////////////
2506 //////////////////////////////////////////////////////////////
2507
Generic_Brent_Lk(phydbl * param,phydbl ax,phydbl cx,phydbl tol,int n_iter_max,int quickdirty,phydbl (* obj_func)(t_edge *,t_tree *,supert_tree *),t_edge * branch,t_tree * tree,supert_tree * stree,int logt)2508 int Generic_Brent_Lk(phydbl *param, phydbl ax, phydbl cx, phydbl tol,
2509 int n_iter_max, int quickdirty,
2510 phydbl (*obj_func)(t_edge *,t_tree *,supert_tree *),
2511 t_edge *branch, t_tree *tree, supert_tree *stree, int logt)
2512 {
2513 int iter;
2514 phydbl a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
2515 phydbl e=0.0;
2516 phydbl old_lnL,init_lnL;
2517 phydbl bx = *param;
2518 int n_opt_step;
2519
2520 n_opt_step = 0;
2521 d=0.0;
2522 a=((ax < cx) ? ax : cx);
2523 b=((ax > cx) ? ax : cx);
2524 x=w=v=bx;
2525 (*param) = bx;
2526 if(logt == YES) (*param) = exp(MIN(1.E+2,*param));
2527 fw=fv=fx=fu=-(*obj_func)(branch,tree,stree);
2528 if(logt == YES) (*param) = log(*param);
2529 init_lnL = old_lnL = fw;
2530
2531 /* PhyML_Printf("\n. %p %p %p init_lnL=%f fu=%f ax=%f cx=%f param=%f",branch,tree,stree,init_lnL,fu,ax,cx,*param); */
2532
2533 for(iter=1;iter<=BRENT_IT_MAX;iter++)
2534 {
2535 xm=0.5*(a+b);
2536 tol2=2.0*(tol1=tol*x+BRENT_ZEPS);
2537
2538 if((fu < init_lnL + tol) && (quickdirty == YES) && (iter > 1))
2539 {
2540 (*param) = x;
2541 if(logt == YES) (*param) = exp(MIN(1.E+2,*param));
2542 fu = (*obj_func)(branch,tree,stree);
2543 if(logt == YES) (*param) = log(*param);
2544 /* printf("\n. return %f [%f] %d",fu,*param,iter); */
2545 return n_opt_step;
2546 }
2547
2548 if((FABS(fu-old_lnL) < tol && iter > 1) || (iter > n_iter_max - 1))
2549 {
2550 (*param) = x;
2551 if(logt == YES) (*param) = exp(MIN(1.E+2,*param));
2552 fu = (*obj_func)(branch,tree,stree);
2553 if(logt == YES) (*param) = log(*param);
2554 /* printf("\n. return %f [%f] %d",*param,fu,iter); */
2555 return n_opt_step;
2556 }
2557
2558 if(FABS(e) > tol1)
2559 {
2560 r=(x-w)*(fx-fv);
2561 q=(x-v)*(fx-fw);
2562 p=(x-v)*q-(x-w)*r;
2563 q=2.0*(q-r);
2564 if(q > 0.0) p = -p;
2565 q=FABS(q);
2566 etemp=e;
2567 e=d;
2568 if(FABS(p) >= FABS(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
2569 {
2570 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
2571 /* PhyML_Printf("\n. Golden section step"); */
2572 }
2573 else
2574 {
2575 d=p/q;
2576 u=x+d;
2577 if (u-a < tol2 || b-u < tol2) d=SIGN(tol1,xm-x);
2578 /* PhyML_Printf("\n. Parabolic step [e=%f]",e); */
2579 }
2580 }
2581 else
2582 {
2583 d=BRENT_CGOLD*(e=(x >= xm ? a-x : b-x));
2584 /* PhyML_Printf("\n. Golden section step (default) [e=%f tol1=%f a=%f b=%f d=%f x=%f]",e,tol1,a,b,d,x); */
2585 }
2586
2587 u=(FABS(d) >= tol1 ? x+d : x+SIGN(tol1,d));
2588 (*param) = u;
2589 n_opt_step++;
2590 old_lnL = fu;
2591 if(logt == YES) (*param) = exp(MIN(1.E+2,*param));
2592 fu = -(*obj_func)(branch,tree,stree);
2593 if(logt == YES) (*param) = log(*param);
2594 /* PhyML_Printf("\n. iter=%d/%d param=%f lnL=%f u: %f x: %f d: %f logt: %d",iter,BRENT_IT_MAX,*param,fu,u,x,d,logt); */
2595
2596 if(fu <= fx)
2597 {
2598 if(u >= x) a=x; else b=x;
2599 SHFT(v,w,x,u)
2600 SHFT(fv,fw,fx,fu)
2601 }
2602 else
2603 {
2604 if (u < x) a=u; else b=u;
2605 if (fu < fw || FABS(w-x) < SMALL)
2606 {
2607 v=w;
2608 w=u;
2609 fv=fw;
2610 fw=fu;
2611 }
2612 /* else if (fu <= fv || v == x || v == w) */
2613 else if (fu < fv || FABS(v-x) < SMALL || FABS(v-w) < SMALL)
2614 {
2615 v=u;
2616 fv=fu;
2617 }
2618 }
2619 }
2620
2621 PhyML_Printf("\n. Too many iterations in Generic_Brent_Lk !");
2622 assert(FALSE);
2623 return(n_opt_step);
2624 /* Not Reached ?? *param=x; */
2625 /* Not Reached ?? return fx; */
2626 }
2627
2628 //////////////////////////////////////////////////////////////
2629 //////////////////////////////////////////////////////////////
2630
2631
2632 /* find ML erstimates of node heights given fixed substitution
2633 rates on branches. Also optimizes the overall substitution
2634 rate */
Round_Optimize_Node_Heights(t_tree * tree)2635 void Round_Optimize_Node_Heights(t_tree *tree)
2636 {
2637 phydbl cur_lnL, new_lnL;
2638 int n_iter;
2639
2640
2641 cur_lnL = UNLIKELY;
2642 new_lnL = Lk(NULL,tree);
2643
2644
2645 n_iter = 0;
2646 while(fabs(new_lnL - cur_lnL) > tree->mod->s_opt->min_diff_lk_local)
2647 {
2648 cur_lnL = tree->c_lnL;
2649
2650 Opt_Node_Heights_Recurr(tree);
2651
2652 Generic_Brent_Lk(&(tree->rates->clock_r),
2653 tree->rates->min_clock,
2654 tree->rates->max_clock,
2655 tree->mod->s_opt->min_diff_lk_local,
2656 tree->mod->s_opt->brent_it_max,
2657 tree->mod->s_opt->quickdirty,
2658 Wrap_Lk,NULL,tree,NULL,NO);
2659
2660 printf("\n. cur_lnL=%f new_lnL=%f clock_r=%G root height=%f",
2661 cur_lnL,new_lnL,tree->rates->clock_r,tree->times->nd_t[tree->n_root->num]);
2662 new_lnL = tree->c_lnL;
2663 n_iter++;
2664 if(n_iter > 100) break;
2665 }
2666 }
2667
2668 //////////////////////////////////////////////////////////////
2669 //////////////////////////////////////////////////////////////
2670
2671
Opt_Node_Heights_Recurr(t_tree * tree)2672 void Opt_Node_Heights_Recurr(t_tree *tree)
2673 {
2674 Opt_Node_Heights_Recurr_Pre(tree->n_root,tree->n_root->v[2],tree);
2675 Opt_Node_Heights_Recurr_Pre(tree->n_root,tree->n_root->v[1],tree);
2676
2677 Generic_Brent_Lk(&(tree->times->nd_t[tree->n_root->num]),
2678 MIN(tree->times->t_prior_max[tree->n_root->num],
2679 MIN(tree->times->nd_t[tree->n_root->v[2]->num],
2680 tree->times->nd_t[tree->n_root->v[1]->num])),
2681 tree->times->t_prior_min[tree->n_root->num],
2682 tree->mod->s_opt->min_diff_lk_local,
2683 tree->mod->s_opt->brent_it_max,
2684 tree->mod->s_opt->quickdirty,
2685 Wrap_Lk,NULL,tree,NULL,NO);
2686 }
2687
2688 //////////////////////////////////////////////////////////////
2689 //////////////////////////////////////////////////////////////
2690
2691
Opt_Node_Heights_Recurr_Pre(t_node * a,t_node * d,t_tree * tree)2692 void Opt_Node_Heights_Recurr_Pre(t_node *a, t_node *d, t_tree *tree)
2693 {
2694 if(d->tax) return;
2695 else
2696 {
2697 int i;
2698 phydbl t0,t2,t3;
2699 phydbl t_min,t_max;
2700 t_node *v2,*v3;
2701
2702 v2 = v3 = NULL;
2703 for(i=0;i<3;i++)
2704 if((d->v[i] != a) && (d->b[i] != tree->e_root))
2705 {
2706 if(!v2) { v2 = d->v[i]; }
2707 else { v3 = d->v[i]; }
2708 }
2709
2710 Opt_Node_Heights_Recurr_Pre(d,v2,tree);
2711 Opt_Node_Heights_Recurr_Pre(d,v3,tree);
2712
2713 t0 = tree->times->nd_t[a->num];
2714 t2 = tree->times->nd_t[v2->num];
2715 t3 = tree->times->nd_t[v3->num];
2716
2717 t_min = t0;
2718 t_max = MIN(t2,t3);
2719
2720 t_min = MAX(t_min,tree->times->t_prior_min[d->num]);
2721 t_max = MIN(t_max,tree->times->t_prior_max[d->num]);
2722
2723 t_min += tree->rates->min_dt;
2724 t_max -= tree->rates->min_dt;
2725
2726 if(t_min > t_max)
2727 {
2728 PhyML_Fprintf(stderr,"\n. Err in file %s at line %d\n",__FILE__,__LINE__);
2729 Exit("\n");
2730 }
2731
2732 Generic_Brent_Lk(&(tree->times->nd_t[d->num]),
2733 t_min,t_max,
2734 tree->mod->s_opt->min_diff_lk_local,
2735 tree->mod->s_opt->brent_it_max,
2736 tree->mod->s_opt->quickdirty,
2737 Wrap_Lk,NULL,tree,NULL,NO);
2738
2739 /* printf("\n. t%d = %f [%f;%f] lnL = %f",d->num,tree->times->nd_t[d->num],t_min,t_max,tree->c_lnL); */
2740
2741 }
2742 }
2743
2744
2745 //////////////////////////////////////////////////////////////
2746 //////////////////////////////////////////////////////////////
2747
Optimize_RR_Params(t_tree * mixt_tree,int verbose)2748 void Optimize_RR_Params(t_tree *mixt_tree, int verbose)
2749 {
2750 t_tree *tree;
2751 t_rmat **r_mat;
2752 int *permut;
2753 int n_r_mat;
2754 int i;
2755 phydbl lk_new,lk_old;
2756 phydbl a,b;
2757 phydbl *opt_val;
2758
2759 Set_Update_Eigen(YES,mixt_tree->mod);
2760
2761 n_r_mat = 0;
2762 tree = mixt_tree;
2763 r_mat = NULL;
2764 permut = NULL;
2765 lk_old = UNLIKELY;
2766 lk_new = UNLIKELY;
2767
2768 do
2769 {
2770 if(tree->is_mixt_tree == YES) tree = tree->next;
2771
2772 for(i=0;i<n_r_mat;i++) if(tree->mod->r_mat == r_mat[i]) break;
2773
2774 if(i == n_r_mat) // tree->mod->r_mat was not found before
2775 {
2776 if(!r_mat) r_mat = (t_rmat **)mCalloc(1,sizeof(t_rmat *));
2777 else r_mat = (t_rmat **)mRealloc(r_mat,n_r_mat+1,sizeof(t_rmat *));
2778 r_mat[n_r_mat] = tree->mod->r_mat;
2779 n_r_mat++;
2780
2781 if((tree->mod->s_opt->opt_rr) &&
2782 ((tree->mod->whichmodel == GTR) ||
2783 ((tree->mod->whichmodel == CUSTOM) &&
2784 (tree->mod->r_mat->n_diff_rr > 1))))
2785 {
2786 int i,iter;
2787
2788 opt_val = (phydbl *)mCalloc(tree->mod->r_mat->n_diff_rr,sizeof(phydbl));
2789
2790 iter = 0;
2791 do
2792 {
2793 lk_old = mixt_tree->c_lnL;
2794 int failed = NO;
2795
2796 if(tree->mod->r_mat->n_diff_rr > 2)
2797 {
2798 for(i=0;i<tree->mod->r_mat->n_diff_rr;i++) opt_val[i] = tree->mod->r_mat->rr_val->v[i];
2799
2800 BFGS(mixt_tree,tree->mod->r_mat->rr_val->v,tree->mod->r_mat->n_diff_rr,1.e-5,tree->mod->s_opt->min_diff_lk_local,1.e-3,NO,NO,
2801 &Return_Abs_Lk,
2802 &Num_Derivative_Several_Param,
2803 &Lnsrch,&failed);
2804
2805
2806 if(failed == YES) for(i=0;i<tree->mod->r_mat->n_diff_rr;i++) tree->mod->r_mat->rr_val->v[i] = opt_val[i];
2807 }
2808
2809 permut = Permutate(tree->mod->r_mat->n_diff_rr);
2810
2811 for(i=0;i<tree->mod->r_mat->n_diff_rr;i++) opt_val[i] = tree->mod->r_mat->rr_val->v[i];
2812
2813 /* for(i=0;i<tree->mod->r_mat->n_diff_rr;i++) tree->mod->r_mat->rr_val->v[i] = 0.0; */
2814
2815
2816 for(i=0;i<tree->mod->r_mat->n_diff_rr;i++)
2817 {
2818 a = tree->mod->r_mat->rr_val->v[permut[i]]/2.;
2819 b = tree->mod->r_mat->rr_val->v[permut[i]]*2.+1.;
2820
2821 Generic_Brent_Lk(&(tree->mod->r_mat->rr_val->v[permut[i]]),
2822 a,b,
2823 tree->mod->s_opt->min_diff_lk_local,
2824 tree->mod->s_opt->brent_it_max,
2825 tree->mod->s_opt->quickdirty,
2826 Wrap_Lk,NULL,mixt_tree,NULL,NO);
2827 }
2828
2829
2830 if(mixt_tree->c_lnL < lk_old)
2831 {
2832 for(i=0;i<tree->mod->r_mat->n_diff_rr;i++) tree->mod->r_mat->rr_val->v[i] = opt_val[i];
2833 Lk(NULL,mixt_tree);
2834 }
2835
2836
2837 if(verbose) Print_Lk(tree->mixt_tree?
2838 tree->mixt_tree:
2839 tree,"[GTR parameters ]");
2840
2841 lk_new = mixt_tree->c_lnL;
2842
2843 Free(permut);
2844
2845
2846 if(lk_new < lk_old - tree->mod->s_opt->min_diff_lk_global)
2847 {
2848 PhyML_Printf("\n. lk_new: %f lk_old: %f",lk_new,lk_old);
2849 assert(FALSE);
2850 }
2851
2852 if(fabs(lk_new-lk_old) < tree->mod->s_opt->min_diff_lk_local) break;
2853 }
2854 /* while(++iter < tree->mod->s_opt->brent_it_max); */
2855 while(++iter < 1);
2856
2857 Free(opt_val);
2858
2859 if(iter == tree->mod->s_opt->brent_it_max)
2860 {
2861 if(tree->verbose > VL0)
2862 {
2863 PhyML_Printf("\n. Failed to optimize GTR parameters this round...");
2864 }
2865 }
2866 }
2867 }
2868
2869 tree = tree->next;
2870 if(!tree) break;
2871
2872 }
2873 while(1);
2874
2875 if(r_mat) Free(r_mat);
2876
2877
2878 Set_Update_Eigen(NO,mixt_tree->mod);
2879
2880 }
2881
2882 //////////////////////////////////////////////////////////////
2883 //////////////////////////////////////////////////////////////
2884
Optimize_TsTv(t_tree * mixt_tree,int verbose)2885 void Optimize_TsTv(t_tree *mixt_tree, int verbose)
2886 {
2887 scalar_dbl **tstv;
2888 int n_tstv;
2889 t_tree *tree;
2890 int i;
2891
2892 Set_Update_Eigen(YES,mixt_tree->mod);
2893
2894 tstv = NULL;
2895 n_tstv = 0;
2896 tree = mixt_tree;
2897
2898 do
2899 {
2900 if(tree->is_mixt_tree == YES) tree = tree->next;
2901
2902 for(i=0;i<n_tstv;i++) if(tree->mod->kappa == tstv[i]) break;
2903
2904 if(i == n_tstv)
2905 {
2906 if(!tstv) tstv = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
2907 else tstv = (scalar_dbl **)mRealloc(tstv,n_tstv+1,sizeof(scalar_dbl *));
2908 tstv[n_tstv] = tree->mod->kappa;
2909 n_tstv++;
2910
2911 if(tree->mod->s_opt->opt_kappa == YES)
2912 {
2913 phydbl a,c;
2914
2915 /* a = tree->mod->kappa->v * .1; */
2916 /* c = tree->mod->kappa->v * 10.; */
2917 a = TSTV_MIN;
2918 c = TSTV_MAX;
2919
2920 Generic_Brent_Lk(&(tree->mod->kappa->v),
2921 a,c,
2922 tree->mod->s_opt->min_diff_lk_local,
2923 tree->mod->s_opt->brent_it_max,
2924 tree->mod->s_opt->quickdirty,
2925 Wrap_Lk,NULL,mixt_tree,NULL,NO);
2926
2927 if(verbose)
2928 {
2929 Print_Lk(mixt_tree,"[Ts/ts ratio ]");
2930 PhyML_Printf("[%10f]",tree->mod->kappa->v);
2931 }
2932 }
2933 }
2934
2935 tree = tree->next;
2936
2937 }
2938 while(tree);
2939
2940 if(tstv) Free(tstv);
2941
2942 Set_Update_Eigen(NO,mixt_tree->mod);
2943
2944 }
2945
2946 //////////////////////////////////////////////////////////////
2947 //////////////////////////////////////////////////////////////
2948
Optimize_Pinv(t_tree * mixt_tree,int verbose)2949 void Optimize_Pinv(t_tree *mixt_tree, int verbose)
2950 {
2951 scalar_dbl **pinv;
2952 int n_pinv;
2953 t_tree *tree;
2954 int i;
2955
2956 Set_Update_Eigen(NO,mixt_tree->mod);
2957
2958 pinv = NULL;
2959 n_pinv = 0;
2960 tree = mixt_tree;
2961
2962 do
2963 {
2964 for(i=0;i<n_pinv;i++) if(tree->mod->ras->pinvar == pinv[i]) break;
2965
2966 if(i == n_pinv)
2967 {
2968 if(!pinv) pinv = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
2969 else pinv = (scalar_dbl **)mRealloc(pinv,n_pinv+1,sizeof(scalar_dbl *));
2970 pinv[n_pinv] = tree->mod->ras->pinvar;
2971 n_pinv++;
2972
2973 if(tree->mod->s_opt->opt_pinvar == YES && (tree->mod->s_opt->opt_alpha == NO || tree->mod->ras->n_catg == 1))
2974 {
2975 Generic_Brent_Lk(&(tree->mod->ras->pinvar->v),
2976 PINV_MIN,PINV_MAX,
2977 tree->mod->s_opt->min_diff_lk_local,
2978 tree->mod->s_opt->brent_it_max,
2979 tree->mod->s_opt->quickdirty,
2980 Wrap_Lk,NULL,mixt_tree,NULL,NO);
2981
2982 if(verbose)
2983 {
2984 Print_Lk(tree,"[P-inv ]");
2985 PhyML_Printf("[%10f]",tree->mod->ras->pinvar->v);
2986 }
2987 }
2988 }
2989
2990 tree = tree->next_mixt;
2991
2992 }
2993 while(tree);
2994
2995 if(pinv) Free(pinv);
2996
2997 }
2998
2999 //////////////////////////////////////////////////////////////
3000 //////////////////////////////////////////////////////////////
3001
Optimize_Alpha(t_tree * mixt_tree,int verbose)3002 void Optimize_Alpha(t_tree *mixt_tree, int verbose)
3003 {
3004 scalar_dbl **alpha;
3005 int n_alpha;
3006 t_tree *tree;
3007 int i;
3008
3009 Set_Update_Eigen(NO,mixt_tree->mod);
3010
3011 alpha = NULL;
3012 n_alpha = 0;
3013 tree = mixt_tree;
3014
3015 do
3016 {
3017
3018 if(tree->mod->s_opt->opt_alpha == YES && tree->mod->ras->n_catg > 1)
3019 {
3020 for(i=0;i<n_alpha;i++) if(tree->mod->ras->alpha == alpha[i]) break;
3021
3022 if(i == n_alpha)
3023 {
3024 if(!alpha) alpha = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
3025 else alpha = (scalar_dbl **)mRealloc(alpha,n_alpha+1,sizeof(scalar_dbl *));
3026 alpha[n_alpha] = tree->mod->ras->alpha;
3027 n_alpha++;
3028
3029 if(tree->mod->s_opt->opt_alpha == YES &&
3030 tree->mod->ras->free_mixt_rates == NO &&
3031 tree->mod->s_opt->opt_pinvar == NO)
3032 {
3033
3034 if(tree->mod->ras->n_catg > 1)
3035 {
3036 Generic_Brent_Lk(&(tree->mod->ras->alpha->v),
3037 ALPHA_MIN,ALPHA_MAX,
3038 tree->mod->s_opt->min_diff_lk_local,
3039 tree->mod->s_opt->brent_it_max,
3040 tree->mod->s_opt->quickdirty,
3041 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3042 }
3043 if(verbose)
3044 {
3045 Print_Lk(mixt_tree,"[Alpha ]");
3046 PhyML_Printf("[%10f]",tree->mod->ras->alpha->v);
3047 }
3048 }
3049 }
3050 }
3051 tree = tree->next_mixt;
3052 }
3053 while(tree);
3054
3055 if(alpha) Free(alpha);
3056
3057 }
3058
3059 //////////////////////////////////////////////////////////////
3060 //////////////////////////////////////////////////////////////
3061
Optimize_Free_Rate(t_tree * mixt_tree,int verbose)3062 void Optimize_Free_Rate(t_tree *mixt_tree, int verbose)
3063 {
3064 t_tree *tree;
3065 int fast;
3066 int i,pos,failed;
3067 int *permut;
3068 phydbl lk_before, lk_after;
3069 tree = mixt_tree;
3070
3071 lk_before = lk_after = UNLIKELY;
3072
3073 do
3074 {
3075 if((tree->mod->s_opt->opt_free_mixt_rates) && (tree->mod->ras->free_mixt_rates == YES) && (tree->mod->ras->n_catg > 1))
3076 {
3077 if(tree->mod->s_opt->serial_free_rates == YES)
3078 {
3079 fast = YES;
3080 lk_before = tree->c_lnL;
3081 Optimize_Free_Rate_Weights(tree,fast,verbose);
3082 lk_after = tree->c_lnL;
3083 Optimize_Free_Rate_Rr(tree,fast,verbose);
3084 lk_after = tree->c_lnL;
3085
3086 if(lk_after < lk_before - tree->mod->s_opt->min_diff_lk_global)
3087 {
3088 PhyML_Fprintf(stderr,"\n. lk_before: %f lk_after: %f diff: %G",lk_before,lk_after,lk_before-lk_after);
3089 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3090 Exit("");
3091 }
3092 }
3093
3094 if(FABS(lk_before - lk_after) > 0.001)
3095 {
3096 phydbl **x,*cpy;
3097
3098 x = (phydbl **)mCalloc(2*tree->mod->ras->n_catg,sizeof(phydbl *));
3099 cpy = (phydbl *)mCalloc(2*tree->mod->ras->n_catg,sizeof(phydbl));
3100
3101
3102 lk_before = tree->c_lnL;
3103
3104 pos = 0;
3105 for(i=0;i<tree->mod->ras->n_catg;i++) x[pos++] = tree->mod->ras->gamma_rr_unscaled->v+i;
3106 for(i=0;i<tree->mod->ras->n_catg;i++) x[pos++] = tree->mod->ras->gamma_r_proba_unscaled->v+i;
3107
3108 pos = 0;
3109 for(i=0;i<tree->mod->ras->n_catg;i++) cpy[pos++] = tree->mod->ras->gamma_rr_unscaled->v[i];
3110 for(i=0;i<tree->mod->ras->n_catg;i++) cpy[pos++] = tree->mod->ras->gamma_r_proba_unscaled->v[i];
3111
3112
3113 for(i=0;i<2*tree->mod->ras->n_catg;++i) *(x[i]) = log(MAX(1.E-10,*(x[i])));
3114
3115 failed = YES;
3116 BFGS_Nonaligned(tree,x,2*tree->mod->ras->n_catg,1.e-5,tree->mod->s_opt->min_diff_lk_global,1.e-5,YES,NO,
3117 &Return_Abs_Lk,
3118 &Num_Derivative_Several_Param_Nonaligned,
3119 &Lnsrch_Nonaligned,&failed);
3120
3121
3122 for(i=0;i<2*tree->mod->ras->n_catg;++i) *(x[i]) = exp(MIN(1.E+2,*(x[i])));
3123
3124
3125 if(failed == YES)
3126 {
3127 pos = 0;
3128 for(i=0;i<tree->mod->ras->n_catg;i++) tree->mod->ras->gamma_rr_unscaled->v[i] = cpy[pos++];
3129 for(i=0;i<tree->mod->ras->n_catg;i++) tree->mod->ras->gamma_r_proba_unscaled->v[i] = cpy[pos++];
3130
3131 permut = Permutate(tree->mod->ras->n_catg);
3132
3133
3134 for(i=0;i<tree->mod->ras->n_catg;++i)
3135 {
3136 phydbl a,c;
3137
3138 a = tree->mod->ras->gamma_rr_unscaled->v[permut[i]] / 2.;
3139 c = tree->mod->ras->gamma_rr_unscaled->v[permut[i]] * 2.+1;
3140
3141 Generic_Brent_Lk(&(tree->mod->ras->gamma_rr_unscaled->v[permut[i]]),
3142 a,c,
3143 tree->mod->s_opt->min_diff_lk_local,
3144 tree->mod->s_opt->brent_it_max,
3145 tree->mod->s_opt->quickdirty,
3146 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3147 /* PhyML_Printf("\n. ++ lk=%f",mixt_tree->c_lnL); */
3148 }
3149
3150 for(i=0;i<tree->mod->ras->n_catg;++i)
3151 {
3152 phydbl a,c;
3153
3154 a = tree->mod->ras->gamma_r_proba_unscaled->v[permut[i]] / 2.;
3155 c = tree->mod->ras->gamma_r_proba_unscaled->v[permut[i]] * 2.;
3156
3157 Generic_Brent_Lk(&(tree->mod->ras->gamma_r_proba_unscaled->v[permut[i]]),
3158 a,c,
3159 tree->mod->s_opt->min_diff_lk_local,
3160 tree->mod->s_opt->brent_it_max,
3161 tree->mod->s_opt->quickdirty,
3162 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3163 /* PhyML_Printf("\n. ++ lk=%f",mixt_tree->c_lnL); */
3164 }
3165
3166 Free(permut);
3167
3168 }
3169
3170
3171 lk_after = tree->c_lnL;
3172
3173 if(lk_after < lk_before - tree->mod->s_opt->min_diff_lk_global)
3174 {
3175 PhyML_Fprintf(stderr,"\n. lk_before: %f lk_after: %f diff: %G",lk_before,lk_after,lk_before-lk_after);
3176 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3177 Exit("");
3178 }
3179
3180 Free(x);
3181 Free(cpy);
3182 }
3183
3184
3185 }
3186 tree = tree->next_mixt;
3187 }
3188 while(tree);
3189 }
3190
3191 //////////////////////////////////////////////////////////////
3192 //////////////////////////////////////////////////////////////
3193
Optimize_Free_Rate_Rr(t_tree * tree,int fast,int verbose)3194 void Optimize_Free_Rate_Rr(t_tree *tree, int fast, int verbose)
3195 {
3196 phydbl lk_before, lk_after;
3197
3198 lk_before = tree->c_lnL;
3199
3200
3201 if(tree->prev == NULL && tree->next == NULL)
3202 {
3203 int i;
3204 phydbl wm;
3205
3206 /* tree->mod->s_opt->curr_opt_free_rates = YES; */
3207
3208 if(fast == YES)
3209 {
3210 for(i=0;i<tree->mod->ras->n_catg;i++) tree->mod->ras->skip_rate_cat[i] = YES;
3211 tree->mod->ras->normalise_rr = NO;
3212
3213 wm = Weighted_Mean(tree->mod->ras->gamma_rr_unscaled->v,
3214 tree->mod->ras->gamma_r_proba->v,
3215 tree->mod->ras->n_catg);
3216
3217 tree->mod->ras->free_rate_mr->v = 100.;
3218 For(i,2*tree->n_otu-1) tree->a_edges[i]->l->v /= (wm * tree->mod->ras->free_rate_mr->v);
3219 }
3220
3221
3222 for(i=0;i<tree->mod->ras->n_catg-1;i++)
3223 {
3224 if(fast == YES) tree->mod->ras->skip_rate_cat[i] = NO;
3225
3226 phydbl a,c;
3227
3228 a = tree->mod->ras->gamma_rr_unscaled->v[i] * .1;
3229 c = tree->mod->ras->gamma_rr_unscaled->v[i] * 10.+1;
3230
3231 Generic_Brent_Lk(&(tree->mod->ras->gamma_rr_unscaled->v[i]),
3232 a,c,
3233 tree->mod->s_opt->min_diff_lk_local,
3234 tree->mod->s_opt->brent_it_max,
3235 tree->mod->s_opt->quickdirty,
3236 Wrap_Lk,NULL,tree,NULL,NO);
3237
3238 if(fast == YES) tree->mod->ras->skip_rate_cat[i] = YES;
3239
3240 lk_after = tree->c_lnL;
3241
3242 if(lk_after < lk_before - tree->mod->s_opt->min_diff_lk_global)
3243 {
3244 PhyML_Fprintf(stderr,"\n. lk_before: %f lk_after: %f diff: %G",lk_before,lk_after,lk_before-lk_after);
3245 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3246 Exit("");
3247 }
3248 }
3249
3250 if(fast == YES)
3251 {
3252 for(i=0;i<tree->mod->ras->n_catg;i++) tree->mod->ras->skip_rate_cat[i] = NO;
3253 tree->mod->ras->normalise_rr = YES;
3254
3255 wm = Weighted_Mean(tree->mod->ras->gamma_rr_unscaled->v,
3256 tree->mod->ras->gamma_r_proba->v,
3257 tree->mod->ras->n_catg);
3258
3259 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->l->v *= (wm * tree->mod->ras->free_rate_mr->v);
3260 }
3261
3262 /* tree->mod->s_opt->curr_opt_free_rates = NO; */
3263
3264 }
3265 else
3266 {
3267 int i;
3268 for(i=0;i<tree->mod->ras->n_catg-1;i++)
3269 {
3270 phydbl a,c;
3271
3272 a = tree->mod->ras->gamma_rr_unscaled->v[i] * .1;
3273 c = tree->mod->ras->gamma_rr_unscaled->v[i] * 10.+1.;
3274
3275 Generic_Brent_Lk(&(tree->mod->ras->gamma_rr_unscaled->v[i]),
3276 a,c,
3277 tree->mod->s_opt->min_diff_lk_local,
3278 tree->mod->s_opt->brent_it_max,
3279 tree->mod->s_opt->quickdirty,
3280 Wrap_Lk,NULL,tree,NULL,NO);
3281 }
3282 }
3283
3284 lk_after = tree->c_lnL;
3285
3286 if(lk_after < lk_before - tree->mod->s_opt->min_diff_lk_global)
3287 {
3288 PhyML_Fprintf(stderr,"\n. lk_before: %f lk_after: %f diff: %G",lk_before,lk_after,lk_before-lk_after);
3289 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3290 Exit("");
3291 }
3292
3293 if(verbose) Print_Lk(tree,"[Rate class values ]");
3294
3295 /* int i; */
3296 /* for(i=0;i<tree->mod->ras->n_catg;i++) */
3297 /* { */
3298 /* printf("\n+ c %2d p: %15f r: %15f up: %15f ur: %5f", */
3299 /* i+1, */
3300 /* tree->mod->ras->gamma_r_proba->v[i], */
3301 /* tree->mod->ras->gamma_rr->v[i], */
3302 /* tree->mod->ras->gamma_r_proba_unscaled->v[i], */
3303 /* tree->mod->ras->gamma_rr_unscaled->v[i]); */
3304 /* } */
3305 /* fflush(NULL); */
3306
3307 /* printf("\n. LK: %f",Lk(NULL,tree)); */
3308
3309 /* int i; */
3310 /* printf("\n"); */
3311 /* printf("X*X %f ",tree->c_lnL); */
3312 /* /\* for(i=0;i<tree->mod->ras->n_catg;i++) printf("%f ",tree->mod->ras->gamma_rr_unscaled->v[i]); *\/ */
3313 /* for(i=0;i<tree->mod->ras->n_catg;i++) printf("%f ",tree->mod->ras->gamma_rr->v[i]); */
3314 /* for(i=0;i<tree->mod->ras->n_catg;i++) printf("%f ",tree->mod->ras->gamma_r_proba->v[i]); */
3315 /* For(i,2*tree->n_otu-3) printf("%f ",tree->a_edges[i]->l->v); */
3316 }
3317
3318 //////////////////////////////////////////////////////////////
3319 //////////////////////////////////////////////////////////////
3320
Optimize_Free_Rate_Weights(t_tree * tree,int fast,int verbose)3321 void Optimize_Free_Rate_Weights(t_tree *tree, int fast, int verbose)
3322 {
3323 int i;
3324 phydbl wm;
3325 phydbl lk_before, lk_after;
3326
3327
3328 lk_before = tree->c_lnL;
3329
3330 /*! Only skip tree traversal when data is not partitionned */
3331 if(tree->prev == NULL && tree->next == NULL && fast == YES)
3332 {
3333 tree->mod->s_opt->skip_tree_traversal = YES;
3334 tree->mod->ras->normalise_rr = NO;
3335
3336 wm = Weighted_Mean(tree->mod->ras->gamma_rr_unscaled->v,
3337 tree->mod->ras->gamma_r_proba->v,
3338 tree->mod->ras->n_catg);
3339
3340 tree->mod->ras->free_rate_mr->v = 100.;
3341 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->l->v /= (wm * tree->mod->ras->free_rate_mr->v);
3342 }
3343
3344 for(i=0;i<tree->mod->ras->n_catg-1;i++)
3345 {
3346 phydbl a,c;
3347
3348 a = tree->mod->ras->gamma_r_proba_unscaled->v[i] * .1;
3349 c = tree->mod->ras->gamma_r_proba_unscaled->v[i] * 10.+1;
3350
3351 Generic_Brent_Lk(&(tree->mod->ras->gamma_r_proba_unscaled->v[i]),
3352 a,c,
3353 tree->mod->s_opt->min_diff_lk_local,
3354 tree->mod->s_opt->brent_it_max,
3355 tree->mod->s_opt->quickdirty,
3356 Wrap_Lk,NULL,tree,NULL,NO);
3357 }
3358
3359 if(tree->mod->s_opt->skip_tree_traversal == YES && fast == YES)
3360 {
3361 tree->mod->s_opt->skip_tree_traversal = NO;
3362 tree->mod->ras->normalise_rr = YES;
3363
3364 wm = Weighted_Mean(tree->mod->ras->gamma_rr_unscaled->v,
3365 tree->mod->ras->gamma_r_proba->v,
3366 tree->mod->ras->n_catg);
3367
3368 for(i=0;i<2*tree->n_otu-1;++i) tree->a_edges[i]->l->v *= (wm * tree->mod->ras->free_rate_mr->v);
3369 }
3370
3371 lk_after = tree->c_lnL;
3372
3373 if(lk_after < lk_before - tree->mod->s_opt->min_diff_lk_global)
3374 {
3375 PhyML_Fprintf(stderr,"\n. lk_before: %f lk_after: %f diff: %G",lk_before,lk_after,lk_before-lk_after);
3376 PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
3377 Exit("");
3378 }
3379
3380 if(verbose) Print_Lk(tree,"[Rate class freqs. ]");
3381 }
3382
3383 //////////////////////////////////////////////////////////////
3384 //////////////////////////////////////////////////////////////
3385
Optimize_State_Freqs(t_tree * mixt_tree,int verbose)3386 void Optimize_State_Freqs(t_tree *mixt_tree, int verbose)
3387 {
3388 /* vect_dbl **freqs; */
3389 /* int n_freqs; */
3390 t_tree *tree;
3391 int i;
3392 phydbl lk_new,lk_old;
3393 int *permut;
3394 phydbl *opt_val;
3395
3396 Set_Update_Eigen(YES,mixt_tree->mod);
3397
3398 /* freqs = NULL; */
3399 /* n_freqs = 0; */
3400 tree = mixt_tree;
3401 lk_old = UNLIKELY;
3402 lk_new = UNLIKELY;
3403 opt_val = NULL;
3404
3405 do
3406 {
3407 if(tree->next) tree = tree->next;
3408
3409 if((tree->mod->s_opt->opt_state_freq) && (tree->io->datatype == NT))
3410 {
3411 int iter;
3412
3413 opt_val = (phydbl *)mCalloc(tree->mod->ns,sizeof(phydbl));
3414
3415 iter = 0;
3416 do
3417 {
3418 lk_old = mixt_tree->c_lnL;
3419 int failed = NO;
3420
3421 for(i=0;i<tree->mod->ns;i++) opt_val[i] = tree->mod->e_frq->pi_unscaled->v[i];
3422
3423 /* PhyML_Printf("\n. -- BFGS lk=%f",mixt_tree->c_lnL); */
3424 BFGS(mixt_tree,tree->mod->e_frq->pi_unscaled->v,tree->mod->ns,1.e-5,tree->mod->s_opt->min_diff_lk_local,1.e-2,NO,NO,
3425 &Return_Abs_Lk,
3426 &Num_Derivative_Several_Param,
3427 &Lnsrch,&failed);
3428 /* PhyML_Printf("\n. ++ BFGS lk=%f",mixt_tree->c_lnL); */
3429
3430
3431 if(failed == YES) for(i=0;i<tree->mod->ns;i++) tree->mod->e_frq->pi_unscaled->v[i] = opt_val[i];
3432
3433
3434 permut = Permutate(tree->mod->ns);
3435
3436 for(i=0;i<tree->mod->ns;++i)
3437 {
3438 phydbl a,c;
3439
3440 a = tree->mod->e_frq->pi_unscaled->v[permut[i]] / 2.;
3441 c = tree->mod->e_frq->pi_unscaled->v[permut[i]] * 2.+1;
3442
3443 /* PhyML_Printf("\n. -- lk=%f",mixt_tree->c_lnL); */
3444 Generic_Brent_Lk(&(tree->mod->e_frq->pi_unscaled->v[permut[i]]),
3445 /* UNSCALED_E_FRQ_MIN,UNSCALED_E_FRQ_MAX, */
3446 a,c,
3447 tree->mod->s_opt->min_diff_lk_local,
3448 tree->mod->s_opt->brent_it_max,
3449 tree->mod->s_opt->quickdirty,
3450 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3451 /* PhyML_Printf("\n. ++ lk=%f",mixt_tree->c_lnL); */
3452 }
3453
3454 Free(permut);
3455
3456 if(verbose)
3457 {
3458 Print_Lk(mixt_tree,"[Nucleotide freqs. ]");
3459 }
3460
3461 lk_new = mixt_tree->c_lnL;
3462
3463 /* PhyML_Printf("\n. lk_new=%f lk_old=%f",lk_new,lk_old); */
3464 assert(lk_new > lk_old - tree->mod->s_opt->min_diff_lk_local);
3465 if(fabs(lk_new-lk_old) < tree->mod->s_opt->min_diff_lk_local) break;
3466 }
3467 /* while(++iter < tree->mod->s_opt->brent_it_max); */
3468 while(++iter < 1);
3469
3470 Free(opt_val);
3471
3472 if(iter == tree->mod->s_opt->brent_it_max)
3473 {
3474 if(tree->verbose > VL0)
3475 {
3476 PhyML_Printf("\n. Failed to optimize frequency parameters this round...");
3477 }
3478 }
3479 }
3480
3481
3482 tree = tree->next;
3483
3484 }
3485 while(tree);
3486
3487
3488 Set_Update_Eigen(NO,mixt_tree->mod);
3489
3490 }
3491
3492 //////////////////////////////////////////////////////////////
3493 //////////////////////////////////////////////////////////////
3494
Optimize_Rmat_Weights(t_tree * mixt_tree,int verbose)3495 void Optimize_Rmat_Weights(t_tree *mixt_tree, int verbose)
3496 {
3497 scalar_dbl *r_mat_weight;
3498
3499 Set_Update_Eigen(NO,mixt_tree->mod);
3500
3501 if(mixt_tree->is_mixt_tree == NO) return;
3502
3503 r_mat_weight = mixt_tree->next->mod->r_mat_weight;
3504
3505 if(mixt_tree->next->mod->s_opt->opt_rmat_weight == YES)
3506 {
3507 do
3508 {
3509 phydbl a,c;
3510
3511 a = r_mat_weight->v * .1;
3512 c = r_mat_weight->v * 10.;
3513
3514 Generic_Brent_Lk(&(r_mat_weight->v),
3515 a,c,
3516 mixt_tree->mod->s_opt->min_diff_lk_local,
3517 mixt_tree->mod->s_opt->brent_it_max,
3518 mixt_tree->mod->s_opt->quickdirty,
3519 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3520
3521 if(verbose)
3522 {
3523 Print_Lk(mixt_tree,"[Rate mat. weights ]");
3524 }
3525
3526 r_mat_weight = r_mat_weight->next;
3527 }
3528 while(r_mat_weight);
3529 }
3530
3531 Set_Update_Eigen(NO,mixt_tree->mod);
3532
3533 }
3534
3535 //////////////////////////////////////////////////////////////
3536 //////////////////////////////////////////////////////////////
3537
Optimize_Efrq_Weights(t_tree * mixt_tree,int verbose)3538 void Optimize_Efrq_Weights(t_tree *mixt_tree, int verbose)
3539 {
3540 scalar_dbl *e_frq_weight;
3541
3542 Set_Update_Eigen(NO,mixt_tree->mod);
3543
3544 if(mixt_tree->is_mixt_tree == NO) return;
3545
3546 e_frq_weight = mixt_tree->next->mod->e_frq_weight;
3547
3548
3549 if(mixt_tree->next->mod->s_opt->opt_efrq_weight == YES)
3550 {
3551 do
3552 {
3553 phydbl a,c;
3554
3555 a = e_frq_weight->v * .1;
3556 c = e_frq_weight->v * 10.+1;
3557
3558 Generic_Brent_Lk(&(e_frq_weight->v),
3559 a,c,
3560 mixt_tree->mod->s_opt->min_diff_lk_local,
3561 mixt_tree->mod->s_opt->brent_it_max,
3562 mixt_tree->mod->s_opt->quickdirty,
3563 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3564
3565 if(verbose)
3566 {
3567 Print_Lk(mixt_tree,"[Equ. frq. weights ]");
3568 }
3569
3570 e_frq_weight = e_frq_weight->next;
3571 }
3572 while(e_frq_weight);
3573 }
3574
3575 Set_Update_Eigen(NO,mixt_tree->mod);
3576
3577 }
3578
3579
3580 //////////////////////////////////////////////////////////////
3581 //////////////////////////////////////////////////////////////
3582
3583
Optimize_Lambda(t_tree * mixt_tree,int verbose)3584 void Optimize_Lambda(t_tree *mixt_tree, int verbose)
3585 {
3586 scalar_dbl **lambda;
3587 int n_lambda;
3588 t_tree *tree;
3589 int i;
3590
3591 Set_Update_Eigen(YES,mixt_tree->mod);
3592
3593 lambda = NULL;
3594 n_lambda = 0;
3595 tree = mixt_tree;
3596
3597 do
3598 {
3599 if(tree->next) tree = tree->next;
3600
3601 for(i=0;i<n_lambda;i++) if(tree->mod->lambda == lambda[i]) break;
3602
3603 if(i == n_lambda)
3604 {
3605 if(!lambda) lambda = (scalar_dbl **)mCalloc(1,sizeof(scalar_dbl *));
3606 else lambda = (scalar_dbl **)mRealloc(lambda,n_lambda+1,sizeof(scalar_dbl *));
3607 lambda[n_lambda] = tree->mod->lambda;
3608 n_lambda++;
3609
3610 if(tree->mod->s_opt->opt_lambda)
3611 {
3612 Generic_Brent_Lk(&(tree->mod->lambda->v),
3613 0.001,100.,
3614 tree->mod->s_opt->min_diff_lk_local,
3615 tree->mod->s_opt->brent_it_max,
3616 tree->mod->s_opt->quickdirty,
3617 Wrap_Lk,NULL,mixt_tree,NULL,NO);
3618
3619 if(verbose)
3620 {
3621 Print_Lk(mixt_tree,"[Lambda ]");
3622 PhyML_Printf("[%10f]",tree->mod->lambda->v);
3623 }
3624 }
3625 }
3626
3627 tree = tree->next;
3628
3629 }
3630 while(tree);
3631
3632 if(lambda) Free(lambda);
3633
3634 Set_Update_Eigen(NO,mixt_tree->mod);
3635
3636 }
3637
3638 //////////////////////////////////////////////////////////////
3639 //////////////////////////////////////////////////////////////
3640
Least_Square_Node_Ages(t_tree * tree)3641 void Least_Square_Node_Ages(t_tree *tree)
3642 {
3643
3644 phydbl new_error,cur_error,sum_error;
3645 int i,j;
3646 phydbl young,old;
3647 int dir1,dir2;
3648 t_node *n;
3649
3650 Update_Ancestors(tree->n_root,tree->n_root->v[2],tree);
3651 Update_Ancestors(tree->n_root,tree->n_root->v[1],tree);
3652
3653 TIMES_Randomize_Node_Ages(tree);
3654
3655 assert(fabs(tree->times->nd_t[tree->n_root->num]) > SMALL);
3656
3657 cur_error = BIG;
3658 new_error = BIG;
3659 do
3660 {
3661 cur_error = new_error;
3662
3663 sum_error = 0.0;
3664 for(i=0;i<2*tree->n_otu-1;++i)
3665 {
3666 if(tree->a_nodes[i]->tax == NO)
3667 {
3668 n = tree->a_nodes[i];
3669 old = -1.;
3670
3671 dir1 = dir2 = -1;
3672 for(j=0;j<3;++j)
3673 {
3674 if(n->v[j] != n->anc && n->b[j] != tree->e_root)
3675 {
3676 if(dir1 < 0) dir1 = j;
3677 else dir2 = j;
3678 }
3679 else
3680 {
3681 if(n != tree->n_root) old = tree->times->nd_t[n->anc->num];
3682 else old = 2.*tree->times->nd_t[tree->n_root->num];
3683 }
3684 }
3685
3686 young = MIN(tree->times->nd_t[n->v[dir1]->num],
3687 tree->times->nd_t[n->v[dir2]->num]);
3688
3689 sum_error += Generic_Brent(tree->times->nd_t + i,
3690 young,old,1.E-10,10000,
3691 TIMES_Least_Square_Criterion,
3692 tree);
3693
3694 /* PhyML_Printf("\n. Node %3d%c time: %15f err: %15f young: %15f old: %15f %15f %15f", */
3695 /* i, */
3696 /* (n == tree->n_root)?'*':' ', */
3697 /* tree->times->nd_t[i], */
3698 /* sum_error, */
3699 /* young,old, */
3700 /* tree->times->nd_t[n->v[dir1]->num], */
3701 /* tree->times->nd_t[n->v[dir2]->num]); */
3702 }
3703 if(RATES_Check_Node_Times(tree)) Exit("\n");
3704 }
3705
3706 sum_error += Generic_Brent(&(tree->rates->clock_r),
3707 tree->rates->min_clock,
3708 tree->rates->max_clock,
3709 1.E-10,10000,
3710 TIMES_Least_Square_Criterion,
3711 tree);
3712
3713 /* PhyML_Printf("\n clock_r: %g",tree->rates->clock_r); */
3714
3715 new_error = sum_error;
3716 /* assert(new_error < cur_error+SMALL); */
3717 }
3718 while(fabs(new_error-cur_error) > 1.E-10);
3719
3720 }
3721
3722 //////////////////////////////////////////////////////////////
3723 //////////////////////////////////////////////////////////////
3724
3725 //////////////////////////////////////////////////////////////
3726 //////////////////////////////////////////////////////////////
3727
3728 //////////////////////////////////////////////////////////////
3729 //////////////////////////////////////////////////////////////
3730
3731 //////////////////////////////////////////////////////////////
3732 //////////////////////////////////////////////////////////////
3733
3734 //////////////////////////////////////////////////////////////
3735 //////////////////////////////////////////////////////////////
3736
3737 //////////////////////////////////////////////////////////////
3738 //////////////////////////////////////////////////////////////
3739
3740 //////////////////////////////////////////////////////////////
3741 //////////////////////////////////////////////////////////////
3742
3743 //////////////////////////////////////////////////////////////
3744 //////////////////////////////////////////////////////////////
3745