Context Navigation

lk.c

Visit:

Last change on this file was 10307, checked in by aboeckma, 11 years ago
added most recent version of phyml
File size: 116.8 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 "lk.h"
14
15	int n_sec2 = 0;
16	/* int LIM_SCALE; */
17	/* phydbl LIM_SCALE_VAL; */
18	/* phydbl BIG; */
19	/* phydbl SMALL; */
20
21	//////////////////////////////////////////////////////////////
22	//////////////////////////////////////////////////////////////
23
24
25	void Init_Tips_At_One_Site_Nucleotides_Float(char state, int pos, phydbl *p_lk)
26	{
27	switch(state)
28	{
29	case 'A' : p_lk[pos+0]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=.0;
30	break;
31	case 'C' : p_lk[pos+1]=1.; p_lk[pos+0]=p_lk[pos+2]=p_lk[pos+3]=.0;
32	break;
33	case 'G' : p_lk[pos+2]=1.; p_lk[pos+1]=p_lk[pos+0]=p_lk[pos+3]=.0;
34	break;
35	case 'T' : p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+0]=.0;
36	break;
37	case 'U' : p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+0]=.0;
38	break;
39	case 'M' : p_lk[pos+0]=p_lk[pos+1]=1.; p_lk[pos+2]=p_lk[pos+3]=.0;
40	break;
41	case 'R' : p_lk[pos+0]=p_lk[pos+2]=1.; p_lk[pos+1]=p_lk[pos+3]=.0;
42	break;
43	case 'W' : p_lk[pos+0]=p_lk[pos+3]=1.; p_lk[pos+1]=p_lk[pos+2]=.0;
44	break;
45	case 'S' : p_lk[pos+1]=p_lk[pos+2]=1.; p_lk[pos+0]=p_lk[pos+3]=.0;
46	break;
47	case 'Y' : p_lk[pos+1]=p_lk[pos+3]=1.; p_lk[pos+0]=p_lk[pos+2]=.0;
48	break;
49	case 'K' : p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+0]=p_lk[pos+1]=.0;
50	break;
51	case 'B' : p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+0]=.0;
52	break;
53	case 'D' : p_lk[pos+0]=p_lk[pos+2]=p_lk[pos+3]=1.; p_lk[pos+1]=.0;
54	break;
55	case 'H' : p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+3]=1.; p_lk[pos+2]=.0;
56	break;
57	case 'V' : p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+2]=1.; p_lk[pos+3]=.0;
58	break;
59	case 'N' : case 'X' : case '?' : case 'O' : case '-' :
60	p_lk[pos+0]=p_lk[pos+1]=p_lk[pos+2]=p_lk[pos+3]=1.;break;
61	default :
62	{
63	PhyML_Printf("\n. Unknown character state : %c\n",state);
64	Exit("\n. Init failed (check the data type)\n");
65	break;
66	}
67	}
68	}
69
70	//////////////////////////////////////////////////////////////
71	//////////////////////////////////////////////////////////////
72
73
74	void Init_Tips_At_One_Site_Nucleotides_Int(char state, int pos, short int *p_pars)
75	{
76	switch(state)
77	{
78	case 'A' : p_pars[pos+0]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=0;
79	break;
80	case 'C' : p_pars[pos+1]=1; p_pars[pos+0]=p_pars[pos+2]=p_pars[pos+3]=0;
81	break;
82	case 'G' : p_pars[pos+2]=1; p_pars[pos+1]=p_pars[pos+0]=p_pars[pos+3]=0;
83	break;
84	case 'T' : p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+0]=0;
85	break;
86	case 'U' : p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+0]=0;
87	break;
88	case 'M' : p_pars[pos+0]=p_pars[pos+1]=1; p_pars[pos+2]=p_pars[pos+3]=0;
89	break;
90	case 'R' : p_pars[pos+0]=p_pars[pos+2]=1; p_pars[pos+1]=p_pars[pos+3]=0;
91	break;
92	case 'W' : p_pars[pos+0]=p_pars[pos+3]=1; p_pars[pos+1]=p_pars[pos+2]=0;
93	break;
94	case 'S' : p_pars[pos+1]=p_pars[pos+2]=1; p_pars[pos+0]=p_pars[pos+3]=0;
95	break;
96	case 'Y' : p_pars[pos+1]=p_pars[pos+3]=1; p_pars[pos+0]=p_pars[pos+2]=0;
97	break;
98	case 'K' : p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+0]=p_pars[pos+1]=0;
99	break;
100	case 'B' : p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+0]=0;
101	break;
102	case 'D' : p_pars[pos+0]=p_pars[pos+2]=p_pars[pos+3]=1; p_pars[pos+1]=0;
103	break;
104	case 'H' : p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+3]=1; p_pars[pos+2]=0;
105	break;
106	case 'V' : p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+2]=1; p_pars[pos+3]=0;
107	break;
108	case 'N' : case 'X' : case '?' : case 'O' : case '-' :
109	p_pars[pos+0]=p_pars[pos+1]=p_pars[pos+2]=p_pars[pos+3]=1;break;
110	default :
111	{
112	PhyML_Printf("\n. Unknown character state : %c\n",state);
113	Exit("\n. Init failed (check the data type)\n");
114	break;
115	}
116	}
117	}
118
119	//////////////////////////////////////////////////////////////
120	//////////////////////////////////////////////////////////////
121
122
123	void Init_Tips_At_One_Site_AA_Float(char aa, int pos, phydbl *p_lk)
124	{
125	int i;
126
127	For(i,20) p_lk[pos+i] = .0;
128
129	switch(aa){
130	case 'A' : p_lk[pos+0]= 1.; break;/* Alanine */
131	case 'R' : p_lk[pos+1]= 1.; break;/* Arginine */
132	case 'N' : p_lk[pos+2]= 1.; break;/* Asparagine */
133	case 'D' : p_lk[pos+3]= 1.; break;/* Aspartic acid */
134	case 'C' : p_lk[pos+4]= 1.; break;/* Cysteine */
135	case 'Q' : p_lk[pos+5]= 1.; break;/* Glutamine */
136	case 'E' : p_lk[pos+6]= 1.; break;/* Glutamic acid */
137	case 'G' : p_lk[pos+7]= 1.; break;/* Glycine */
138	case 'H' : p_lk[pos+8]= 1.; break;/* Histidine */
139	case 'I' : p_lk[pos+9]= 1.; break;/* Isoleucine */
140	case 'L' : p_lk[pos+10]=1.; break;/* Leucine */
141	case 'K' : p_lk[pos+11]=1.; break;/* Lysine */
142	case 'M' : p_lk[pos+12]=1.; break;/* Methionine */
143	case 'F' : p_lk[pos+13]=1.; break;/* Phenylalanin */
144	case 'P' : p_lk[pos+14]=1.; break;/* Proline */
145	case 'S' : p_lk[pos+15]=1.; break;/* Serine */
146	case 'T' : p_lk[pos+16]=1.; break;/* Threonine */
147	case 'W' : p_lk[pos+17]=1.; break;/* Tryptophan */
148	case 'Y' : p_lk[pos+18]=1.; break;/* Tyrosine */
149	case 'V' : p_lk[pos+19]=1.; break;/* Valine */
150
151	case 'B' : p_lk[pos+2]= 1.; break;/* Asparagine */
152	case 'Z' : p_lk[pos+5]= 1.; break;/* Glutamine */
153
154	case 'X' : case '?' : case '-' : For(i,20) p_lk[pos+i] = 1.; break;
155	default :
156	{
157	PhyML_Printf("\n. Unknown character state : %c\n",aa);
158	Exit("\n. Init failed (check the data type)\n");
159	break;
160	}
161	}
162	}
163
164	//////////////////////////////////////////////////////////////
165	//////////////////////////////////////////////////////////////
166
167
168	void Init_Tips_At_One_Site_AA_Int(char aa, int pos, short int *p_pars)
169	{
170	int i;
171
172	For(i,20) p_pars[pos+i] = .0;
173
174	switch(aa){
175	case 'A' : p_pars[pos+0] = 1; break;/* Alanine */
176	case 'R' : p_pars[pos+1] = 1; break;/* Arginine */
177	case 'N' : p_pars[pos+2] = 1; break;/* Asparagine */
178	case 'D' : p_pars[pos+3] = 1; break;/* Aspartic acid */
179	case 'C' : p_pars[pos+4] = 1; break;/* Cysteine */
180	case 'Q' : p_pars[pos+5] = 1; break;/* Glutamine */
181	case 'E' : p_pars[pos+6] = 1; break;/* Glutamic acid */
182	case 'G' : p_pars[pos+7] = 1; break;/* Glycine */
183	case 'H' : p_pars[pos+8] = 1; break;/* Histidine */
184	case 'I' : p_pars[pos+9] = 1; break;/* Isoleucine */
185	case 'L' : p_pars[pos+10] = 1; break;/* Leucine */
186	case 'K' : p_pars[pos+11] = 1; break;/* Lysine */
187	case 'M' : p_pars[pos+12] = 1; break;/* Methionine */
188	case 'F' : p_pars[pos+13] = 1; break;/* Phenylalanin */
189	case 'P' : p_pars[pos+14] = 1; break;/* Proline */
190	case 'S' : p_pars[pos+15] = 1; break;/* Serine */
191	case 'T' : p_pars[pos+16] = 1; break;/* Threonine */
192	case 'W' : p_pars[pos+17] = 1; break;/* Tryptophan */
193	case 'Y' : p_pars[pos+18] = 1; break;/* Tyrosine */
194	case 'V' : p_pars[pos+19] = 1; break;/* Valine */
195
196	case 'B' : p_pars[pos+2] = 1; break;/* Asparagine */
197	case 'Z' : p_pars[pos+5] = 1; break;/* Glutamine */
198
199	case 'X' : case '?' : case '-' : For(i,20) p_pars[pos+i] = 1; break;
200	default :
201	{
202	PhyML_Printf("\n. Unknown character state : %c\n",aa);
203	Exit("\n. Init failed (check the data type)\n");
204	break;
205	}
206	}
207	}
208
209	//////////////////////////////////////////////////////////////
210	//////////////////////////////////////////////////////////////
211
212
213	void Init_Tips_At_One_Site_Generic_Float(char state, int ns, int state_len, int pos, phydbl p_lk)
214	{
215	int i;
216	int state_int;
217
218	For(i,ns) p_lk[pos+i] = 0.;
219
220	if(Is_Ambigu(state,GENERIC,state_len)) For(i,ns) p_lk[pos+i] = 1.;
221	else
222	{
223	char format[6];
224	sprintf(format,"%%%dd",state_len);
225	if(!sscanf(state,format,&state_int))
226	{
227	PhyML_Printf("\n. state='%c'",state);
228	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
229	Warn_And_Exit("");
230	}
231	if(state_int > ns)
232	{
233	PhyML_Printf("\n. %s %d cstate: %.2s istate: %d state_len: %d.\n",__FILE__,__LINE__,state,state_int,state_len);
234	PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
235	Warn_And_Exit("");
236	}
237	p_lk[pos+state_int] = 1.;
238	/* PhyML_Printf("\n. %s %d cstate: %.2s istate: %d state_len: %d ns: %d pos: %d",__FILE__,__LINE__,state,state_int,state_len,ns,pos); */
239	}
240	}
241
242	//////////////////////////////////////////////////////////////
243	//////////////////////////////////////////////////////////////
244
245
246	void Init_Tips_At_One_Site_Generic_Int(char state, int ns, int state_len, int pos, short int p_pars)
247	{
248	int i;
249	int state_int;
250
251	For(i,ns) p_pars[pos+i] = 0;
252
253	if(Is_Ambigu(state,GENERIC,state_len)) For(i,ns) p_pars[pos+i] = 1;
254	else
255	{
256	char format[6];
257	sprintf(format,"%%%dd",state_len);
258	if(!sscanf(state,format,&state_int))
259	{
260	PhyML_Printf("\n. state='%c'",state);
261	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
262	Warn_And_Exit("");
263	}
264	if(state_int > ns)
265	{
266	PhyML_Printf("\n. %s %d cstate: %.2s istate: %d state_len: %d.\n",__FILE__,__LINE__,state,state_int,state_len);
267	PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
268	Warn_And_Exit("");
269	}
270	p_pars[pos+state_int] = 1;
271	/* PhyML_Printf("\n* %s %d cstate: %.2s istate: %d state_len: %d ns: %d pos: %d",__FILE__,__LINE__,state,state_int,state_len,ns,pos); */
272	}
273	}
274
275	//////////////////////////////////////////////////////////////
276	//////////////////////////////////////////////////////////////
277
278
279	void Get_All_Partial_Lk_Scale(t_tree tree, t_edge b_fcus, t_node a, t_node d)
280	{
281	Update_P_Lk(tree,b_fcus,d);
282	}
283
284	//////////////////////////////////////////////////////////////
285	//////////////////////////////////////////////////////////////
286
287	/* void Post_Order_Lk(t_node a, t_node d, t_tree tree) /
288	/* { */
289	/* int i,dir; */
290
291	/* if(tree->is_mixt_tree == YES) */
292	/* { */
293	/* MIXT_Post_Order_Lk(a,d,tree); */
294	/* return; */
295	/* } */
296
297	/* dir = -1; */
298
299	/* if(d->tax) */
300	/* { */
301	/* Get_All_Partial_Lk_Scale(tree,d->b[0],a,d); */
302	/* return; */
303	/* } */
304	/* else */
305	/* { */
306	/* For(i,3) */
307	/* { */
308	/* if(d->v[i] != a) */
309	/* Post_Order_Lk(d,d->v[i],tree); */
310	/* else dir = i; */
311	/* } */
312	/* Get_All_Partial_Lk_Scale(tree,d->b[dir],a,d); */
313	/* } */
314	/* } */
315
316	//////////////////////////////////////////////////////////////
317	//////////////////////////////////////////////////////////////
318
319	void Post_Order_Lk(t_node a, t_node d, t_tree *tree)
320	{
321	int i,dir;
322
323	dir = -1;
324
325	if(d->tax) return;
326	else
327	{
328	if(tree->is_mixt_tree)
329	{
330	MIXT_Post_Order_Lk(a,d,tree);
331	return;
332	}
333
334	For(i,3)
335	{
336	if(d->v[i] != a && d->b[i] != tree->e_root)
337	Post_Order_Lk(d,d->v[i],tree);
338	else dir = i;
339	}
340
341	if(d->b[dir] != tree->e_root)
342	Get_All_Partial_Lk_Scale(tree,d->b[dir],a,d);
343	else
344	{
345	if(d == tree->n_root->v[1]) Get_All_Partial_Lk_Scale(tree,tree->n_root->b[1],tree->n_root,d);
346	else Get_All_Partial_Lk_Scale(tree,tree->n_root->b[2],tree->n_root,d);
347	}
348	}
349	}
350
351	//////////////////////////////////////////////////////////////
352	//////////////////////////////////////////////////////////////
353
354	void Pre_Order_Lk(t_node a, t_node d, t_tree *tree)
355	{
356	int i;
357
358	if(d->tax) return;
359	else
360	{
361	if(tree->is_mixt_tree)
362	{
363	MIXT_Pre_Order_Lk(a,d,tree);
364	return;
365	}
366
367	For(i,3)
368	{
369	if(d->v[i] != a && d->b[i] != tree->e_root)
370	{
371	Get_All_Partial_Lk_Scale(tree,d->b[i],d->v[i],d);
372	Pre_Order_Lk(d,d->v[i],tree);
373	}
374	}
375	}
376	}
377
378	//////////////////////////////////////////////////////////////
379	//////////////////////////////////////////////////////////////
380
381	phydbl Lk(t_edge b, t_tree tree)
382	{
383	int br;
384	int n_patterns,ambiguity_check,state;
385
386	if(b == NULL && tree->mod->s_opt->curr_opt_free_rates == YES)
387	{
388	tree->mod->s_opt->curr_opt_free_rates = NO;
389	Optimize_Free_Rate_Weights(tree,YES,YES);
390	tree->mod->s_opt->curr_opt_free_rates = YES;
391	}
392
393
394	if(tree->is_mixt_tree) return MIXT_Lk(b,tree);
395
396	tree->old_lnL = tree->c_lnL;
397
398	#if (defined PHYTIME \|\| defined SERGEII)
399	if((tree->rates) && (tree->rates->bl_from_rt)) RATES_Update_Cur_Bl(tree);
400	#endif
401
402
403	if(tree->rates && tree->io->lk_approx == NORMAL)
404	{
405	tree->c_lnL = Lk_Normal_Approx(tree);
406	return tree->c_lnL;
407	}
408
409	n_patterns = tree->n_pattern;
410
411	if(!b) Set_Model_Parameters(tree->mod);
412
413	Set_Br_Len_Var(tree);
414
415	if(tree->mod->s_opt->skip_tree_traversal == NO)
416	{
417	if(!b)
418	{
419	For(br,2*tree->n_otu-3) Update_PMat_At_Given_Edge(tree->a_edges[br],tree);
420	if(tree->n_root)
421	{
422	Update_PMat_At_Given_Edge(tree->n_root->b[1],tree);
423	Update_PMat_At_Given_Edge(tree->n_root->b[2],tree);
424	}
425	}
426	else
427	{
428	Update_PMat_At_Given_Edge(b,tree);
429	}
430
431	if(!b)
432	{
433	if(tree->n_root)
434	{
435	Post_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
436	Post_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
437
438	Update_P_Lk(tree,tree->n_root->b[1],tree->n_root);
439	Update_P_Lk(tree,tree->n_root->b[2],tree->n_root);
440
441	if(tree->both_sides == YES)
442	Pre_Order_Lk(tree->n_root,tree->n_root->v[1],tree);
443
444	if(tree->both_sides == YES)
445	Pre_Order_Lk(tree->n_root,tree->n_root->v[2],tree);
446	}
447	else
448	{
449	Post_Order_Lk(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
450	if(tree->both_sides == YES)
451	Pre_Order_Lk(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
452	}
453	}
454	}
455
456	if(!b)
457	{
458	if(tree->n_root) b = (tree->n_root->v[1]->tax == NO)?(tree->n_root->b[2]):(tree->n_root->b[1]);
459	else b = tree->a_nodes[0]->b[0];
460	}
461
462	tree->c_lnL = .0;
463	tree->sum_min_sum_scale = .0;
464
465	For(tree->curr_site,n_patterns)
466	{
467	ambiguity_check = -1;
468	state = -1;
469
470	if(tree->data->wght[tree->curr_site] > SMALL)
471	{
472	if((b->rght->tax) && (!tree->mod->s_opt->greedy))
473	{
474	ambiguity_check = b->rght->c_seq->is_ambigu[tree->curr_site];
475	if(!ambiguity_check)
476	{
477	state = b->rght->c_seq->d_state[tree->curr_site];
478	}
479	}
480
481	if(tree->mod->use_m4mod) ambiguity_check = YES;
482
483	Lk_Core(state,ambiguity_check,b,tree);
484	}
485	}
486
487	/* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); */
488
489	/* tree->c_lnL = .0; */
490	/* For(tree->curr_site,n_patterns) */
491	/* { */
492	/* tree->c_lnL += tree->c_lnL_sorted[tree->curr_site]; */
493	/* } */
494
495
496	Adjust_Min_Diff_Lk(tree);
497
498	return tree->c_lnL;
499	}
500
501	//////////////////////////////////////////////////////////////
502	//////////////////////////////////////////////////////////////
503
504	/* Core of the likelihood calculcation. Assume that the partial likelihoods on both
505	sides of t_edge *b are up-to-date. Calculate the log-likelihood at one site.
506	*/
507
508	phydbl Lk_Core(int state, int ambiguity_check, t_edge b, t_tree tree)
509	{
510	phydbl log_site_lk;
511	phydbl site_lk_cat, site_lk, inv_site_lk;
512	int fact_sum_scale;
513	phydbl max_sum_scale,min_sum_scale;
514	phydbl sum;
515	int catg,ns,k,l,site;
516	int dim1,dim2,dim3;
517	int sum_scale_left_cat,sum_scale_rght_cat;
518	int exponent;
519	int num_prec_issue;
520	phydbl tmp;
521
522	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
523	dim2 = tree->mod->ns;
524	dim3 = tree->mod->ns * tree->mod->ns;
525
526	log_site_lk = .0;
527	site_lk = .0;
528	site_lk_cat = .0;
529	site = tree->curr_site;
530	ns = tree->mod->ns;
531
532	sum_scale_left_cat = b->sum_scale_left_cat;
533	sum_scale_rght_cat = b->sum_scale_rght_cat;
534
535	/* Skip this if no tree traveral was required, i.e. likelihood is already up to date */
536	if(tree->mod->s_opt->skip_tree_traversal == NO)
537	{
538	/* Actual likelihood calculation */
539	/* For all classes of rates */
540	For(catg,tree->mod->ras->n_catg)
541	{
542	site_lk_cat = .0;
543
544	/* b is an external edge */
545	if((b->rght->tax) && (!tree->mod->s_opt->greedy))
546	{
547	/* If the character observed at the tip is NOT ambiguous: ns x 1 terms to consider */
548	if(!ambiguity_check)
549	{
550	sum = .0;
551	For(l,ns)
552	{
553	sum +=
554	b->Pij_rr[catgdim3+statedim2+l] *
555	b->p_lk_left[sitedim1+catgdim2+l];
556	}
557
558	site_lk_cat += sum * tree->mod->e_frq->pi->v[state];
559	}
560	/* If the character observed at the tip is ambiguous: ns x ns terms to consider */
561	else
562	{
563	For(k,ns)
564	{
565	sum = .0;
566	if(b->p_lk_tip_r[site*dim2+k] > .0)
567	{
568	For(l,ns)
569	{
570	sum +=
571	b->Pij_rr[catgdim3+kdim2+l] *
572	b->p_lk_left[sitedim1+catgdim2+l];
573	}
574
575	site_lk_cat +=
576	sum *
577	tree->mod->e_frq->pi->v[k] *
578	b->p_lk_tip_r[site*dim2+k];
579	}
580	}
581	}
582	}
583	/* b is an internal edge: ns x ns terms to consider */
584	else
585	{
586	For(k,ns)
587	{
588	sum = .0;
589	if(b->p_lk_rght[sitedim1+catgdim2+k] > .0)
590	{
591	For(l,ns)
592	{
593	sum +=
594	b->Pij_rr[catgdim3+kdim2+l] *
595	b->p_lk_left[sitedim1+catgdim2+l];
596	}
597
598	site_lk_cat +=
599	sum *
600	tree->mod->e_frq->pi->v[k] *
601	b->p_lk_rght[sitedim1+catgdim2+k];
602	}
603	}
604	}
605	tree->site_lk_cat[catg] = site_lk_cat;
606	}
607
608	if(tree->apply_lk_scaling == YES)
609	{
610	max_sum_scale = (phydbl)BIG;
611	min_sum_scale = -(phydbl)BIG;
612
613	For(catg,tree->mod->ras->n_catg)
614	{
615	sum_scale_left_cat[catg] =
616	(b->sum_scale_left)?
617	(b->sum_scale_left[catg*tree->n_pattern+site]):
618	(0.0);
619
620	sum_scale_rght_cat[catg] =
621	(b->sum_scale_rght)?
622	(b->sum_scale_rght[catg*tree->n_pattern+site]):
623	(0.0);
624
625	sum = sum_scale_left_cat[catg] + sum_scale_rght_cat[catg];
626
627	if(sum < .0)
628	{
629	PhyML_Printf("\n== b->num = %d sum = %G",sum,b->num);
630	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
631	Exit("\n");
632	}
633
634	tmp = sum + ((phydbl)LOGBIG - LOG(tree->site_lk_cat[catg]))/(phydbl)LOG2;
635	if(tmp < max_sum_scale) max_sum_scale = tmp; /* min of the maxs */
636
637	tmp = sum + ((phydbl)LOGSMALL - LOG(tree->site_lk_cat[catg]))/(phydbl)LOG2;
638	if(tmp > min_sum_scale) min_sum_scale = tmp; /* max of the mins */
639
640	}
641
642	if(min_sum_scale > max_sum_scale)
643	{
644	/* PhyML_Printf("\n== Numerical precision issue alert."); */
645	/* PhyML_Printf("\n== min_sum_scale = %G max_sum_scale = %G",min_sum_scale,max_sum_scale); */
646	/* PhyML_Printf("\n== Err in file %s at line %d\n\n",__FILE__,__LINE__); */
647	/* Warn_And_Exit("\n"); */
648	min_sum_scale = max_sum_scale;
649	}
650
651	fact_sum_scale = (int)((max_sum_scale + min_sum_scale) / 2);
652
653
654	/* fact_sum_scale = (int)(max_sum_scale / 2); */
655
656	/* Apply scaling factors */
657	For(catg,tree->mod->ras->n_catg)
658	{
659	exponent = -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale;
660	site_lk_cat = tree->site_lk_cat[catg];
661	Rate_Correction(exponent,&site_lk_cat,tree);
662	tree->site_lk_cat[catg] = site_lk_cat;
663	}
664	}
665	else // No scaling of lk
666	{
667	fact_sum_scale = 0;
668	}
669
670
671	site_lk = .0;
672	For(catg,tree->mod->ras->n_catg)
673	{
674	site_lk +=
675	tree->site_lk_cat[catg] *
676	tree->mod->ras->gamma_r_proba->v[catg];
677
678	}
679
680	if(tree->mod->ras->invar == YES)
681	{
682	num_prec_issue = NO;
683	inv_site_lk = Invariant_Lk(&fact_sum_scale,site,&num_prec_issue,tree);
684
685	if(num_prec_issue == YES) // inv_site_lk >> site_lk
686	{
687	site_lk = inv_site_lk * tree->mod->ras->pinvar->v;
688	}
689	else
690	{
691	site_lk = site_lk * (1. - tree->mod->ras->pinvar->v) + inv_site_lk * tree->mod->ras->pinvar->v;
692	}
693	}
694
695	log_site_lk = LOG(site_lk) - (phydbl)LOG2 * fact_sum_scale;
696
697	int piecewise_exponent;
698	phydbl multiplier;
699	if(fact_sum_scale >= 0)
700	{
701	tree->cur_site_lk[site] = site_lk;
702	exponent = -fact_sum_scale;
703	do
704	{
705	piecewise_exponent = MAX(exponent,-63);
706	multiplier = 1. / (phydbl)((unsigned long long)(1) << -piecewise_exponent);
707	tree->cur_site_lk[site] *= multiplier;
708	exponent -= piecewise_exponent;
709	}
710	while(exponent != 0);
711	}
712	else
713	{
714	tree->cur_site_lk[site] = site_lk;
715	exponent = fact_sum_scale;
716	do
717	{
718	piecewise_exponent = MIN(exponent,63);
719	multiplier = (phydbl)((unsigned long long)(1) << piecewise_exponent);
720	tree->cur_site_lk[site] *= multiplier;
721	exponent -= piecewise_exponent;
722	}
723	while(exponent != 0);
724	}
725
726	/* tree->cur_site_lk[site] = EXP(log_site_lk); */
727
728	For(catg,tree->mod->ras->n_catg) tree->log_site_lk_cat[catg][site] = LOG(tree->site_lk_cat[catg]) - (phydbl)LOG2 * fact_sum_scale;
729
730	if(isinf(log_site_lk) \|\| isnan(log_site_lk))
731	{
732	PhyML_Printf("\n== Site = %d",site);
733	PhyML_Printf("\n== Invar = %d",tree->data->invar[site]);
734	PhyML_Printf("\n== Mixt = %d",tree->is_mixt_tree);
735	PhyML_Printf("\n== Lk = %G LOG(Lk) = %f < %G",site_lk,log_site_lk,-BIG);
736	For(catg,tree->mod->ras->n_catg) PhyML_Printf("\n== rr=%f p=%f",tree->mod->ras->gamma_rr->v[catg],tree->mod->ras->gamma_r_proba->v[catg]);
737	PhyML_Printf("\n== Pinv = %G",tree->mod->ras->pinvar->v);
738	PhyML_Printf("\n== Bl mult = %G",tree->mod->br_len_multiplier->v);
739
740	/* int i; */
741	/* For(i,2tree->n_otu-3) /
742	/* { */
743	/* PhyML_Printf("\n. b%3d->l->v = %f %f [%G] %f %f %f %f [%s]",i, */
744	/* tree->a_edges[i]->l->v, */
745	/* tree->a_edges[i]->gamma_prior_mean, */
746	/* tree->a_edges[i]->gamma_prior_var, */
747	/* tree->rates->nd_t[tree->a_edges[i]->left->num], */
748	/* tree->rates->nd_t[tree->a_edges[i]->rght->num], */
749	/* tree->rates->br_r[tree->a_edges[i]->left->num], */
750	/* tree->rates->br_r[tree->a_edges[i]->rght->num], */
751	/* tree->a_edges[i] == tree->e_root ? "YES" : "NO"); */
752	/* fflush(NULL); */
753
754	/* } */
755
756	PhyML_Printf("\n== Err. in file %s at line %d",__FILE__,__LINE__);
757	Exit("\n");
758	}
759	}
760	else
761	{
762	site_lk = .0;
763	For(catg,tree->mod->ras->n_catg)
764	site_lk +=
765	EXP(tree->log_site_lk_cat[catg][site])*
766	tree->mod->ras->gamma_r_proba->v[catg];
767
768	tree->cur_site_lk[site] = site_lk;
769	log_site_lk = LOG(site_lk);
770	}
771
772	/* Multiply log likelihood by the number of times this site pattern is found in the data */
773	tree->c_lnL_sorted[site] = tree->data->wght[site]*log_site_lk;
774
775	tree->c_lnL += tree->data->wght[site]*log_site_lk;
776
777	return log_site_lk;
778	}
779
780	//////////////////////////////////////////////////////////////
781	//////////////////////////////////////////////////////////////
782
783	/* phydbl Lk_At_Given_Edge(t_edge b_fcus, t_tree tree) */
784	/* { */
785	/* int n_patterns; */
786
787	/* if(tree->is_mixt_tree) return MIXT_Lk_At_Given_Edge(tree); */
788
789	/* n_patterns = tree->n_pattern; */
790
791	/* #ifdef PHYTIME */
792	/* if((tree->rates) && (tree->rates->bl_from_rt)) RATES_Update_Cur_Bl(tree); */
793	/* #endif */
794
795	/* Check_Br_Len_Bounds(tree); */
796
797	/* if(tree->rates && tree->io->lk_approx == NORMAL) */
798	/* { */
799	/* tree->c_lnL = Lk_Normal_Approx(tree); */
800	/* return tree->c_lnL; */
801	/* } */
802
803	/* Update_PMat_At_Given_Edge(b_fcus,tree); */
804
805	/* if(b_fcus->left->tax) */
806	/* { */
807	/* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
808	/* Warn_And_Exit(""); */
809	/* } */
810
811	/* tree->c_lnL = .0; */
812	/* tree->sum_min_sum_scale = .0; */
813	/* For(tree->curr_site,n_patterns) */
814	/* if(tree->data->wght[tree->curr_site] > SMALL) */
815	/* Lk_Core(b_fcus,tree); */
816
817	/* /\* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); \/ /
818
819	/* /\* tree->c_lnL = .0; \/ /
820	/* /\* For(tree->curr_site,n_patterns) \/ /
821	/* /\* { \/ /
822	/* /\* tree->c_lnL += tree->c_lnL_sorted[tree->curr_site]; \/ /
823	/* /\* } \/ /
824
825
826	/* Adjust_Min_Diff_Lk(tree); */
827
828	/* return tree->c_lnL; */
829	/* } */
830
831	//////////////////////////////////////////////////////////////
832	//////////////////////////////////////////////////////////////
833
834	/////////////////////////////////////////////////////////////
835	//////////////////////////////////////////////////////////////
836
837	void Rate_Correction(int exponent, phydbl site_lk_cat, t_tree tree)
838	{
839	int piecewise_exponent;
840	phydbl multiplier;
841
842	if(exponent >= 0)
843	{
844	/! Multiply by 2^exponent /
845	do
846	{
847	piecewise_exponent = MIN(exponent,63);
848	multiplier = (phydbl)((unsigned long long)(1) << piecewise_exponent);
849	(site_lk_cat) = multiplier;
850	exponent -= piecewise_exponent;
851	}
852	while(exponent != 0);
853	}
854	else
855	{
856	do
857	{
858	piecewise_exponent = MAX(exponent,-63);
859	multiplier = 1. / (phydbl)((unsigned long long)(1) << -piecewise_exponent);
860	(site_lk_cat) = multiplier;
861	exponent -= piecewise_exponent;
862	}
863	while(exponent != 0);
864	}
865
866	if(isinf(*site_lk_cat))
867	{
868	PhyML_Printf("\n== Numerical precision issue alert.");
869	PhyML_Printf("\n== exponent: %d site_lk_cat: %f", exponent,*site_lk_cat);
870	PhyML_Printf("\n== File %s at line %d\n\n",__FILE__,__LINE__);
871	(*site_lk_cat) = BIG / 10;
872	}
873
874	if(*site_lk_cat < SMALL)
875	{
876	if(tree->mod->ras->n_catg == 1)
877	{
878	PhyML_Printf("\n== site_lk_cat = %G",*site_lk_cat);
879	PhyML_Printf("\n== exponent: %d", exponent);
880	PhyML_Printf("\n== Numerical precision issue alert.");
881	PhyML_Printf("\n== File %s at line %d\n\n",__FILE__,__LINE__);
882	Exit("\n");
883	}
884	(*site_lk_cat) = .0;
885	}
886	}
887	//////////////////////////////////////////////////////////////
888	//////////////////////////////////////////////////////////////
889
890	// Returns the scaled likelihood for invariable sites
891	phydbl Invariant_Lk(int fact_sum_scale, int site, int num_prec_issue, t_tree *tree)
892	{
893	int exponent,piecewise_exponent;
894	phydbl multiplier;
895	phydbl inv_site_lk = 0.;
896
897	(*num_prec_issue) = NO;
898
899	/* The substitution model does include invariable sites */
900	if(tree->mod->ras->invar == YES)
901	{
902	/* The site is invariant */
903	if(tree->data->invar[site] > -0.5)
904	{
905	inv_site_lk = tree->mod->e_frq->pi->v[tree->data->invar[site]];
906
907	/* printf("\n. inv_site_lk = %f [%c] [%d]",inv_site_lk,tree->data->c_seq[0]->state[site],tree->data->invar[site]); */
908
909	if(tree->apply_lk_scaling == YES)
910	{
911	exponent = (*fact_sum_scale);
912	do
913	{
914	piecewise_exponent = MIN(exponent,63);
915	multiplier = (phydbl)((unsigned long long)(1) << piecewise_exponent);
916	inv_site_lk *= multiplier;
917	exponent -= piecewise_exponent;
918	}
919	while(exponent != 0);
920	}
921
922	/* Update the value of site_lk */
923	if(isinf(inv_site_lk)) // P(D\|r=0) >> P(D\|r>0) => assume P(D) = P(D\|r=0)P(r=0)
924	{
925	PhyML_Printf("\n== Numerical precision issue alert.");
926	PhyML_Printf("\n== File %s at line %d\n\n",__FILE__,__LINE__);
927	(*num_prec_issue) = YES;
928	}
929	}
930	}
931
932	return inv_site_lk;
933
934	}
935
936	//////////////////////////////////////////////////////////////
937	//////////////////////////////////////////////////////////////
938
939	/* Update partial likelihood on edge b on the side of b where
940	node d lies.
941	*/
942
943	void Update_P_Lk(t_tree tree, t_edge b, t_node *d)
944	{
945
946	if(tree->is_mixt_tree)
947	{
948	MIXT_Update_P_Lk(tree,b,d);
949	return;
950	}
951
952	if((tree->io->do_alias_subpatt == YES) &&
953	(tree->update_alias_subpatt == YES))
954	Alias_One_Subpatt((d==b->left)?(b->rght):(b->left),d,tree);
955
956	if(d->tax) return;
957
958	if(tree->mod->use_m4mod == NO)
959	{
960	if(tree->io->datatype == NT)
961	{
962	Update_P_Lk_Nucl(tree,b,d);
963	}
964	else if(tree->io->datatype == AA)
965	{
966	Update_P_Lk_AA(tree,b,d);
967	}
968	else
969	{
970	Update_P_Lk_Generic(tree,b,d);
971	}
972	}
973	else
974	{
975	Update_P_Lk_Generic(tree,b,d);
976	}
977	}
978
979	//////////////////////////////////////////////////////////////
980	//////////////////////////////////////////////////////////////
981
982	void Update_P_Lk_Generic(t_tree tree, t_edge b, t_node *d)
983	{
984	/*
985	\|
986	\|<- b
987	\|
988	d
989	/ \
990	/ \
991	/ \
992	n_v1 n_v2
993	*/
994	t_node n_v1, n_v2;
995	phydbl p1_lk1,p2_lk2;
996	phydbl p_lk,p_lk_v1,*p_lk_v2;
997	phydbl Pij1,Pij2;
998	int sum_scale, sum_scale_v1, *sum_scale_v2;
999	int sum_scale_v1_val, sum_scale_v2_val;
1000	int i,j;
1001	int catg,site;
1002	int n_patterns;
1003	short int ambiguity_check_v1,ambiguity_check_v2;
1004	int state_v1,state_v2;
1005	int NsNg, Ns, NsNs;
1006	phydbl curr_scaler;
1007	int curr_scaler_pow, piecewise_scaler_pow;
1008	phydbl p_lk_lim_inf;
1009	phydbl smallest_p_lk;
1010	int *p_lk_loc;
1011
1012	p_lk_lim_inf = (phydbl)P_LK_LIM_INF;
1013
1014	NsNg = tree->mod->ras->n_catg * tree->mod->ns;
1015	Ns = tree->mod->ns;
1016	NsNs = tree->mod->ns * tree->mod->ns;
1017
1018	state_v1 = state_v2 = -1;
1019	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1020	sum_scale_v1_val = sum_scale_v2_val = 0;
1021	p1_lk1 = p2_lk2 = .0;
1022
1023	if(d->tax)
1024	{
1025	PhyML_Printf("\n== t_node %d is a leaf...",d->num);
1026	PhyML_Printf("\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
1027	Warn_And_Exit("\n");
1028	}
1029
1030	n_patterns = tree->n_pattern;
1031
1032	n_v1 = n_v2 = NULL;
1033	p_lk = p_lk_v1 = p_lk_v2 = NULL;
1034	Pij1 = Pij2 = NULL;
1035	sum_scale_v1 = sum_scale_v2 = NULL;
1036	p_lk_loc = NULL;
1037
1038	Set_All_P_Lk(&n_v1,&n_v2,
1039	&p_lk,&sum_scale,&p_lk_loc,
1040	&Pij1,&p_lk_v1,&sum_scale_v1,
1041	&Pij2,&p_lk_v2,&sum_scale_v2,
1042	d,b,tree);
1043
1044	/* For every site in the alignment */
1045	For(site,n_patterns)
1046	{
1047	state_v1 = state_v2 = -1;
1048	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1049
1050	if(!tree->mod->s_opt->greedy)
1051	{
1052	/* n_v1 and n_v2 are tip nodes */
1053	if(n_v1 && n_v1->tax)
1054	{
1055	/* Is the state at this tip ambiguous? */
1056	ambiguity_check_v1 = n_v1->c_seq->is_ambigu[site];
1057	/* if(ambiguity_check_v1 == NO) state_v1 = Get_State_From_P_Pars(n_v1->b[0]->p_lk_tip_r,siteNs,tree); /
1058	if(ambiguity_check_v1 == NO) state_v1 = n_v1->c_seq->d_state[site];
1059	}
1060
1061	if(n_v2 && n_v2->tax)
1062	{
1063	/* Is the state at this tip ambiguous? */
1064	ambiguity_check_v2 = n_v2->c_seq->is_ambigu[site];
1065	/* ambiguity_check_v2 = tree->data->c_seq[n_v2->num]->is_ambigu[site]; */
1066	/* if(ambiguity_check_v2 == NO) state_v2 = Get_State_From_P_Pars(n_v2->b[0]->p_lk_tip_r,siteNs,tree); /
1067	if(ambiguity_check_v2 == NO) state_v2 = n_v2->c_seq->d_state[site];
1068	}
1069	}
1070
1071	if(tree->mod->use_m4mod)
1072	{
1073	ambiguity_check_v1 = YES;
1074	ambiguity_check_v2 = YES;
1075	}
1076
1077	if(p_lk_loc[site] < site)
1078	{
1079	Copy_P_Lk(p_lk,p_lk_loc[site],site,tree);
1080	Copy_Scale(sum_scale,p_lk_loc[site],site,tree);
1081	}
1082	else
1083	{
1084	/* For all the rate classes */
1085	For(catg,tree->mod->ras->n_catg)
1086	{
1087	if(tree->mod->ras->skip_rate_cat[catg] == YES) continue;
1088
1089	smallest_p_lk = BIG;
1090
1091	/* For all the states at node d */
1092	For(i,tree->mod->ns)
1093	{
1094	p1_lk1 = .0;
1095
1096	if(n_v1)
1097	{
1098	/* n_v1 is a tip */
1099	if((n_v1->tax) && (!tree->mod->s_opt->greedy))
1100	{
1101	if(ambiguity_check_v1 == NO)
1102	{
1103	/* For the (non-ambiguous) state at node n_v1 */
1104	p1_lk1 = Pij1[catgNsNs+iNs+state_v1];
1105	}
1106	else
1107	{
1108	/* For all the states at node n_v1 */
1109	For(j,tree->mod->ns)
1110	{
1111	p1_lk1 += Pij1[catgNsNs+iNs+j] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*Ns+j];
1112	}
1113	}
1114	}
1115	/* n_v1 is an internal node */
1116	else
1117	{
1118	/* For the states at node n_v1 */
1119	For(j,tree->mod->ns)
1120	{
1121	p1_lk1 += Pij1[catgNsNs+iNs+j] * p_lk_v1[siteNsNg+catgNs+j];
1122	}
1123	}
1124	}
1125	else
1126	{
1127	p1_lk1 = 1.0;
1128	}
1129
1130	p2_lk2 = .0;
1131
1132	/* We do exactly the same as for node n_v1 but for node n_v2 this time.*/
1133	if(n_v2)
1134	{
1135	/* n_v2 is a tip */
1136	if((n_v2->tax) && (!tree->mod->s_opt->greedy))
1137	{
1138	if(ambiguity_check_v2 == NO)
1139	{
1140	/* For the (non-ambiguous) state at node n_v2 */
1141	p2_lk2 = Pij2[catgNsNs+iNs+state_v2];
1142	}
1143	else
1144	{
1145	/* For all the states at node n_v2 */
1146	For(j,tree->mod->ns)
1147	{
1148	p2_lk2 += Pij2[catgNsNs+iNs+j] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*Ns+j];
1149	}
1150	}
1151	}
1152	/* n_v2 is an internal node */
1153	else
1154	{
1155	/* For all the states at node n_v2 */
1156	For(j,tree->mod->ns)
1157	{
1158	p2_lk2 += Pij2[catgNsNs+iNs+j] * p_lk_v2[siteNsNg+catgNs+j];
1159	}
1160	}
1161	}
1162	else
1163	{
1164	p2_lk2 = 1.0;
1165	}
1166
1167	p_lk[siteNsNg+catgNs+i] = p1_lk1 * p2_lk2;
1168
1169	/* PhyML_Printf("\n+ %G",p_lk[siteNsNg+catgNs+i]); */
1170
1171	if(p_lk[siteNsNg+catgNs+i] < smallest_p_lk) smallest_p_lk = p_lk[siteNsNg+catgNs+i] ;
1172	}
1173
1174	/* Current scaling values at that site */
1175	sum_scale_v1_val = (sum_scale_v1)?(sum_scale_v1[catg*n_patterns+site]):(0);
1176	sum_scale_v2_val = (sum_scale_v2)?(sum_scale_v2[catg*n_patterns+site]):(0);
1177
1178	sum_scale[catg*n_patterns+site] = sum_scale_v1_val + sum_scale_v2_val;
1179
1180	/* PhyML_Printf("\n@ %d",sum_scale[catgn_patterns+site]); /
1181
1182	/* Scaling */
1183	if(smallest_p_lk < p_lk_lim_inf)
1184	{
1185	curr_scaler_pow = (int)(LOG(p_lk_lim_inf)-LOG(smallest_p_lk))/LOG2;
1186	curr_scaler = (phydbl)((unsigned long long)(1) << curr_scaler_pow);
1187
1188	/* if(fabs(curr_scaler_pow) > 63 \|\| fabs(curr_scaler_pow) > 63) */
1189	/* { */
1190	/* PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk); */
1191	/* PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow); */
1192	/* PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__); */
1193	/* Warn_And_Exit("\n"); */
1194	/* } */
1195
1196	sum_scale[catg*n_patterns+site] += curr_scaler_pow;
1197
1198	do
1199	{
1200	piecewise_scaler_pow = MIN(curr_scaler_pow,63);
1201	curr_scaler = (phydbl)((unsigned long long)(1) << piecewise_scaler_pow);
1202	For(i,tree->mod->ns)
1203	{
1204	p_lk[siteNsNg+catgNs+i] *= curr_scaler;
1205
1206	if(p_lk[siteNsNg+catgNs+i] > BIG)
1207	{
1208	PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk);
1209	PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow);
1210	PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__);
1211	Warn_And_Exit("\n");
1212	}
1213	}
1214	curr_scaler_pow -= piecewise_scaler_pow;
1215	}
1216	while(curr_scaler_pow != 0);
1217	}
1218	}
1219	}
1220	}
1221	}
1222
1223
1224	//////////////////////////////////////////////////////////////
1225	//////////////////////////////////////////////////////////////
1226
1227	void Update_P_Lk_Nucl(t_tree tree, t_edge b, t_node *d)
1228	{
1229	/*
1230	\|
1231	\|<- b
1232	\|
1233	d
1234	/ \
1235	/ \
1236	/ \
1237	n_v1 n_v2
1238	*/
1239	t_node n_v1, n_v2;
1240	phydbl p1_lk1,p2_lk2;
1241	phydbl p_lk,p_lk_v1,*p_lk_v2;
1242	phydbl Pij1,Pij2;
1243	int sum_scale, sum_scale_v1, *sum_scale_v2;
1244	int sum_scale_v1_val, sum_scale_v2_val;
1245	int i;
1246	int catg,site;
1247	int n_patterns;
1248	short int ambiguity_check_v1,ambiguity_check_v2;
1249	int state_v1,state_v2;
1250	int dim1, dim2, dim3;
1251	phydbl curr_scaler;
1252	int curr_scaler_pow, piecewise_scaler_pow;
1253	phydbl p_lk_lim_inf;
1254	phydbl smallest_p_lk;
1255	phydbl p0,p1,p2,p3;
1256	int *p_lk_loc;
1257
1258
1259	p_lk_lim_inf = (phydbl)P_LK_LIM_INF;
1260
1261	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
1262	dim2 = tree->mod->ns;
1263	dim3 = tree->mod->ns * tree->mod->ns;
1264
1265	state_v1 = state_v2 = -1;
1266	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1267	sum_scale_v1_val = sum_scale_v2_val = 0;
1268	p1_lk1 = p2_lk2 = .0;
1269
1270	if(d->tax)
1271	{
1272	PhyML_Printf("\n== t_node %d is a leaf...",d->num);
1273	PhyML_Printf("\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
1274	Warn_And_Exit("\n");
1275	}
1276
1277	n_patterns = tree->n_pattern;
1278
1279	n_v1 = n_v2 = NULL;
1280	p_lk = p_lk_v1 = p_lk_v2 = NULL;
1281	Pij1 = Pij2 = NULL;
1282	sum_scale_v1 = sum_scale_v2 = NULL;
1283	p_lk_loc = NULL;
1284
1285	Set_All_P_Lk(&n_v1,&n_v2,
1286	&p_lk,&sum_scale,&p_lk_loc,
1287	&Pij1,&p_lk_v1,&sum_scale_v1,
1288	&Pij2,&p_lk_v2,&sum_scale_v2,
1289	d,b,tree);
1290
1291	/* printf("\n. Tree: %p",(void )tree); /
1292	/* printf("\n. Update on edge %d node %d [%d %d] l: %f",b->num,d->num,n_v1?n_v1->num:-1,n_v2?n_v2->num:-1,b->l->v); */
1293	/* printf("\n. p_lk: %p p_lk_v1: %p p_lk_v2: %p",(void )p_lk,(void )p_lk_v1,(void )p_lk_v2); /
1294	/* printf("\n. Pij1: %p Pij2: %p",(void )Pij1,(void )Pij2); */
1295
1296
1297	/* For every site in the alignment */
1298	For(site,n_patterns)
1299	{
1300	state_v1 = state_v2 = -1;
1301	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1302
1303	if(!tree->mod->s_opt->greedy)
1304	{
1305	/* n_v1 and n_v2 are tip nodes */
1306	if(n_v1 && n_v1->tax)
1307	{
1308	/* Is the state at this tip ambiguous? */
1309	ambiguity_check_v1 = n_v1->c_seq->is_ambigu[site];
1310	if(ambiguity_check_v1 == NO) state_v1 = n_v1->c_seq->d_state[site];
1311	}
1312
1313	if(n_v2 && n_v2->tax)
1314	{
1315	/* Is the state at this tip ambiguous? */
1316	ambiguity_check_v2 = n_v2->c_seq->is_ambigu[site];
1317	if(ambiguity_check_v2 == NO) state_v2 = n_v2->c_seq->d_state[site];
1318	}
1319	}
1320
1321
1322	if(tree->mod->use_m4mod)
1323	{
1324	ambiguity_check_v1 = YES;
1325	ambiguity_check_v2 = YES;
1326	}
1327
1328
1329	if(p_lk_loc[site] < site)
1330	{
1331	Copy_P_Lk(p_lk,p_lk_loc[site],site,tree);
1332	Copy_Scale(sum_scale,p_lk_loc[site],site,tree);
1333	}
1334	else
1335	{
1336	/* For all the rate classes */
1337	For(catg,tree->mod->ras->n_catg)
1338	{
1339	smallest_p_lk = BIG;
1340
1341	/* For all the state at node d */
1342	For(i,tree->mod->ns)
1343	{
1344	p1_lk1 = .0;
1345
1346	if(n_v1)
1347	{
1348	/* n_v1 is a tip */
1349	if((n_v1->tax) && (!tree->mod->s_opt->greedy))
1350	{
1351	if(ambiguity_check_v1 == NO)
1352	{
1353	/* For the (non-ambiguous) state at node n_v1 */
1354	p1_lk1 = Pij1[catgdim3+idim2+state_v1];
1355	}
1356	else
1357	{
1358	/* For all the states at node n_v1 */
1359	p0=Pij1[catgdim3+idim2+0] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+0];
1360	p1=Pij1[catgdim3+idim2+1] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+1];
1361	p2=Pij1[catgdim3+idim2+2] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+2];
1362	p3=Pij1[catgdim3+idim2+3] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+3];
1363	p1_lk1 = p0+p1+p2+p3;
1364
1365
1366	if(isnan(p1_lk1))
1367	{
1368	PhyML_Printf("\n== p0=%f p1=%f p2=%f p3=%f",p0,p1,p2,p3);
1369	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1370	Warn_And_Exit("\n");
1371	}
1372	}
1373	}
1374	/* n_v1 is an internal node */
1375	else
1376	{
1377	/* For the states at node n_v1 */
1378	/* For(j,tree->mod->ns) */
1379	/* { */
1380	/* p1_lk1 += Pij1[catgdim3+idim2+j] * p_lk_v1[sitedim1+catgdim2+j]; */
1381	/* } */
1382
1383	/* printf("\n. PIJ1 %p PLK!:%p %d",Pij1,p_lk_v1,n_v1->num); fflush(NULL); */
1384
1385	p0=Pij1[catgdim3+idim2+0] * p_lk_v1[sitedim1+catgdim2+0];
1386	p1=Pij1[catgdim3+idim2+1] * p_lk_v1[sitedim1+catgdim2+1];
1387	p2=Pij1[catgdim3+idim2+2] * p_lk_v1[sitedim1+catgdim2+2];
1388	p3=Pij1[catgdim3+idim2+3] * p_lk_v1[sitedim1+catgdim2+3];
1389	p1_lk1 = p0+p1+p2+p3;
1390
1391	if(isnan(p1_lk1))
1392	{
1393	PhyML_Printf("\n== p0=%f p1=%f p2=%f p3=%f",p0,p1,p2,p3);
1394	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1395	Warn_And_Exit("\n");
1396	}
1397	}
1398	}
1399	else
1400	{
1401	p1_lk1 = 1.0;
1402	}
1403
1404	p2_lk2 = .0;
1405
1406	/* We do exactly the same as for node n_v1 but for node n_v2 this time.*/
1407	if(n_v2)
1408	{
1409	/* n_v2 is a tip */
1410	if((n_v2->tax) && (!tree->mod->s_opt->greedy))
1411	{
1412	if(ambiguity_check_v2 == NO)
1413	{
1414	/* For the (non-ambiguous) state at node n_v2 */
1415	p2_lk2 = Pij2[catgdim3+idim2+state_v2];
1416	if(isnan(p2_lk2))
1417	{
1418	PhyML_Printf("\n== Tree %d",tree->tree_num);
1419	PhyML_Printf("\n== catg=%d dim3=%d dim2=%d i=%d state_v2=%d",catg,dim3,dim2,i,state_v2);
1420	PhyML_Printf("\n== Pij2[0] = %f",Pij2[0]);
1421	PhyML_Printf("\n== q[0]=%f",tree->mod->eigen->q[0]);
1422	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1423	Warn_And_Exit("\n");
1424	}
1425	}
1426	else
1427	{
1428	/* For all the states at node n_v2 */
1429	p0=Pij2[catgdim3+idim2+0] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+0];
1430	p1=Pij2[catgdim3+idim2+1] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+1];
1431	p2=Pij2[catgdim3+idim2+2] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+2];
1432	p3=Pij2[catgdim3+idim2+3] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+3];
1433	p2_lk2 = p0+p1+p2+p3;
1434
1435
1436	if(isnan(p2_lk2))
1437	{
1438	PhyML_Printf("\n== p0=%f p1=%f p2=%f p3=%f",p0,p1,p2,p3);
1439	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1440	Warn_And_Exit("\n");
1441	}
1442	}
1443	}
1444	/* n_v2 is an internal node */
1445	else
1446	{
1447	/* For all the states at node n_v2 */
1448	p0=Pij2[catgdim3+idim2+0] * p_lk_v2[sitedim1+catgdim2+0];
1449	p1=Pij2[catgdim3+idim2+1] * p_lk_v2[sitedim1+catgdim2+1];
1450	p2=Pij2[catgdim3+idim2+2] * p_lk_v2[sitedim1+catgdim2+2];
1451	p3=Pij2[catgdim3+idim2+3] * p_lk_v2[sitedim1+catgdim2+3];
1452	p2_lk2 = p0+p1+p2+p3;
1453	if(isnan(p2_lk2))
1454	{
1455	PhyML_Printf("\n== %f %f",b->l->v,b->l->vb->l->vtree->mod->l_var_sigma);
1456	PhyML_Printf("\n== p0=%f p1=%f p2=%f p3=%f",p0,p1,p2,p3);
1457	PhyML_Printf("\n== Pij2[0]=%f Pij2[1]=%f Pij2[2]=%f Pij2[3]=%f",Pij2[catgdim3+idim2+0],Pij2[catgdim3+idim2+1],Pij2[catgdim3+idim2+2],Pij2[catgdim3+idim2+3]);
1458	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1459	Warn_And_Exit("\n");
1460	}
1461	}
1462	}
1463	else
1464	{
1465	p2_lk2 = 1.0;
1466	}
1467
1468	p_lk[sitedim1+catgdim2+i] = p1_lk1 * p2_lk2;
1469
1470	/* printf("\n.<< site %3d catg: %3d i: %3d p_lk %f p1_lk1: %f p2_lk2: %f", */
1471	/* site, */
1472	/* catg, */
1473	/* i, */
1474	/* p_lk[sitedim1+catgdim2+i], */
1475	/* p1_lk1, */
1476	/* p2_lk2); */
1477
1478	if(isnan(p2_lk2))
1479	{
1480	PhyML_Printf("\n== Err. in file %s at line %d.",__FILE__,__LINE__);
1481	Warn_And_Exit("\n");
1482	}
1483
1484	if(p_lk[sitedim1+catgdim2+i] < smallest_p_lk) smallest_p_lk = p_lk[sitedim1+catgdim2+i] ;
1485	}
1486
1487	/* if(n_v1->tax == NO && n_v2->tax == NO) */
1488	/* printf("\n<< site %3d node %3d %3d %3d %f %f %f %f v1 : %f %f %f %f v2: %f %f %f %f>>", */
1489	/* site,d->num,n_v1->num,n_v2->num,p_lk[0],p_lk[1],p_lk[2],p_lk[3], */
1490	/* (phydbl)p_lk_v1[sitedim1+catgdim2+0], */
1491	/* (phydbl)p_lk_v1[sitedim1+catgdim2+1], */
1492	/* (phydbl)p_lk_v1[sitedim1+catgdim2+2], */
1493	/* (phydbl)p_lk_v1[sitedim1+catgdim2+3], */
1494	/* (phydbl)p_lk_v2[sitedim1+catgdim2+0], */
1495	/* (phydbl)p_lk_v2[sitedim1+catgdim2+1], */
1496	/* (phydbl)p_lk_v2[sitedim1+catgdim2+2], */
1497	/* (phydbl)p_lk_v2[sitedim1+catgdim2+3]); */
1498	/* else if(n_v1->tax == NO && n_v2->tax == YES) */
1499	/* printf("\n<< site %3d node %3d %3d %3d %f %f %f %f v1 : %f %f %f %f v2: %f %f %f %f>>", /
1500	/* site,d->num,n_v1->num,n_v2->num,p_lk[0],p_lk[1],p_lk[2],p_lk[3], */
1501	/* (phydbl)p_lk_v1[sitedim1+catgdim2+0], */
1502	/* (phydbl)p_lk_v1[sitedim1+catgdim2+1], */
1503	/* (phydbl)p_lk_v1[sitedim1+catgdim2+2], */
1504	/* (phydbl)p_lk_v1[sitedim1+catgdim2+3], */
1505	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+0], /
1506	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+1], /
1507	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+2], /
1508	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+3]); /
1509	/* else if(n_v1->tax == YES && n_v2->tax == NO) */
1510	/* printf("\n<< site %3d node %3d %3d %3d %f %f %f %f v1: %f %f %f %f v2 : %f %f %f %f>>", /
1511	/* site,d->num,n_v1->num,n_v2->num,p_lk[0],p_lk[1],p_lk[2],p_lk[3], */
1512	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+0], /
1513	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+1], /
1514	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+2], /
1515	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+3], /
1516	/* (phydbl)p_lk_v2[sitedim1+catgdim2+0], */
1517	/* (phydbl)p_lk_v2[sitedim1+catgdim2+1], */
1518	/* (phydbl)p_lk_v2[sitedim1+catgdim2+2], */
1519	/* (phydbl)p_lk_v2[sitedim1+catgdim2+3]); */
1520	/* else if(n_v1->tax == YES && n_v2->tax == YES) */
1521	/* printf("\n<< site %3d node %3d %3d %3d %f %f %f %f v1: %f %f %f %f v2: %f %f %f %f>>", */
1522	/* site,d->num,n_v1->num,n_v2->num,p_lk[0],p_lk[1],p_lk[2],p_lk[3], */
1523	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+0], /
1524	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+1], /
1525	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+2], /
1526	/* (phydbl)n_v1->b[0]->p_lk_tip_r[sitedim2+3], /
1527	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+0], /
1528	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+1], /
1529	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+2], /
1530	/* (phydbl)n_v2->b[0]->p_lk_tip_r[sitedim2+3]); /
1531	/* Print_Square_Matrix_Generic(4,Pij1); */
1532	/* Print_Square_Matrix_Generic(4,Pij2); */
1533
1534
1535	/* Current scaling values at that site */
1536	sum_scale_v1_val = (sum_scale_v1)?(sum_scale_v1[catg*n_patterns+site]):(0);
1537	sum_scale_v2_val = (sum_scale_v2)?(sum_scale_v2[catg*n_patterns+site]):(0);
1538
1539	sum_scale[catg*n_patterns+site] = sum_scale_v1_val + sum_scale_v2_val;
1540
1541	/* Scaling */
1542	if(smallest_p_lk < p_lk_lim_inf)
1543	{
1544	curr_scaler_pow = (int)(LOG(p_lk_lim_inf)-LOG(smallest_p_lk))/LOG2;
1545	curr_scaler = (phydbl)((unsigned long long)(1) << curr_scaler_pow);
1546
1547	/* if(fabs(curr_scaler_pow) > 63 \|\| fabs(curr_scaler_pow) > 63) */
1548	/* { */
1549	/* PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk); */
1550	/* PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow); */
1551	/* PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__); */
1552	/* Warn_And_Exit("\n"); */
1553	/* } */
1554
1555	sum_scale[catg*n_patterns+site] += curr_scaler_pow;
1556
1557	do
1558	{
1559	piecewise_scaler_pow = MIN(curr_scaler_pow,63);
1560	curr_scaler = (phydbl)((unsigned long long)(1) << piecewise_scaler_pow);
1561	For(i,tree->mod->ns)
1562	{
1563	p_lk[sitedim1+catgdim2+i] *= curr_scaler;
1564
1565	if(p_lk[sitedim1+catgdim2+i] > BIG)
1566	{
1567	PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk);
1568	PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow);
1569	PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__);
1570	Warn_And_Exit("\n");
1571	}
1572	}
1573	curr_scaler_pow -= piecewise_scaler_pow;
1574	}
1575	while(curr_scaler_pow != 0);
1576	}
1577	}
1578	}
1579	}
1580	}
1581
1582	//////////////////////////////////////////////////////////////
1583	//////////////////////////////////////////////////////////////
1584
1585	void Update_P_Lk_AA(t_tree tree, t_edge b, t_node *d)
1586	{
1587	/*
1588	\|
1589	\|<- b
1590	\|
1591	d
1592	/ \
1593	/ \
1594	/ \
1595	n_v1 n_v2
1596	*/
1597	t_node n_v1, n_v2;
1598	phydbl p1_lk1,p2_lk2;
1599	phydbl p_lk,p_lk_v1,*p_lk_v2;
1600	phydbl Pij1,Pij2;
1601	int sum_scale, sum_scale_v1, *sum_scale_v2;
1602	int sum_scale_v1_val, sum_scale_v2_val;
1603	int i;
1604	int catg,site;
1605	int n_patterns;
1606	short int ambiguity_check_v1,ambiguity_check_v2;
1607	int state_v1,state_v2;
1608	int dim1, dim2, dim3;
1609	phydbl curr_scaler;
1610	int curr_scaler_pow, piecewise_scaler_pow;
1611	phydbl p_lk_lim_inf;
1612	phydbl smallest_p_lk;
1613	phydbl p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19;
1614	int *p_lk_loc;
1615
1616	p_lk_lim_inf = (phydbl)P_LK_LIM_INF;
1617
1618	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
1619	dim2 = tree->mod->ns;
1620	dim3 = tree->mod->ns * tree->mod->ns;
1621
1622	state_v1 = state_v2 = -1;
1623	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1624	sum_scale_v1_val = sum_scale_v2_val = 0;
1625	p1_lk1 = p2_lk2 = .0;
1626
1627	if(d->tax)
1628	{
1629	PhyML_Printf("\n== t_node %d is a leaf...",d->num);
1630	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
1631	Warn_And_Exit("\n");
1632	}
1633
1634	n_patterns = tree->n_pattern;
1635
1636	n_v1 = n_v2 = NULL;
1637	p_lk = p_lk_v1 = p_lk_v2 = NULL;
1638	Pij1 = Pij2 = NULL;
1639	sum_scale_v1 = sum_scale_v2 = NULL;
1640	p_lk_loc = NULL;
1641
1642	Set_All_P_Lk(&n_v1,&n_v2,
1643	&p_lk,&sum_scale,&p_lk_loc,
1644	&Pij1,&p_lk_v1,&sum_scale_v1,
1645	&Pij2,&p_lk_v2,&sum_scale_v2,
1646	d,b,tree);
1647
1648
1649	/* For every site in the alignment */
1650	For(site,n_patterns)
1651	{
1652	state_v1 = state_v2 = -1;
1653	ambiguity_check_v1 = ambiguity_check_v2 = NO;
1654
1655	if(!tree->mod->s_opt->greedy)
1656	{
1657	/* n_v1 and n_v2 are tip nodes */
1658	if(n_v1 && n_v1->tax)
1659	{
1660	/* Is the state at this tip ambiguous? */
1661	ambiguity_check_v1 = n_v1->c_seq->is_ambigu[site];
1662	/* if(ambiguity_check_v1 == NO) state_v1 = Get_State_From_P_Pars(n_v1->b[0]->p_lk_tip_r,sitedim2,tree); /
1663	if(ambiguity_check_v1 == NO) state_v1 = n_v1->c_seq->d_state[site];
1664	}
1665
1666	if(n_v2 && n_v2->tax)
1667	{
1668	/* Is the state at this tip ambiguous? */
1669	ambiguity_check_v2 = n_v2->c_seq->is_ambigu[site];
1670	/* if(ambiguity_check_v2 == NO) state_v2 = Get_State_From_P_Pars(n_v2->b[0]->p_lk_tip_r,sitedim2,tree); /
1671	if(ambiguity_check_v2 == NO) state_v2 = n_v2->c_seq->d_state[site];
1672	}
1673	}
1674
1675	if(tree->mod->use_m4mod)
1676	{
1677	ambiguity_check_v1 = YES;
1678	ambiguity_check_v2 = YES;
1679	}
1680
1681	if(p_lk_loc[site] < site)
1682	{
1683	Copy_P_Lk(p_lk,p_lk_loc[site],site,tree);
1684	Copy_Scale(sum_scale,p_lk_loc[site],site,tree);
1685	}
1686	else
1687	{
1688	/* For all the rate classes */
1689	For(catg,tree->mod->ras->n_catg)
1690	{
1691	smallest_p_lk = BIG;
1692
1693	/* For all the state at node d */
1694	For(i,tree->mod->ns)
1695	{
1696	p1_lk1 = .0;
1697
1698	if(n_v1)
1699	{
1700	/* n_v1 is a tip */
1701	if((n_v1->tax) && (!tree->mod->s_opt->greedy))
1702	{
1703	if(ambiguity_check_v1 == NO)
1704	{
1705	/* For the (non-ambiguous) state at node n_v1 */
1706	p1_lk1 = Pij1[catgdim3+idim2+state_v1];
1707	}
1708	else
1709	{
1710	/* For all the states at node n_v1 */
1711	p0 = Pij1[catgdim3+idim2+ 0] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 0];
1712	p1 = Pij1[catgdim3+idim2+ 1] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 1];
1713	p2 = Pij1[catgdim3+idim2+ 2] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 2];
1714	p3 = Pij1[catgdim3+idim2+ 3] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 3];
1715	p4 = Pij1[catgdim3+idim2+ 4] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 4];
1716	p5 = Pij1[catgdim3+idim2+ 5] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 5];
1717	p6 = Pij1[catgdim3+idim2+ 6] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 6];
1718	p7 = Pij1[catgdim3+idim2+ 7] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 7];
1719	p8 = Pij1[catgdim3+idim2+ 8] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 8];
1720	p9 = Pij1[catgdim3+idim2+ 9] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+ 9];
1721	p10 = Pij1[catgdim3+idim2+10] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+10];
1722	p11 = Pij1[catgdim3+idim2+11] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+11];
1723	p12 = Pij1[catgdim3+idim2+12] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+12];
1724	p13 = Pij1[catgdim3+idim2+13] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+13];
1725	p14 = Pij1[catgdim3+idim2+14] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+14];
1726	p15 = Pij1[catgdim3+idim2+15] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+15];
1727	p16 = Pij1[catgdim3+idim2+16] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+16];
1728	p17 = Pij1[catgdim3+idim2+17] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+17];
1729	p18 = Pij1[catgdim3+idim2+18] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+18];
1730	p19 = Pij1[catgdim3+idim2+19] * (phydbl)n_v1->b[0]->p_lk_tip_r[site*dim2+19];
1731	p1_lk1 = p0+p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12+p13+p14+p15+p16+p17+p18+p19;
1732	}
1733	}
1734	/* n_v1 is an internal node */
1735	else
1736	{
1737	/* For the states at node n_v1 */
1738	p0 = Pij1[catgdim3+idim2+ 0] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 0];
1739	p1 = Pij1[catgdim3+idim2+ 1] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 1];
1740	p2 = Pij1[catgdim3+idim2+ 2] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 2];
1741	p3 = Pij1[catgdim3+idim2+ 3] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 3];
1742	p4 = Pij1[catgdim3+idim2+ 4] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 4];
1743	p5 = Pij1[catgdim3+idim2+ 5] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 5];
1744	p6 = Pij1[catgdim3+idim2+ 6] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 6];
1745	p7 = Pij1[catgdim3+idim2+ 7] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 7];
1746	p8 = Pij1[catgdim3+idim2+ 8] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 8];
1747	p9 = Pij1[catgdim3+idim2+ 9] * (phydbl)p_lk_v1[sitedim1+catgdim2+ 9];
1748	p10 = Pij1[catgdim3+idim2+10] * (phydbl)p_lk_v1[sitedim1+catgdim2+10];
1749	p11 = Pij1[catgdim3+idim2+11] * (phydbl)p_lk_v1[sitedim1+catgdim2+11];
1750	p12 = Pij1[catgdim3+idim2+12] * (phydbl)p_lk_v1[sitedim1+catgdim2+12];
1751	p13 = Pij1[catgdim3+idim2+13] * (phydbl)p_lk_v1[sitedim1+catgdim2+13];
1752	p14 = Pij1[catgdim3+idim2+14] * (phydbl)p_lk_v1[sitedim1+catgdim2+14];
1753	p15 = Pij1[catgdim3+idim2+15] * (phydbl)p_lk_v1[sitedim1+catgdim2+15];
1754	p16 = Pij1[catgdim3+idim2+16] * (phydbl)p_lk_v1[sitedim1+catgdim2+16];
1755	p17 = Pij1[catgdim3+idim2+17] * (phydbl)p_lk_v1[sitedim1+catgdim2+17];
1756	p18 = Pij1[catgdim3+idim2+18] * (phydbl)p_lk_v1[sitedim1+catgdim2+18];
1757	p19 = Pij1[catgdim3+idim2+19] * (phydbl)p_lk_v1[sitedim1+catgdim2+19];
1758	p1_lk1 = p0+p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12+p13+p14+p15+p16+p17+p18+p19;
1759	}
1760	}
1761	else
1762	{
1763	p1_lk1 = 1.0;
1764	}
1765
1766	p2_lk2 = .0;
1767
1768	/* We do exactly the same as for node n_v1 but for node n_v2 this time.*/
1769	if(n_v2)
1770	{
1771	/* n_v2 is a tip */
1772	if((n_v2->tax) && (!tree->mod->s_opt->greedy))
1773	{
1774	if(ambiguity_check_v2 == NO)
1775	{
1776	/* For the (non-ambiguous) state at node n_v2 */
1777	p2_lk2 = Pij2[catgdim3+idim2+state_v2];
1778	}
1779	else
1780	{
1781	/* For all the states at node n_v2 */
1782	p0 = Pij2[catgdim3+idim2+ 0] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 0];
1783	p1 = Pij2[catgdim3+idim2+ 1] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 1];
1784	p2 = Pij2[catgdim3+idim2+ 2] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 2];
1785	p3 = Pij2[catgdim3+idim2+ 3] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 3];
1786	p4 = Pij2[catgdim3+idim2+ 4] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 4];
1787	p5 = Pij2[catgdim3+idim2+ 5] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 5];
1788	p6 = Pij2[catgdim3+idim2+ 6] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 6];
1789	p7 = Pij2[catgdim3+idim2+ 7] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 7];
1790	p8 = Pij2[catgdim3+idim2+ 8] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 8];
1791	p9 = Pij2[catgdim3+idim2+ 9] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+ 9];
1792	p10 = Pij2[catgdim3+idim2+10] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+10];
1793	p11 = Pij2[catgdim3+idim2+11] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+11];
1794	p12 = Pij2[catgdim3+idim2+12] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+12];
1795	p13 = Pij2[catgdim3+idim2+13] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+13];
1796	p14 = Pij2[catgdim3+idim2+14] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+14];
1797	p15 = Pij2[catgdim3+idim2+15] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+15];
1798	p16 = Pij2[catgdim3+idim2+16] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+16];
1799	p17 = Pij2[catgdim3+idim2+17] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+17];
1800	p18 = Pij2[catgdim3+idim2+18] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+18];
1801	p19 = Pij2[catgdim3+idim2+19] * (phydbl)n_v2->b[0]->p_lk_tip_r[site*dim2+19];
1802	p2_lk2 = p0+p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12+p13+p14+p15+p16+p17+p18+p19;
1803
1804	}
1805	}
1806	/* n_v2 is an internal node */
1807	else
1808	{
1809	/* For all the states at node n_v2 */
1810	p0 = Pij2[catgdim3+idim2+ 0] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 0];
1811	p1 = Pij2[catgdim3+idim2+ 1] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 1];
1812	p2 = Pij2[catgdim3+idim2+ 2] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 2];
1813	p3 = Pij2[catgdim3+idim2+ 3] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 3];
1814	p4 = Pij2[catgdim3+idim2+ 4] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 4];
1815	p5 = Pij2[catgdim3+idim2+ 5] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 5];
1816	p6 = Pij2[catgdim3+idim2+ 6] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 6];
1817	p7 = Pij2[catgdim3+idim2+ 7] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 7];
1818	p8 = Pij2[catgdim3+idim2+ 8] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 8];
1819	p9 = Pij2[catgdim3+idim2+ 9] * (phydbl)p_lk_v2[sitedim1+catgdim2+ 9];
1820	p10 = Pij2[catgdim3+idim2+10] * (phydbl)p_lk_v2[sitedim1+catgdim2+10];
1821	p11 = Pij2[catgdim3+idim2+11] * (phydbl)p_lk_v2[sitedim1+catgdim2+11];
1822	p12 = Pij2[catgdim3+idim2+12] * (phydbl)p_lk_v2[sitedim1+catgdim2+12];
1823	p13 = Pij2[catgdim3+idim2+13] * (phydbl)p_lk_v2[sitedim1+catgdim2+13];
1824	p14 = Pij2[catgdim3+idim2+14] * (phydbl)p_lk_v2[sitedim1+catgdim2+14];
1825	p15 = Pij2[catgdim3+idim2+15] * (phydbl)p_lk_v2[sitedim1+catgdim2+15];
1826	p16 = Pij2[catgdim3+idim2+16] * (phydbl)p_lk_v2[sitedim1+catgdim2+16];
1827	p17 = Pij2[catgdim3+idim2+17] * (phydbl)p_lk_v2[sitedim1+catgdim2+17];
1828	p18 = Pij2[catgdim3+idim2+18] * (phydbl)p_lk_v2[sitedim1+catgdim2+18];
1829	p19 = Pij2[catgdim3+idim2+19] * (phydbl)p_lk_v2[sitedim1+catgdim2+19];
1830	p2_lk2 = p0+p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12+p13+p14+p15+p16+p17+p18+p19;
1831	}
1832	}
1833	else
1834	{
1835	p2_lk2 = 1.0;
1836	}
1837
1838	p_lk[sitedim1+catgdim2+i] = p1_lk1 * p2_lk2;
1839
1840	if(p_lk[sitedim1+catgdim2+i] < smallest_p_lk) smallest_p_lk = p_lk[sitedim1+catgdim2+i] ;
1841	}
1842
1843	/* Current scaling values at that site */
1844	sum_scale_v1_val = (sum_scale_v1)?(sum_scale_v1[catg*n_patterns+site]):(0);
1845	sum_scale_v2_val = (sum_scale_v2)?(sum_scale_v2[catg*n_patterns+site]):(0);
1846
1847	sum_scale[catg*n_patterns+site] = sum_scale_v1_val + sum_scale_v2_val;
1848
1849	/* Scaling */
1850	if(smallest_p_lk < p_lk_lim_inf)
1851	{
1852	curr_scaler_pow = (int)(LOG(p_lk_lim_inf)-LOG(smallest_p_lk))/LOG2;
1853	curr_scaler = (phydbl)((unsigned long long)(1) << curr_scaler_pow);
1854
1855	/* if(fabs(curr_scaler_pow) > 63 \|\| fabs(curr_scaler_pow) > 63) */
1856	/* { */
1857	/* PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk); */
1858	/* PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow); */
1859	/* PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__); */
1860	/* Warn_And_Exit("\n"); */
1861	/* } */
1862
1863	sum_scale[catg*n_patterns+site] += curr_scaler_pow;
1864
1865	do
1866	{
1867	piecewise_scaler_pow = MIN(curr_scaler_pow,63);
1868	curr_scaler = (phydbl)((unsigned long long)(1) << piecewise_scaler_pow);
1869	For(i,tree->mod->ns)
1870	{
1871	p_lk[sitedim1+catgdim2+i] *= curr_scaler;
1872
1873	if(p_lk[sitedim1+catgdim2+i] > BIG)
1874	{
1875	PhyML_Printf("\n. p_lk_lim_inf = %G smallest_p_lk = %G",p_lk_lim_inf,smallest_p_lk);
1876	PhyML_Printf("\n. curr_scaler_pow = %d",curr_scaler_pow);
1877	PhyML_Printf("\n. Err in file %s at line %d.",__FILE__,__LINE__);
1878	Warn_And_Exit("\n");
1879	}
1880	}
1881	curr_scaler_pow -= piecewise_scaler_pow;
1882	}
1883	while(curr_scaler_pow != 0);
1884	}
1885	}
1886	}
1887	}
1888	}
1889
1890
1891	//////////////////////////////////////////////////////////////
1892	//////////////////////////////////////////////////////////////
1893
1894
1895	phydbl Return_Abs_Lk(t_tree *tree)
1896	{
1897	Lk(NULL,tree);
1898	return FABS(tree->c_lnL);
1899	}
1900
1901	//////////////////////////////////////////////////////////////
1902	//////////////////////////////////////////////////////////////
1903
1904	matrix ML_Dist(calign data, t_mod *mod)
1905	{
1906	int i,j,k,l;
1907	phydbl init;
1908	int n_catg;
1909	phydbl d_max,sum;
1910	matrix *mat;
1911	calign twodata,tmpdata;
1912	int state0, state1,len;
1913	phydbl *F;
1914	eigen *eigen_struct;
1915
1916	tmpdata = (calign *)mCalloc(1,sizeof(calign));
1917	tmpdata->c_seq = (align *)mCalloc(2,sizeof(align ));
1918	tmpdata->b_frq = (phydbl *)mCalloc(mod->ns,sizeof(phydbl));
1919	tmpdata->ambigu = (short int *)mCalloc(data->crunch_len,sizeof(short int));
1920	F = (phydbl )mCalloc(mod->nsmod->ns,sizeof(phydbl ));
1921	eigen_struct = (eigen *)Make_Eigen_Struct(mod->ns);
1922
1923	tmpdata->n_otu = 2;
1924
1925	tmpdata->crunch_len = data->crunch_len;
1926	tmpdata->init_len = data->init_len;
1927
1928	mat = NULL;
1929	if(mod->io->datatype == NT) mat = (mod->whichmodel < 10)?(K80_dist(data,1E+6)):(JC69_Dist(data,mod));
1930	else if(mod->io->datatype == AA) mat = JC69_Dist(data,mod);
1931	else if(mod->io->datatype == GENERIC) mat = JC69_Dist(data,mod);
1932
1933	For(i,mod->ras->n_catg) /* Don't use the discrete gamma distribution */
1934	{
1935	mod->ras->gamma_rr->v[i] = 1.0;
1936	mod->ras->gamma_r_proba->v[i] = 1.0;
1937	}
1938
1939	n_catg = mod->ras->n_catg;
1940	mod->ras->n_catg = 1;
1941
1942	For(j,data->n_otu-1)
1943	{
1944	tmpdata->c_seq[0] = data->c_seq[j];
1945	tmpdata->c_seq[0]->name = data->c_seq[j]->name;
1946	tmpdata->wght = data->wght;
1947
1948
1949	for(k=j+1;k<data->n_otu;k++)
1950	{
1951
1952	tmpdata->c_seq[1] = data->c_seq[k];
1953	tmpdata->c_seq[1]->name = data->c_seq[k]->name;
1954
1955	twodata = Compact_Cdata(tmpdata,mod->io);
1956
1957	Check_Ambiguities(twodata,mod->io->datatype,mod->io->state_len);
1958
1959	Hide_Ambiguities(twodata);
1960
1961	init = mat->dist[j][k];
1962
1963	if((init > DIST_MAX-SMALL) \|\| (init < .0)) init = 0.1;
1964
1965	d_max = init;
1966
1967	For(i,mod->ns*mod->ns) F[i]=.0;
1968	len = 0;
1969	For(l,twodata->c_seq[0]->len)
1970	{
1971	state0 = Assign_State(twodata->c_seq[0]->state+l*mod->io->state_len,mod->io->datatype,mod->io->state_len);
1972	state1 = Assign_State(twodata->c_seq[1]->state+l*mod->io->state_len,mod->io->datatype,mod->io->state_len);
1973
1974	if((state0 > -1) && (state1 > -1))
1975	{
1976	F[mod->ns*state0+state1] += twodata->wght[l];
1977	len += (int)twodata->wght[l];
1978	}
1979	}
1980	if(len > .0) {For(i,mod->ns*mod->ns) F[i] /= (phydbl)len;}
1981
1982	sum = 0.;
1983	For(i,mod->ns*mod->ns) sum += F[i];
1984
1985
1986	/* if(sum < .001) d_max = -1.; */
1987	if(sum < .001) d_max = init;
1988	else if((sum > 1. - .001) && (sum < 1. + .001)) Opt_Dist_F(&(d_max),F,mod);
1989	else
1990	{
1991	PhyML_Printf("\n== sum = %f\n",sum);
1992	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
1993	Exit("");
1994	}
1995
1996	if(d_max >= DIST_MAX) d_max = DIST_MAX;
1997
1998
1999	/* Do not correct for dist < BL_MIN, otherwise Fill_Missing_Dist
2000	* will not be called
2001	*/
2002	mat->dist[j][k] = d_max;
2003	mat->dist[k][j] = mat->dist[j][k];
2004	Free_Cseq(twodata);
2005	}
2006	}
2007
2008	mod->ras->n_catg = n_catg;
2009
2010
2011	Free(tmpdata->ambigu);
2012	Free(tmpdata->b_frq);
2013	Free(tmpdata->c_seq);
2014	free(tmpdata);
2015	Free_Eigen(eigen_struct);
2016	Free(F);
2017
2018	return mat;
2019	}
2020
2021	//////////////////////////////////////////////////////////////
2022	//////////////////////////////////////////////////////////////
2023
2024
2025	phydbl Lk_Given_Two_Seq(calign data, int numseq1, int numseq2, phydbl dist, t_mod mod, phydbl *loglk)
2026	{
2027	align seq1,seq2;
2028	phydbl site_lk,log_site_lk;
2029	int i,j,k,l;
2030	/* phydbl p_lk_l,p_lk_r; */
2031	phydbl p_lk_l,p_lk_r;
2032	phydbl len;
2033	int dim1,dim2;
2034
2035	dim1 = mod->ns;
2036	dim2 = mod->ns * mod->ns;
2037
2038	DiscreteGamma(mod->ras->gamma_r_proba->v, mod->ras->gamma_rr->v, mod->ras->alpha->v,
2039	mod->ras->alpha->v,mod->ras->n_catg,mod->ras->gamma_median);
2040
2041	seq1 = data->c_seq[numseq1];
2042	seq2 = data->c_seq[numseq2];
2043
2044
2045	p_lk_l = (phydbl )mCalloc(data->c_seq[0]->len mod->ns,sizeof(phydbl));
2046	p_lk_r = (phydbl )mCalloc(data->c_seq[0]->len mod->ns,sizeof(phydbl));
2047
2048
2049	For(i,mod->ras->n_catg)
2050	{
2051	len = dist*mod->ras->gamma_rr->v[i];
2052	if(len < mod->l_min) len = mod->l_min;
2053	else if(len > mod->l_max) len = mod->l_max;
2054	PMat(len,mod,dim2*i,mod->Pij_rr->v);
2055	}
2056
2057	if(mod->io->datatype == NT)
2058	{
2059	For(i,data->c_seq[0]->len)
2060	{
2061	Init_Tips_At_One_Site_Nucleotides_Float(seq1->state[i],i*mod->ns,p_lk_l);
2062	Init_Tips_At_One_Site_Nucleotides_Float(seq2->state[i],i*mod->ns,p_lk_r);
2063	}
2064	}
2065	else if(mod->io->datatype == AA)
2066	{
2067	For(i,data->c_seq[0]->len)
2068	{
2069	Init_Tips_At_One_Site_AA_Float(seq1->state[i],i*mod->ns,p_lk_l);
2070	Init_Tips_At_One_Site_AA_Float(seq2->state[i],i*mod->ns,p_lk_r);
2071	}
2072	}
2073	else
2074	{
2075	PhyML_Printf("\n. Not implemented yet...");
2076	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2077	Warn_And_Exit("\n");
2078	}
2079
2080
2081	site_lk = .0;
2082	*loglk = 0;
2083
2084	For(i,data->c_seq[0]->len)
2085	{
2086	if(data->wght[i])
2087	{
2088	site_lk = log_site_lk = .0;
2089	if(!data->ambigu[i])
2090	{
2091	For(k,mod->ns) {if(p_lk_l[i*mod->ns+k] > .0001) break;}
2092	For(l,mod->ns) {if(p_lk_r[i*mod->ns+l] > .0001) break;}
2093	For(j,mod->ras->n_catg)
2094	{
2095	site_lk +=
2096	mod->ras->gamma_r_proba->v[j] *
2097	mod->e_frq->pi->v[k] *
2098	p_lk_l[idim1+k]
2099	mod->Pij_rr->v[jdim2+kdim1+l] *
2100	p_lk_r[i*dim1+l];
2101	}
2102	}
2103	else
2104	{
2105	For(j,mod->ras->n_catg)
2106	{
2107	For(k,mod->ns) /sort sum terms ? No global effect/
2108	{
2109	For(l,mod->ns)
2110	{
2111	site_lk +=
2112	mod->ras->gamma_r_proba->v[j] *
2113	mod->e_frq->pi->v[k] *
2114	p_lk_l[idim1+k]
2115	mod->Pij_rr->v[jdim2+kdim1+l] *
2116	p_lk_r[i*dim1+l];
2117	}
2118	}
2119	}
2120	}
2121
2122	if(site_lk <= .0)
2123	{
2124	PhyML_Printf("'%c' '%c'\n",seq1->state[i],seq2->state[i]);
2125	Exit("\n. Err: site lk <= 0\n");
2126	}
2127
2128	log_site_lk += (phydbl)LOG(site_lk);
2129
2130	loglk += data->wght[i] log_site_lk;/* sort sum terms ? No global effect*/
2131	}
2132	}
2133
2134	/* For(i,data->c_seq[0]->len) */
2135	/* { */
2136	/* Free(p_lk_l[i]); */
2137	/* Free(p_lk_r[i]); */
2138	/* } */
2139
2140	Free(p_lk_l); Free(p_lk_r);
2141	return *loglk;
2142	}
2143
2144	//////////////////////////////////////////////////////////////
2145	//////////////////////////////////////////////////////////////
2146
2147
2148	void Unconstraint_Lk(t_tree *tree)
2149	{
2150	int i;
2151
2152	tree->unconstraint_lk = .0;
2153
2154	For(i,tree->data->crunch_len)
2155	{
2156	tree->unconstraint_lk +=
2157	tree->data->wght[i]*(phydbl)LOG(tree->data->wght[i]);
2158	}
2159	tree->unconstraint_lk -=
2160	tree->data->init_len*(phydbl)LOG(tree->data->init_len);
2161	}
2162
2163	//////////////////////////////////////////////////////////////
2164	//////////////////////////////////////////////////////////////
2165
2166	void Init_P_Lk_Tips_Double(t_tree *tree)
2167	{
2168	int curr_site,i,j,k,dim1,dim2;
2169
2170	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
2171	dim2 = tree->mod->ns;
2172
2173
2174	For(i,tree->n_otu)
2175	{
2176	if(!tree->a_nodes[i]->c_seq \|\|
2177	strcmp(tree->a_nodes[i]->c_seq->name,tree->a_nodes[i]->name))
2178	{
2179	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
2180	Exit("");
2181	}
2182	}
2183
2184	For(curr_site,tree->data->crunch_len)
2185	{
2186	For(i,tree->n_otu)
2187	{
2188	if (tree->io->datatype == NT)
2189	/* Init_Tips_At_One_Site_Nucleotides_Float(tree->data->c_seq[i]->state[curr_site], */
2190	/* curr_sitedim1+0dim2, */
2191	/* tree->a_nodes[i]->b[0]->p_lk_rght); */
2192	Init_Tips_At_One_Site_Nucleotides_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2193	curr_sitedim1+0dim2,
2194	tree->a_nodes[i]->b[0]->p_lk_rght);
2195
2196	else if(tree->io->datatype == AA)
2197	Init_Tips_At_One_Site_AA_Float(tree->a_nodes[i]->c_seq->state[curr_site],
2198	curr_sitedim1+0dim2,
2199	tree->a_nodes[i]->b[0]->p_lk_rght);
2200	/* Init_Tips_At_One_Site_AA_Float(tree->data->c_seq[i]->state[curr_site], */
2201	/* curr_sitedim1+0dim2, */
2202	/* tree->a_nodes[i]->b[0]->p_lk_rght); */
2203
2204	else if(tree->io->datatype == GENERIC)
2205	Init_Tips_At_One_Site_Generic_Float(tree->a_nodes[i]->c_seq->state+curr_site*tree->mod->io->state_len,
2206	tree->mod->ns,
2207	tree->mod->io->state_len,
2208	curr_sitedim1+0dim2,
2209	tree->a_nodes[i]->b[0]->p_lk_rght);
2210	/* Init_Tips_At_One_Site_Generic_Float(tree->data->c_seq[i]->state+curr_sitetree->mod->io->state_len, /
2211	/* tree->mod->ns, */
2212	/* tree->mod->io->state_len, */
2213	/* curr_sitedim1+0dim2, */
2214	/* tree->a_nodes[i]->b[0]->p_lk_rght); */
2215
2216	for(j=1;j<tree->mod->ras->n_catg;j++)
2217	{
2218	For(k,tree->mod->ns)
2219	{
2220	tree->a_nodes[i]->b[0]->p_lk_rght[curr_sitedim1+jdim2+k] =
2221	tree->a_nodes[i]->b[0]->p_lk_rght[curr_sitedim1+0dim2+k];
2222	}
2223	}
2224	}
2225	}
2226
2227	if(tree->mod->use_m4mod)
2228	{
2229	M4_Init_P_Lk_Tips_Double(tree);
2230	}
2231	}
2232
2233	//////////////////////////////////////////////////////////////
2234	//////////////////////////////////////////////////////////////
2235
2236
2237	void Init_P_Lk_Tips_Int(t_tree *tree)
2238	{
2239	int curr_site,i,dim1;
2240
2241	dim1 = tree->mod->ns;
2242
2243	For(i,tree->n_otu)
2244	{
2245	if(!tree->a_nodes[i]->c_seq \|\|
2246	strcmp(tree->a_nodes[i]->c_seq->name,tree->a_nodes[i]->name))
2247	{
2248	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
2249	Exit("");
2250	}
2251	}
2252
2253	For(curr_site,tree->data->crunch_len)
2254	{
2255	For(i,tree->n_otu)
2256	{
2257	/* printf("\n. site: %3d %c",curr_site,tree->a_nodes[i]->c_seq->state[curr_site]); */
2258	if(tree->io->datatype == NT)
2259	{
2260	Init_Tips_At_One_Site_Nucleotides_Int(tree->a_nodes[i]->c_seq->state[curr_site],
2261	curr_site*dim1,
2262	tree->a_nodes[i]->b[0]->p_lk_tip_r);
2263	/* Init_Tips_At_One_Site_Nucleotides_Int(tree->data->c_seq[i]->state[curr_site], */
2264	/* curr_sitedim1, /
2265	/* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2266	}
2267	else if(tree->io->datatype == AA)
2268	Init_Tips_At_One_Site_AA_Int(tree->a_nodes[i]->c_seq->state[curr_site],
2269	curr_site*dim1,
2270	tree->a_nodes[i]->b[0]->p_lk_tip_r);
2271	/* Init_Tips_At_One_Site_AA_Int(tree->data->c_seq[i]->state[curr_site], */
2272	/* curr_sitedim1, /
2273	/* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2274
2275	else if(tree->io->datatype == GENERIC)
2276	{
2277	Init_Tips_At_One_Site_Generic_Int(tree->a_nodes[i]->c_seq->state+curr_site*tree->mod->io->state_len,
2278	tree->mod->ns,
2279	tree->mod->io->state_len,
2280	curr_site*dim1,
2281	tree->a_nodes[i]->b[0]->p_lk_tip_r);
2282
2283	/* Init_Tips_At_One_Site_Generic_Int(tree->data->c_seq[i]->state+curr_sitetree->mod->io->state_len, /
2284	/* tree->mod->ns, */
2285	/* tree->mod->io->state_len, */
2286	/* curr_sitedim1, /
2287	/* tree->a_nodes[i]->b[0]->p_lk_tip_r); */
2288	}
2289	}
2290	}
2291	if(tree->mod->use_m4mod)
2292	{
2293	M4_Init_P_Lk_Tips_Int(tree);
2294	}
2295	}
2296
2297	//////////////////////////////////////////////////////////////
2298	//////////////////////////////////////////////////////////////
2299
2300	void Update_PMat_At_Given_Edge(t_edge b_fcus, t_tree tree)
2301	{
2302	int i;
2303	phydbl len;
2304	phydbl l_min, l_max;
2305	phydbl shape, scale, mean, var;
2306
2307	if(tree->is_mixt_tree)
2308	{
2309	MIXT_Update_PMat_At_Given_Edge(b_fcus,tree);
2310	return;
2311	}
2312
2313	l_min = tree->mod->l_min;
2314	l_max = tree->mod->l_max;
2315
2316	len = -1.0;
2317
2318	if(tree->mod->log_l == YES) b_fcus->l->v = EXP(b_fcus->l->v);
2319
2320	For(i,tree->mod->ras->n_catg)
2321	{
2322	if(tree->mod->ras->skip_rate_cat[i] == YES) continue;
2323
2324	if(b_fcus->has_zero_br_len == YES)
2325	{
2326	len = -1.0;
2327	mean = -1.0;
2328	var = -1.0;
2329	}
2330	else
2331	{
2332	len = MAX(0.0,b_fcus->l->v)*tree->mod->ras->gamma_rr->v[i];
2333	len *= tree->mod->br_len_multiplier->v;
2334	if(tree->mixt_tree) len *= tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number];
2335	if(len < l_min) len = l_min;
2336	else if(len > l_max) len = l_max;
2337
2338	mean = len;
2339	var = MAX(0.0,b_fcus->l_var->v) * POW(tree->mod->ras->gamma_rr->v[i]*tree->mod->br_len_multiplier->v,2);
2340	if(tree->mixt_tree) var *= POW(tree->mixt_tree->mod->ras->gamma_rr->v[tree->mod->ras->parent_class_number],2);
2341
2342	if(var > tree->mod->l_var_max) var = tree->mod->l_var_max;
2343	if(var < tree->mod->l_var_min) var = tree->mod->l_var_min;
2344	}
2345
2346
2347	if(tree->mod->gamma_mgf_bl == NO)
2348	PMat(len,tree->mod,tree->mod->nstree->mod->nsi,b_fcus->Pij_rr);
2349	else
2350	{
2351	shape = mean*mean/var;
2352	scale = var/mean;
2353
2354	PMat_MGF_Gamma(b_fcus->Pij_rr+tree->mod->nstree->mod->nsi,
2355	shape,scale,1.0,tree->mod);
2356
2357	}
2358	}
2359	if(tree->mod->log_l == YES) b_fcus->l->v = LOG(b_fcus->l->v);
2360	}
2361
2362	//////////////////////////////////////////////////////////////
2363	//////////////////////////////////////////////////////////////
2364
2365
2366	/* void Update_P_Lk_On_A_Path(t_node a, t_node d, t_edge b, t_node target_one_side, t_node target_other_side, t_tree tree) */
2367	/* { */
2368
2369
2370	/* /\* */
2371	/* \ / */
2372	/* target\___________b_/ */
2373	/* / \ \ \ \ */
2374	/* / \ \ \ \ */
2375
2376	/* target is the root of the subtree at which we want */
2377	/* the likelihood to be updated */
2378	/* \/ /
2379
2380
2381
2382	/* /\* PhyML_Printf("Update_p_lk on (%d %d) at %d (target=%d %d)\n", \/ /
2383	/* /\* b->left->num, \/ /
2384	/* /\* b->rght->num, \/ /
2385	/* /\* a->num, \/ /
2386	/* /\* target_one_side->num, \/ /
2387	/* /\* target_other_side->num); \/ /
2388
2389	/* Update_P_Lk(tree,b,a); */
2390	/* if((a == target_one_side) && (d == target_other_side)) */
2391	/* return; */
2392	/* else */
2393	/* { */
2394	/* Update_P_Lk_On_A_Path(d, */
2395	/* d->v[tree->t_dir[d->num][target_other_side->num]], */
2396	/* d->b[tree->t_dir[d->num][target_other_side->num]], */
2397	/* target_one_side, */
2398	/* target_other_side, */
2399	/* tree); */
2400	/* } */
2401	/* } */
2402
2403	void Update_P_Lk_Along_A_Path(t_node *path, int path_length, t_tree tree)
2404	{
2405	int i,j;
2406
2407	For(i,path_length-1)
2408	{
2409	For(j,3)
2410	if(path[i]->v[j] == path[i+1])
2411	{
2412	if(path[i] == path[i]->b[j]->left)
2413	{
2414	Update_P_Lk(tree,path[i]->b[j],path[i]->b[j]->left);
2415	}
2416	else if(path[i] == path[i]->b[j]->rght)
2417	{
2418	Update_P_Lk(tree,path[i]->b[j],path[i]->b[j]->rght);
2419	}
2420	else
2421	{
2422	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
2423	Exit("");
2424	}
2425	break;
2426	}
2427	#ifdef DEBUG
2428	if(j == 3)
2429	{
2430	PhyML_Printf("\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
2431	Exit("");
2432	}
2433	#endif
2434	}
2435	}
2436
2437	//////////////////////////////////////////////////////////////
2438	//////////////////////////////////////////////////////////////
2439
2440
2441	phydbl Lk_Dist(phydbl F, phydbl dist, t_mod mod)
2442	{
2443	int i,j,k;
2444	phydbl lnL,len;
2445	int dim1,dim2;
2446
2447	if(mod->log_l == YES) dist = EXP(dist);
2448
2449	For(k,mod->ras->n_catg)
2450	{
2451	len = dist*mod->ras->gamma_rr->v[k];
2452	if(len < mod->l_min) len = mod->l_min;
2453	else if(len > mod->l_max) len = mod->l_max;
2454	PMat(len,mod,mod->nsmod->nsk,mod->Pij_rr->v);
2455	}
2456
2457	dim1 = mod->ns*mod->ns;
2458	dim2 = mod->ns;
2459	lnL = .0;
2460
2461	/* For(i,mod->ns) */
2462	/* { */
2463	/* For(j,mod->ns) */
2464	/* { */
2465	/* For(k,mod->ras->n_catg) */
2466	/* { */
2467	/* lnL += */
2468	/* F[dim1k+dim2i+j] * */
2469	/* LOG(mod->pi[i] * mod->Pij_rr[dim1k+dim2i+j]); */
2470	/* } */
2471	/* } */
2472	/* } */
2473
2474	For(i,mod->ns-1)
2475	{
2476	for(j=i+1;j<mod->ns;j++)
2477	{
2478	For(k,mod->ras->n_catg)
2479	{
2480	lnL +=
2481	(F[dim1k+dim2i+j] + F[dim1k+dim2j+i])*
2482	LOG(mod->e_frq->pi->v[i] * mod->Pij_rr->v[dim1k+dim2i+j]);
2483	/* printf("\n. f: %f Pij:%f F:%f", */
2484	/* mod->e_frq->pi->v[i], */
2485	/* mod->Pij_rr->v[dim1k+dim2i+j], */
2486	/* F[dim1k+dim2j+i]); */
2487	}
2488	}
2489	}
2490
2491	For(i,mod->ns) For(k,mod->ras->n_catg) lnL += F[dim1k+dim2i+i]* LOG(mod->e_frq->pi->v[i] * mod->Pij_rr->v[dim1k+dim2i+i]);
2492
2493	return lnL;
2494	}
2495
2496	//////////////////////////////////////////////////////////////
2497	//////////////////////////////////////////////////////////////
2498
2499
2500	phydbl Update_Lk_At_Given_Edge(t_edge b_fcus, t_tree tree)
2501	{
2502	Update_P_Lk(tree,b_fcus,b_fcus->left);
2503	Update_P_Lk(tree,b_fcus,b_fcus->rght);
2504	tree->c_lnL = Lk(b_fcus,tree);
2505	return tree->c_lnL;
2506	}
2507
2508	//////////////////////////////////////////////////////////////
2509	//////////////////////////////////////////////////////////////
2510
2511
2512	void Print_Lk_Given_Edge_Recurr(t_node a, t_node d, t_edge b, t_tree tree)
2513	{
2514	PhyML_Printf("\n___ Edge %3d (left=%3d rght=%3d) lnL=%f",
2515	b->num,
2516	b->left->num,
2517	b->rght->num,
2518	Lk(b,tree));
2519
2520	if(d->tax) return;
2521	else
2522	{
2523	int i;
2524	For(i,3)
2525	if(d->v[i] != a)
2526	Print_Lk_Given_Edge_Recurr(d,d->v[i],d->b[i],tree);
2527	}
2528	}
2529
2530	//////////////////////////////////////////////////////////////
2531	//////////////////////////////////////////////////////////////
2532
2533	//////////////////////////////////////////////////////////////
2534	//////////////////////////////////////////////////////////////
2535
2536	//////////////////////////////////////////////////////////////
2537	//////////////////////////////////////////////////////////////
2538
2539	//////////////////////////////////////////////////////////////
2540	//////////////////////////////////////////////////////////////
2541
2542	//////////////////////////////////////////////////////////////
2543	//////////////////////////////////////////////////////////////
2544
2545	//////////////////////////////////////////////////////////////
2546	//////////////////////////////////////////////////////////////
2547
2548	//////////////////////////////////////////////////////////////
2549	//////////////////////////////////////////////////////////////
2550
2551	//////////////////////////////////////////////////////////////
2552	//////////////////////////////////////////////////////////////
2553
2554	//////////////////////////////////////////////////////////////
2555	//////////////////////////////////////////////////////////////
2556
2557	//////////////////////////////////////////////////////////////
2558	//////////////////////////////////////////////////////////////
2559
2560
2561	//////////////////////////////////////////////////////////////
2562	//////////////////////////////////////////////////////////////
2563
2564
2565	void Alias_Subpatt(t_tree *tree)
2566	{
2567
2568	if(tree->n_root)
2569	{
2570	Alias_Subpatt_Post(tree->n_root,tree->n_root->v[2],tree);
2571	Alias_Subpatt_Post(tree->n_root,tree->n_root->v[1],tree);
2572	}
2573	else
2574	{
2575	Alias_Subpatt_Post(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
2576	/* if(tree->both_sides) */
2577	Alias_Subpatt_Pre(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree);
2578	}
2579	}
2580
2581	//////////////////////////////////////////////////////////////
2582	//////////////////////////////////////////////////////////////
2583
2584
2585	void Alias_One_Subpatt(t_node a, t_node d, t_tree *tree)
2586	{
2587	int i,j;
2588	int patt_id_v1, patt_id_v2, *patt_id_d;
2589	int p_lk_loc_d, p_lk_loc_v1, *p_lk_loc_v2;
2590	t_node v1, v2;
2591	t_edge b0, b1, *b2;
2592	int curr_patt_id_v1, curr_patt_id_v2;
2593	int curr_p_lk_loc_v1, curr_p_lk_loc_v2;
2594	int num_subpatt;
2595
2596	b0 = b1 = b2 = NULL;
2597
2598	if(d->tax)
2599	{
2600	patt_id_d = (d == d->b[0]->left)?(d->b[0]->patt_id_left):(d->b[0]->patt_id_rght);
2601	p_lk_loc_d = (d == d->b[0]->left)?(d->b[0]->p_lk_loc_left):(d->b[0]->p_lk_loc_rght);
2602
2603	For(i,tree->n_pattern)
2604	{
2605	For(j,tree->n_pattern)
2606	{
2607	if(patt_id_d[i] == patt_id_d[j])
2608	{
2609	p_lk_loc_d[i] = j;
2610	break;
2611	}
2612	if(j > i)
2613	{
2614	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2615	Warn_And_Exit("");
2616	}
2617	}
2618	}
2619	return;
2620	}
2621	else
2622	{
2623	v1 = v2 = NULL;
2624	For(i,3)
2625	{
2626	if(d->v[i] != a && d->b[i] != tree->e_root)
2627	{
2628	if(!v1) { v1=d->v[i]; b1=d->b[i];}
2629	else { v2=d->v[i]; b2=d->b[i];}
2630	}
2631	else
2632	{
2633	b0 = d->b[i];
2634	}
2635	}
2636
2637
2638	patt_id_v1 = (v1 == b1->left)?(b1->patt_id_left):(b1->patt_id_rght);
2639	patt_id_v2 = (v2 == b2->left)?(b2->patt_id_left):(b2->patt_id_rght);
2640	patt_id_d = (d == b0->left)?(b0->patt_id_left):(b0->patt_id_rght);
2641	p_lk_loc_d = (d == b0->left)?(b0->p_lk_loc_left):(b0->p_lk_loc_rght);
2642	p_lk_loc_v1 = (v1 == b1->left)?(b1->p_lk_loc_left):(b1->p_lk_loc_rght);
2643	p_lk_loc_v2 = (v2 == b2->left)?(b2->p_lk_loc_left):(b2->p_lk_loc_rght);
2644
2645	num_subpatt = 0;
2646	For(i,tree->n_pattern)
2647	{
2648	curr_patt_id_v1 = patt_id_v1[i];
2649	curr_patt_id_v2 = patt_id_v2[i];
2650	curr_p_lk_loc_v1 = p_lk_loc_v1[i];
2651	curr_p_lk_loc_v2 = p_lk_loc_v2[i];
2652
2653	p_lk_loc_d[i] = i;
2654
2655	if((curr_p_lk_loc_v1 == i) \|\| (curr_p_lk_loc_v2 == i))
2656	{
2657	p_lk_loc_d[i] = i;
2658	patt_id_d[i] = num_subpatt;
2659	num_subpatt++;
2660	}
2661	else
2662	if(curr_p_lk_loc_v1 == curr_p_lk_loc_v2)
2663	{
2664	p_lk_loc_d[i] = curr_p_lk_loc_v1;
2665	patt_id_d[i] = patt_id_d[curr_p_lk_loc_v1];
2666	}
2667	else
2668	{
2669	for(j=MAX(curr_p_lk_loc_v1,curr_p_lk_loc_v2);j<tree->n_pattern;j++)
2670	{
2671	if((patt_id_v1[j] == curr_patt_id_v1) &&
2672	(patt_id_v2[j] == curr_patt_id_v2))
2673	{
2674	p_lk_loc_d[i] = j;
2675
2676	if(j == i)
2677	{
2678	patt_id_d[i] = num_subpatt;
2679	num_subpatt++;
2680	}
2681	else patt_id_d[i] = patt_id_d[j];
2682	break;
2683	}
2684	if(j > i)
2685	{
2686	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
2687	Warn_And_Exit("");
2688	}
2689	}
2690	}
2691	}
2692	}
2693	}
2694
2695	//////////////////////////////////////////////////////////////
2696	//////////////////////////////////////////////////////////////
2697
2698	void Alias_Subpatt_Post(t_node a, t_node d, t_tree *tree)
2699	{
2700	if(d->tax) return;
2701	else
2702	{
2703	int i;
2704
2705	For(i,3)
2706	{
2707	if(d->v[i] != a && d->b[i] != tree->e_root)
2708	{
2709	Alias_Subpatt_Post(d,d->v[i],tree);
2710	}
2711	}
2712	Alias_One_Subpatt(a, d, tree);
2713	}
2714	}
2715
2716	//////////////////////////////////////////////////////////////
2717	//////////////////////////////////////////////////////////////
2718
2719
2720	void Alias_Subpatt_Pre(t_node a, t_node d, t_tree *tree)
2721	{
2722	if(d->tax) return;
2723	else
2724	{
2725	int i;
2726
2727	For(i,3)
2728	{
2729	if(d->v[i] != a && d->b[i] != tree->e_root)
2730	{
2731	Alias_One_Subpatt(d->v[i],d,tree);
2732	Alias_Subpatt_Pre(d,d->v[i],tree);
2733	}
2734	}
2735	}
2736	}
2737
2738	//////////////////////////////////////////////////////////////
2739	//////////////////////////////////////////////////////////////
2740
2741
2742	void Copy_P_Lk(phydbl p_lk, int site_from, int site_to, t_tree tree)
2743	{
2744	int i,j;
2745	int dim1,dim2;
2746
2747
2748	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
2749	dim2 = tree->mod->ns;
2750
2751	/* PhyML_Printf("\n# %d %d",site_to,site_from); */
2752
2753	For(i,tree->mod->ns) For(j,tree->mod->ras->n_catg)
2754	{
2755	p_lk[site_todim1+jdim2+i] = p_lk[site_fromdim1+jdim2+i];
2756	}
2757	}
2758
2759	//////////////////////////////////////////////////////////////
2760	//////////////////////////////////////////////////////////////
2761
2762
2763	void Copy_Scale(int scale, int site_from, int site_to, t_tree tree)
2764	{
2765	int i;
2766
2767	For(i,tree->mod->ras->n_catg)
2768	{
2769	scale[itree->n_pattern+site_to] = scale[itree->n_pattern+site_from];
2770	/* PhyML_Printf("\n. %d",scale[itree->n_pattern+site_to]); /
2771	}
2772	}
2773
2774	//////////////////////////////////////////////////////////////
2775	//////////////////////////////////////////////////////////////
2776
2777
2778	void Init_P_Lk_Loc(t_tree *tree)
2779	{
2780	int i,j;
2781	t_node *d;
2782	int *patt_id_d;
2783
2784	For(i,2*tree->n_otu-1)
2785	{
2786	For(j,tree->n_pattern)
2787	{
2788	tree->a_edges[i]->p_lk_loc_left[j] = j;
2789	tree->a_edges[i]->p_lk_loc_rght[j] = j;
2790	}
2791	}
2792
2793	For(i,tree->n_otu)
2794	{
2795	d = tree->a_nodes[i];
2796	patt_id_d = (d == d->b[0]->left)?(d->b[0]->patt_id_left):(d->b[0]->patt_id_rght);
2797	For(j,tree->n_pattern)
2798	{
2799	patt_id_d[j] = (int)tree->a_nodes[d->num]->c_seq->state[j];
2800	/* patt_id_d[j] = (int)tree->data->c_seq[d->num]->state[j]; */
2801	}
2802	}
2803	}
2804
2805	//////////////////////////////////////////////////////////////
2806	//////////////////////////////////////////////////////////////
2807
2808
2809	phydbl Lk_Normal_Approx(t_tree *tree)
2810	{
2811	phydbl lnL;
2812	int i;
2813	int dim;
2814	phydbl first_order;
2815
2816	dim = 2*tree->n_otu-3;
2817
2818	lnL = Dnorm_Multi_Given_InvCov_Det(tree->rates->u_cur_l,
2819	tree->rates->mean_l,
2820	tree->rates->invcov,
2821	tree->rates->covdet,
2822	2*tree->n_otu-3,YES);
2823
2824	first_order = 0.0;
2825	For(i,dim) first_order += (tree->rates->u_cur_l[i] - tree->rates->mean_l[i]) * tree->rates->grad_l[i];
2826
2827	lnL += first_order;
2828
2829	/* printf("\n"); */
2830	/* For(i,dim) printf("%f\t",tree->rates->u_cur_l[i]); */
2831	/* printf("\n. Lk=%f %f",lnL,tree->mod->l_min); */
2832
2833
2834	/* err = NO; */
2835	/* dim = 2tree->n_otu-3; /
2836	/* For(i,dim) */
2837	/* { */
2838	/* if((tree->rates->mean_l[i] / tree->mod->l_min < 1.1) && */
2839	/* (tree->rates->mean_l[i] / tree->mod->l_min > 0.9)) */
2840	/* { */
2841	/* lnL -= Log_Dnorm(tree->a_edges[i]->l->v,tree->rates->mean_l[i],SQRT(tree->rates->cov_l[idim+i]),&err); /
2842	/* if(err) */
2843	/* { */
2844	/* PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__); */
2845	/* Warn_And_Exit(""); */
2846	/* } */
2847	/* lambda = 1./SQRT(tree->rates->cov_l[idim+i]); /
2848	/* l = (tree->mod->log_l == YES)?(EXP(tree->a_edges[i]->l->v)):(tree->a_edges[i]->l->v); */
2849	/* lnL += LOG(Dexp(l,lambda)); */
2850	/* /\* printf("\n. lambda = %f",lambda); \/ /
2851	/* } */
2852	/* } */
2853
2854	return(lnL);
2855
2856	}
2857
2858	//////////////////////////////////////////////////////////////
2859	//////////////////////////////////////////////////////////////
2860
2861
2862	phydbl Wrap_Part_Lk_At_Given_Edge(t_edge b, t_tree tree, supert_tree *stree)
2863	{
2864	return PART_Lk_At_Given_Edge(b,stree);;
2865	}
2866
2867	//////////////////////////////////////////////////////////////
2868	//////////////////////////////////////////////////////////////
2869
2870
2871	phydbl Wrap_Part_Lk(t_edge b, t_tree tree, supert_tree *stree)
2872	{
2873	return PART_Lk(stree);
2874	}
2875
2876	//////////////////////////////////////////////////////////////
2877	//////////////////////////////////////////////////////////////
2878
2879
2880	phydbl Wrap_Lk(t_edge b, t_tree tree, supert_tree *stree)
2881	{
2882	Lk(NULL,tree);
2883	return tree->c_lnL;
2884	}
2885
2886	//////////////////////////////////////////////////////////////
2887	//////////////////////////////////////////////////////////////
2888
2889
2890	phydbl Wrap_Geo_Lk(t_edge b, t_tree tree, supert_tree *stree)
2891	{
2892	TIPO_Lk(tree);
2893	return tree->geo_lnL;
2894	}
2895
2896	//////////////////////////////////////////////////////////////
2897	//////////////////////////////////////////////////////////////
2898
2899
2900	phydbl Wrap_Lk_At_Given_Edge(t_edge b, t_tree tree, supert_tree *stree)
2901	{
2902	Lk(b,tree);
2903	return tree->c_lnL;
2904	}
2905
2906	//////////////////////////////////////////////////////////////
2907	//////////////////////////////////////////////////////////////
2908
2909
2910	phydbl Wrap_Lk_Rates(t_edge b, t_tree tree, supert_tree *stree)
2911	{
2912	RATES_Lk_Rates(tree);
2913	return tree->rates->c_lnL_rates;
2914	}
2915
2916	//////////////////////////////////////////////////////////////
2917	//////////////////////////////////////////////////////////////
2918
2919	phydbl Wrap_Lk_Times(t_edge b, t_tree tree, supert_tree *stree)
2920	{
2921	TIMES_Lk_Times(tree);
2922	return tree->rates->c_lnL_times;
2923	}
2924
2925	//////////////////////////////////////////////////////////////
2926	//////////////////////////////////////////////////////////////
2927
2928
2929
2930	phydbl Wrap_Lk_Linreg(t_edge b, t_tree tree, supert_tree *stree)
2931	{
2932	/* RATES_Lk_Linreg(tree); */
2933	return -1.;
2934	/* return tree->rates->c_lnL_linreg; */
2935	}
2936
2937	//////////////////////////////////////////////////////////////
2938	//////////////////////////////////////////////////////////////
2939
2940
2941	phydbl Wrap_Diff_Lk_Norm_At_Given_Edge(t_edge b, t_tree tree, supert_tree *stree)
2942	{
2943	phydbl diff;
2944	diff = Diff_Lk_Norm_At_Given_Edge(b,tree);
2945	return(-diff);
2946
2947	}
2948
2949	//////////////////////////////////////////////////////////////
2950	//////////////////////////////////////////////////////////////
2951
2952	void Sample_Ancestral_Seq(int mutmap, int fromprior, t_tree *tree)
2953	{
2954	int rate_cat;
2955	int i,j,k,l;
2956	phydbl *probs;
2957	phydbl sum;
2958	phydbl u;
2959	int n_mut;
2960	FILE *fp;
2961	phydbl *muttime;
2962	int *muttype;
2963	char *s;
2964	int *ordering;
2965
2966	probs = (phydbl *)mCalloc(tree->mod->ras->n_catg,sizeof(phydbl));
2967
2968	muttime = (phydbl )mCalloc((2tree->n_otu-3)*5, // At most 5 mutations per branch on average
2969	sizeof(phydbl));
2970
2971	muttype = (int )mCalloc((2tree->n_otu-3)*5, // At most 5 mutations per branch on average
2972	sizeof(int));
2973
2974	ordering = (int )mCalloc((2tree->n_otu-3)*5, // At most 5 mutations per branch on average
2975	sizeof(int));
2976
2977	s = (char *)mCalloc(T_MAX_NAME,sizeof(char));
2978
2979	if(fromprior == YES)
2980	{
2981	/* Update P(D_x\|X=i) for each state i and node X */
2982	Set_Both_Sides(YES,tree);
2983	Lk(NULL,tree);
2984	}
2985
2986	For(i,tree->n_pattern)
2987	{
2988	/* Sample the rate class from its posterior density */
2989	For(j,tree->mod->ras->n_catg)
2990	{
2991	if(fromprior == NO)
2992	probs[j] =
2993	EXP(tree->log_site_lk_cat[j][i])*
2994	tree->mod->ras->gamma_r_proba->v[j];
2995	else
2996	probs[j] = tree->mod->ras->gamma_r_proba->v[j];
2997	}
2998
2999
3000	/* Scale probas. */
3001	sum = .0;
3002	For(j,tree->mod->ras->n_catg) sum += probs[j];
3003	For(j,tree->mod->ras->n_catg) probs[j]/=sum;
3004
3005	/* CDF */
3006	for(j=1;j<tree->mod->ras->n_catg;j++) probs[j] += probs[j-1];
3007
3008	/* Sample rate */
3009	u = Uni();
3010	rate_cat = -1;
3011	For(j,tree->mod->ras->n_catg)
3012	if(probs[j] > u)
3013	{
3014	rate_cat = j;
3015	break;
3016	}
3017
3018	n_mut = 0;
3019	Sample_Ancestral_Seq_Pre(tree->a_nodes[0],tree->a_nodes[0]->v[0],tree->a_nodes[0]->b[0],
3020	i,rate_cat,
3021	muttype,muttime,&n_mut,
3022	mutmap,fromprior,tree);
3023
3024
3025	For(j,n_mut) ordering[j] = 0;
3026
3027	For(j,n_mut-1)
3028	{
3029	for(k=j+1;k<n_mut;k++)
3030	{
3031	if(muttime[k] > muttime[j]) ordering[k]++;
3032	else ordering[j]++;
3033	}
3034	}
3035
3036	strcpy(s,"rosetta.");
3037	sprintf(s+strlen(s),"%d",i);
3038	fp = fopen(s,"a");
3039	PhyML_Fprintf(fp,"\n-1 -1 -1.0 -1");
3040
3041	For(j,n_mut)
3042	{
3043	For(k,n_mut)
3044	{
3045	if(ordering[k] == j)
3046	{
3047	For(l,tree->data->init_len) if(tree->data->sitepatt[l] == i) break;
3048	PhyML_Fprintf(fp,"\n%4d %4d %12f %4d",j,muttype[k],muttime[k],l);
3049	/* PhyML_Fprintf(fp,"\n%d",muttype[ordering[j]]); */
3050	break;
3051	}
3052	}
3053	}
3054
3055
3056	For(j,n_mut)
3057	{
3058	muttype[j] = -2;
3059	muttime[j] = +1.;
3060	}
3061
3062	fclose(fp);
3063	}
3064
3065	Free(s);
3066	Free(muttype);
3067	Free(muttime);
3068	Free(ordering);
3069	Free(probs);
3070	}
3071
3072	//////////////////////////////////////////////////////////////
3073	//////////////////////////////////////////////////////////////
3074
3075	void Sample_Ancestral_Seq_Pre(t_node a, t_node d, t_edge *b,
3076	int site, int rate_cat,
3077	int muttype, phydbl muttime, int *n_mut,
3078	int mutmap, int fromprior, t_tree *tree)
3079	{
3080
3081	int i,j;
3082	int sa,sd;
3083	phydbl *Pij;
3084	phydbl *p_lk;
3085	int dim1, dim2, dim3;
3086	phydbl sum;
3087	phydbl u;
3088	char *c;
3089	phydbl *probs;
3090
3091	probs = (phydbl *)mCalloc(tree->mod->ns,sizeof(phydbl));
3092
3093	if(a->tax)
3094	c = a->c_seq->state+site*tree->mod->io->state_len;
3095	/* c = tree->data->c_seq[a->num]->state+sitetree->mod->io->state_len; /
3096	else
3097	c = a->c_seq_anc->state+site*tree->mod->io->state_len;
3098	/* c = tree->anc_data->c_seq[a->num-tree->n_otu]->state+sitetree->mod->io->state_len; /
3099
3100	sa = Assign_State(c,
3101	tree->mod->io->datatype,
3102	tree->mod->io->state_len);
3103
3104	if(sa == -1) // c is an indel
3105	{
3106	For(j,tree->mod->ns) probs[j] = tree->mod->e_frq->pi->v[j];
3107
3108	for(j=1;j<tree->mod->ns;j++) probs[j] += probs[j-1];
3109
3110	u = Uni();
3111	For(j,tree->mod->ns)
3112	if(probs[j] > u)
3113	{
3114	sa = j;
3115	break;
3116	}
3117	}
3118
3119	if(d->tax == NO) // Need to sample state at node d
3120	{
3121
3122	dim1 = tree->mod->ras->n_catg * tree->mod->ns;
3123	dim2 = tree->mod->ns;
3124	dim3 = tree->mod->ns * tree->mod->ns;
3125	sum = 0.0;
3126
3127
3128	Pij = b->Pij_rr;
3129
3130	if(d == b->left)
3131	p_lk = b->p_lk_left;
3132	else
3133	p_lk = b->p_lk_rght;
3134
3135	For(j,tree->mod->ns) probs[j] = 0.0;
3136
3137	/* Formula (10) in Nielsen's Mutation Maping paper, e.g. */
3138	For(j,tree->mod->ns)
3139	{
3140	if(fromprior == NO)
3141	probs[j] =
3142	p_lk[sitedim1+rate_catdim2+j] *
3143	Pij[rate_catdim3+sadim2+j];
3144	else
3145	probs[j] = Pij[rate_catdim3+sadim2+j];
3146	}
3147
3148
3149	/* if(site == 92) */
3150	/* printf("\n. l=%f pr[%f %f %f %f] Pij[%f %f %f %f] plk[%f %f %f %f]", */
3151	/* b->l->v * tree->mod->ras->gamma_rr->v[rate_cat], */
3152	/* probs[0],probs[1],probs[2],probs[3], */
3153	/* Pij[rate_catdim3+sadim2+0], */
3154	/* Pij[rate_catdim3+sadim2+1], */
3155	/* Pij[rate_catdim3+sadim2+2], */
3156	/* Pij[rate_catdim3+sadim2+3], */
3157	/* p_lk[sitedim1+rate_catdim2+0], */
3158	/* p_lk[sitedim1+rate_catdim2+1], */
3159	/* p_lk[sitedim1+rate_catdim2+2], */
3160	/* p_lk[sitedim1+rate_catdim2+3]); */
3161
3162
3163	/* Scale the probabilities */
3164	sum = 0.0;
3165	For(j,tree->mod->ns) sum += probs[j];
3166	For(j,tree->mod->ns) probs[j] /= sum;
3167
3168	/* CDF */
3169	for(j=1;j<tree->mod->ns;j++) probs[j] += probs[j-1];
3170
3171	/* Sample state according to their posterior probas. */
3172	sd = -1;
3173	u = Uni();
3174	For(j,tree->mod->ns)
3175	if(probs[j] > u)
3176	{
3177	sd = j;
3178	break;
3179	}
3180
3181	/* Assign state */
3182	/* tree->anc_data->c_seq[d->num-tree->n_otu]->state[site] = Reciproc_Assign_State(sd,tree->io->datatype); */
3183	d->c_seq_anc->state[site] = Reciproc_Assign_State(sd,tree->io->datatype);
3184	}
3185	else
3186	{
3187	c = d->c_seq->state+site*tree->mod->io->state_len;
3188	/* c = tree->data->c_seq[d->num]->state+sitetree->mod->io->state_len; /
3189
3190	sd = Assign_State(c,
3191	tree->mod->io->datatype,
3192	tree->mod->io->state_len);
3193
3194	if(sd == -1) // c is an indel
3195	{
3196	For(j,tree->mod->ns) probs[j] = tree->mod->e_frq->pi->v[j];
3197
3198	for(j=1;j<tree->mod->ns;j++) probs[j] += probs[j-1];
3199
3200	u = Uni();
3201	For(j,tree->mod->ns)
3202	if(probs[j] > u)
3203	{
3204	sd = j;
3205	break;
3206	}
3207	}
3208	}
3209
3210	/* if(site == 92) */
3211	/* { */
3212	/* printf("\n. sa=%d (%s,%d) sd=%d (%s,%d) b->l->v=%f", */
3213	/* sa,a->tax?a->name:"",a->num,sd,d->tax?d->name:"",d->num,b->l->v); */
3214	/* } */
3215
3216	if(mutmap == YES) Map_Mutations(a,d,sa,sd,b,site,rate_cat,muttype,muttime,n_mut,tree);
3217
3218	Free(probs);
3219
3220	if(d->tax) return;
3221	else
3222	{
3223	For(i,3)
3224	{
3225	if(d->v[i] != a)
3226	{
3227	Sample_Ancestral_Seq_Pre(d,d->v[i],d->b[i],site,rate_cat,muttype,muttime,n_mut,mutmap,fromprior,tree);
3228	}
3229	}
3230	}
3231	}
3232
3233	//////////////////////////////////////////////////////////////
3234	//////////////////////////////////////////////////////////////
3235
3236	void Map_Mutations(t_node a, t_node d, int sa, int sd, t_edge b, int site, int rate_cat, int muttype, phydbl muttime, int n_mut, t_tree *tree)
3237	{
3238	int i,j;
3239	phydbl probs,all_probs;
3240	int slast; // Last state visited
3241	phydbl tlast, td;
3242	phydbl *Q;
3243	phydbl u,sum;
3244	int *mut; // Array of mutations
3245	int thismut;
3246	int n_mut_branch;
3247	phydbl br,cr,ta,gr;
3248	int n_iter;
3249	int first_mut;
3250
3251	all_probs = (phydbl )mCalloc(tree->mod->nstree->mod->ns,sizeof(phydbl));
3252	mut = (int )mCalloc(tree->mod->nstree->mod->ns,sizeof(int));
3253
3254	// Edge rate
3255	br =
3256	(tree->rates->model_log_rates == YES)?
3257	EXP(tree->rates->br_r[d->num]):
3258	tree->rates->br_r[d->num];
3259
3260	// Clock (i.e., overall) rate
3261	cr = tree->rates->clock_r;
3262
3263	// Age of node a
3264	ta = tree->rates->nd_t[a->num];
3265
3266	// Site relative rate
3267	gr = tree->mod->ras->gamma_rr->v[rate_cat];
3268
3269
3270	// Rate matrix
3271	Q = tree->mod->r_mat->qmat->v;
3272
3273	// Length of the 'time' interval considered: product of the branch length by the
3274	// current relative rate at that site (set when sampling ancestral sequences)
3275	td = b->l->v * gr;
3276
3277	// Matrix of change probabilities
3278	For(i,tree->mod->ns)
3279	{
3280	// We only care about the non-diagonal elements here
3281	For(j,tree->mod->ns) all_probs[itree->mod->ns+j] = -Q[itree->mod->ns+j] / Q[i*tree->mod->ns+i];
3282
3283	// Set the diagonal to 0
3284	all_probs[i*tree->mod->ns+i] = 0.0;
3285
3286	// \sum_{j != i} -q_{ij}/q_{ii}
3287	sum = 0;
3288	For(j,tree->mod->ns) sum += all_probs[i*tree->mod->ns+j];
3289
3290	// Diagonal: 1 - \sum_{j != i} -q_{ij}/q_{ii}
3291	all_probs[i*tree->mod->ns+i] = 1.-sum;
3292
3293	// Get the cumulative probas
3294	for(j=1;j<tree->mod->ns;j++) all_probs[itree->mod->ns+j] += all_probs[itree->mod->ns+j-1];
3295	}
3296
3297	For(i,tree->mod->ns*tree->mod->ns) mut[i] = 0;
3298	tlast = .0;
3299	slast = sa;
3300	probs = NULL;
3301	n_mut_branch = 0;
3302	n_iter = 0;
3303	first_mut = YES;
3304
3305	do
3306	{
3307	if((sa != sd) && (first_mut == YES)) // ancestral and descendant states are distinct
3308	{
3309	// Sample a time for the first mutation conditional on at least one mutation
3310	// occurring (see formula A2 in Nielsen's Genetics paper (2001))
3311	u = Uni();
3312	tlast = -LOG(1. - u(1.-EXP(Q[satree->mod->ns+sa]td)))/-Q[satree->mod->ns+sa];
3313	}
3314	else
3315	{
3316	// Sample a time for the next mutation
3317	tlast = tlast + Rexp(-Q[slast*tree->mod->ns+slast]);
3318	}
3319
3320	// Select the appropriate vector of change probabilities
3321	probs = all_probs+slast*tree->mod->ns;
3322
3323	/* printf("\n. slast=%2d sd=%2d tlast=%12G td=%12G p=%12f rcat=%12f site=%4d", */
3324	/* slast,sd,tlast,td,-Q[slasttree->mod->ns+slast], /
3325	/* tree->mod->ras->gamma_rr->v[rate_cat],site); */
3326
3327	// The time for the next mutation does not exceed the length
3328	// of the time interval -> sample a new mutation event
3329	if(tlast < td)
3330	{
3331	first_mut = NO;
3332
3333	n_mut_branch++;
3334
3335	u = Uni();
3336	For(i,tree->mod->ns)
3337	if(probs[i] > u)
3338	{
3339	// Record mutation type
3340	mut[slast*tree->mod->ns+i]++;
3341
3342	// Record mutation type in the site mutation array
3343	thismut = MIN(i,slast) * tree->mod->ns + MAX(i,slast) - (MIN(i,slast)+1+(int)POW(MIN(i,slast)+1,2))/2;
3344	muttype[(*n_mut)+n_mut_branch-1] = thismut;
3345
3346	if((thismut > 5) \|\| (thismut < 0))
3347	{
3348	PhyML_Printf("\n. thismut = %d",thismut);
3349	PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
3350	Warn_And_Exit("");
3351	}
3352
3353	// Record time of mutation
3354	muttime[(n_mut)+n_mut_branch-1] = ta + brcr*gr;
3355
3356	// Update the last state
3357	slast = i;
3358	break;
3359	}
3360	}
3361	else
3362	{
3363	if(slast == sd) break;
3364	else
3365	{
3366	// Restart from the beginning
3367	For(i,tree->mod->ns*tree->mod->ns) mut[i] = 0;
3368	For(i,n_mut_branch) muttype[(*n_mut)+n_mut_branch-1] = -2;
3369	For(i,n_mut_branch) muttime[(*n_mut)+n_mut_branch-1] = +1.;
3370	tlast = 0.0;
3371	slast = sa;
3372	n_mut_branch = 0;
3373	first_mut = YES;
3374	n_iter++;
3375	}
3376	}
3377	}
3378	while(1);
3379
3380	(*n_mut) += n_mut_branch;
3381
3382
3383	For(i,tree->mod->ns)
3384	{
3385	for(j=i+1;j<tree->mod->ns;j++)
3386	{
3387	if(mut[itree->mod->ns+j] + mut[jtree->mod->ns+i] > 0)
3388	{
3389	thismut = MIN(i,j) * tree->mod->ns + MAX(i,j) - (MIN(i,j)+1+(int)POW(MIN(i,j)+1,2))/2;
3390	tree->mutmap[thismut(tree->n_pattern)(2tree->n_otu-3) + b->num(tree->n_pattern) + site]++;
3391	/* if(site == 92) */
3392	/* { */
3393	/* printf("\nx sa=%d sd=%d mut=%d",sa,sd,thismut); */
3394	/* } */
3395	}
3396	}
3397	}
3398
3399	Free(all_probs);
3400	Free(mut);
3401	}
3402
3403	//////////////////////////////////////////////////////////////
3404	//////////////////////////////////////////////////////////////
3405
3406	int Check_Lk_At_Given_Edge(int verbose, t_tree *tree)
3407	{
3408	int res;
3409	int i;
3410	phydbl *lk;
3411
3412	lk = (phydbl )mCalloc(2tree->n_otu-3,sizeof(phydbl));
3413
3414	res = 0;
3415	For(i,2*tree->n_otu-3)
3416	{
3417	lk[i] = Lk(tree->a_edges[i],tree);
3418	if(verbose == YES) PhyML_Printf("\n. Edge %3d %13G %f %13G",
3419	tree->a_edges[i]->num,tree->a_edges[i]->l->v,lk[i],
3420	tree->a_edges[i]->l_var->v);
3421	}
3422
3423	if(tree->n_root)
3424	{
3425	Lk(tree->n_root->b[1],tree);
3426	if(verbose == YES) PhyML_Printf("\n. Edge %3d %13G %f %13G",
3427	tree->n_root->b[1]->num,tree->n_root->b[1]->l->v,tree->c_lnL,
3428	tree->n_root->b[1]->l_var->v
3429	);
3430
3431	Lk(tree->n_root->b[2],tree);
3432	if(verbose == YES) PhyML_Printf("\n. Edge %3d %13G %f %13G",
3433	tree->n_root->b[2]->num,tree->n_root->b[2]->l->v,tree->c_lnL,
3434	tree->n_root->b[2]->l_var->v
3435	);
3436
3437	}
3438
3439	res=1;
3440	for(i=1;i<2*tree->n_otu-3;i++)
3441	{
3442	if(FABS(lk[i]-lk[i-1]) > 1.E-2) res=0;
3443	}
3444	Free(lk);
3445
3446	return res;
3447	}
3448
3449	//////////////////////////////////////////////////////////////
3450	//////////////////////////////////////////////////////////////
3451
3452	void ML_Ancestral_Sequences(t_tree *tree)
3453	{
3454	int i;
3455
3456	For(i,2*tree->n_otu-1)
3457	if(tree->a_nodes[i]->tax == NO)
3458	ML_Ancestral_Sequences_One_Node(tree->a_nodes[i],tree);
3459
3460	}
3461
3462	//////////////////////////////////////////////////////////////
3463	//////////////////////////////////////////////////////////////
3464
3465	void ML_Ancestral_Sequences_One_Node(t_node mixt_d, t_tree mixt_tree)
3466	{
3467	if(mixt_d->tax) return;
3468	else
3469	{
3470	t_node v0,v1,*v2; // three neighbours of d
3471	t_edge b0,b1,*b2;
3472	int i,j;
3473	int catg;
3474	phydbl p0, p1, p2;
3475	phydbl *p;
3476	t_node d,curr_mixt_d;
3477	t_tree tree, curr_mixt_tree;
3478	int site;
3479	phydbl p_lk0, p_lk1, *p_lk2;
3480	int sum_scale0, sum_scale1, *sum_scale2;
3481	phydbl r_mat_weight_sum, e_frq_weight_sum, sum_probas;
3482	phydbl Pij0, Pij1, *Pij2;
3483	int NsNs, Ns, NsNg;
3484
3485	curr_mixt_tree = mixt_tree;
3486	curr_mixt_d = mixt_d;
3487
3488	do // For each partition element
3489	{
3490
3491	if(curr_mixt_tree->next)
3492	{
3493	r_mat_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(curr_mixt_tree->next->mod->r_mat_weight);
3494	e_frq_weight_sum = MIXT_Get_Sum_Chained_Scalar_Dbl(curr_mixt_tree->next->mod->e_frq_weight);
3495	sum_probas = MIXT_Get_Sum_Of_Probas_Across_Mixtures(r_mat_weight_sum, e_frq_weight_sum, curr_mixt_tree);
3496	}
3497	else
3498	{
3499	r_mat_weight_sum = 1.;
3500	e_frq_weight_sum = 1.;
3501	sum_probas = 1.;
3502	}
3503
3504	Ns = curr_mixt_tree->next ? curr_mixt_tree->next->mod->ns : curr_mixt_tree->mod->ns;
3505	NsNs = Ns*Ns;
3506	NsNg = Ns*curr_mixt_tree->mod->ras->n_catg;
3507
3508	p = (phydbl *)mCalloc(Ns,sizeof(phydbl));
3509
3510	For(site,curr_mixt_tree->n_pattern) // For each site in the current partition element
3511	{
3512
3513	d = curr_mixt_d->next ? curr_mixt_d->next : curr_mixt_d;
3514	tree = curr_mixt_tree->next ? curr_mixt_tree->next : curr_mixt_tree;
3515
3516	For(i,tree->mod->ns) p[i] = .0;
3517
3518	do // For each class of the mixture model that applies to the current partition element
3519	{
3520	if(tree->is_mixt_tree == YES)
3521	{
3522	tree = tree->next;
3523	d = d->next;
3524	}
3525
3526	v0 = d->v[0];
3527	v1 = d->v[1];
3528	v2 = d->v[2];
3529
3530	b0 = d->b[0];
3531	b1 = d->b[1];
3532	b2 = d->b[2];
3533
3534	Pij0 = b0->Pij_rr;
3535	Pij1 = b1->Pij_rr;
3536	Pij2 = b2->Pij_rr;
3537
3538	if(v0 == b0->left)
3539	{
3540	p_lk0 = b0->p_lk_left;
3541	sum_scale0 = b0->sum_scale_left;
3542	}
3543	else
3544	{
3545	p_lk0 = b0->p_lk_rght;
3546	sum_scale0 = b0->sum_scale_rght;
3547	}
3548
3549	if(v1 == b1->left)
3550	{
3551	p_lk1 = b1->p_lk_left;
3552	sum_scale1 = b1->sum_scale_left;
3553	}
3554	else
3555	{
3556	p_lk1 = b1->p_lk_rght;
3557	sum_scale1 = b1->sum_scale_rght;
3558	}
3559
3560	if(v2 == b2->left)
3561	{
3562	p_lk2 = b2->p_lk_left;
3563	sum_scale2 = b2->sum_scale_left;
3564	}
3565	else
3566	{
3567	p_lk2 = b2->p_lk_rght;
3568	sum_scale2 = b2->sum_scale_rght;
3569	}
3570
3571
3572	For(catg,tree->mod->ras->n_catg)
3573	{
3574	For(i,Ns)
3575	{
3576	p0 = .0;
3577	if(v0->tax)
3578	For(j,tree->mod->ns)
3579	{
3580	p0 += v0->b[0]->p_lk_tip_r[siteNs+j] Pij0[catgNsNs+iNs+j];
3581
3582	/* printf("\n. p0 %d %f", */
3583	/* v0->b[0]->p_lk_tip_r[siteNs+j], /
3584	/* Pij0[catgNsNs+iNs+j]); */
3585	}
3586	else
3587	For(j,tree->mod->ns)
3588	{
3589	p0 += p_lk0[siteNsNg+catgNs+j] * Pij0[catgNsNs+iNs+j] / (phydbl)POW(2,sum_scale0[catg*curr_mixt_tree->n_pattern+site]);
3590
3591	/* p0 += p_lk0[siteNsNg+catgNs+j] * Pij0[catgNsNs+iNs+j]; */
3592
3593	/* printf("\n. p0 %f %f", */
3594	/* p_lk0[siteNsNg+catgNs+j], */
3595	/* Pij0[catgNsNs+iNs+j]); */
3596	}
3597	p1 = .0;
3598	if(v1->tax)
3599	For(j,tree->mod->ns)
3600	{
3601	p1 += v1->b[0]->p_lk_tip_r[siteNs+j] Pij1[catgNsNs+iNs+j];
3602
3603	/* printf("\n. p1 %d %f", */
3604	/* v1->b[0]->p_lk_tip_r[siteNs+j], /
3605	/* Pij1[catgNsNs+iNs+j]); */
3606	}
3607
3608	else
3609	For(j,tree->mod->ns)
3610	{
3611	p1 += p_lk1[siteNsNg+catgNs+j] * Pij1[catgNsNs+iNs+j] / (phydbl)POW(2,sum_scale1[catg*curr_mixt_tree->n_pattern+site]);
3612
3613	/* p1 += p_lk1[siteNsNg+catgNs+j] * Pij1[catgNsNs+iNs+j]; */
3614
3615	/* printf("\n. p1 %f %f", */
3616	/* p_lk1[siteNsNg+catgNs+j], */
3617	/* Pij1[catgNsNs+iNs+j]); */
3618	}
3619
3620
3621	p2 = .0;
3622	if(v2->tax)
3623	For(j,tree->mod->ns)
3624	{
3625	p2 += v2->b[0]->p_lk_tip_r[siteNs+j] Pij2[catgNsNs+iNs+j];
3626	/* printf("\n. p2 %d %f", */
3627	/* v2->b[0]->p_lk_tip_r[siteNs+j], /
3628	/* Pij2[catgNsNs+iNs+j]); */
3629	}
3630	else
3631	For(j,tree->mod->ns)
3632	{
3633	p2 += p_lk2[siteNsNg+catgNs+j] * Pij2[catgNsNs+iNs+j] / (phydbl)POW(2,sum_scale2[catg*curr_mixt_tree->n_pattern+site]);
3634
3635	/* p2 += p_lk2[siteNsNg+catgNs+j] * Pij2[catgNsNs+iNs+j]; */
3636
3637	/* printf("\n. p2 %f %f", */
3638	/* p_lk2[siteNsNg+catgNs+j], */
3639	/* Pij2[catgNsNs+iNs+j]); */
3640	}
3641
3642	p[i] +=
3643	p0p1p2*
3644	tree->mod->e_frq->pi->v[i] /
3645	tree->cur_site_lk[site] *
3646	curr_mixt_tree->mod->ras->gamma_r_proba->v[tree->mod->ras->parent_class_number] *
3647	tree->mod->r_mat_weight->v / r_mat_weight_sum *
3648	tree->mod->e_frq_weight->v / e_frq_weight_sum /
3649	sum_probas;
3650
3651	}
3652	}
3653
3654	For(i,Ns)
3655	{
3656	printf("\n. p%d: %f",i,p[i]);
3657	}
3658	Exit("\n");
3659	}
3660	while(tree->is_mixt_tree == NO);
3661	}
3662
3663	Free(p);
3664	curr_mixt_tree = curr_mixt_tree->next_mixt;
3665	curr_mixt_d = curr_mixt_d->next_mixt;
3666	}
3667	while(curr_mixt_tree != NULL);
3668	}
3669	}
3670
3671	//////////////////////////////////////////////////////////////
3672	//////////////////////////////////////////////////////////////
3673
3674	//////////////////////////////////////////////////////////////
3675	//////////////////////////////////////////////////////////////
3676
3677	//////////////////////////////////////////////////////////////
3678	//////////////////////////////////////////////////////////////
3679
3680	//////////////////////////////////////////////////////////////
3681	//////////////////////////////////////////////////////////////
3682

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

Context Navigation

source: branches/ali/GDE/PHYML20130708/phyml/src/lk.c

Download in other formats: