Context Navigation

optimiz.c

Visit:

Last change on this file was 10307, checked in by aboeckma, 12 years ago
added most recent version of phyml
File size: 88.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/GDE/PHYML20130708/phyml/src/optimiz.c

Download in other formats: