Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Visit:

Last change on this file was 4073, checked in by westram, 18 years ago
phyml 2.4.5
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 111.3 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 "utilities.h"
14	#include "lk.h"
15	#include "optimiz.h"
16	#include "models.h"
17	#include "free.h"
18	#include "bionj.h"
19	#include "simu.h"
20
21
22	/*********************************************************/
23	/* NUMERICAL RECIPES ROUTINES FOR COMPUTING C(n,k) */
24	/*********************************************************/
25
26	double bico(int n, int k)
27	{
28	return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
29	}
30
31	double factln(int n)
32	{
33	static double a[101];
34
35	if (n < 0){ Exit("Err: negative factorial in routine FACTLN"); }
36	if (n <= 1) return 0.0;
37	if (n <= 100) return a[n] ? a[n] : (a[n]=gammln(n+1.0));
38	else return gammln(n+1.0);
39	}
40
41	double gammln(double xx)
42	{
43	double x,tmp,ser;
44	static double cof[6]={76.18009173,-86.50532033,24.01409822,
45	-1.231739516,0.120858003e-2,-0.536382e-5};
46	int j;
47
48	x=xx-1.0;
49	tmp=x+5.5;
50	tmp -= (x+0.5)*log(tmp);
51	ser=1.0;
52	for (j=0;j<=5;j++) {
53	x += 1.0;
54	ser += cof[j]/x;
55	}
56	return -tmp+log(2.50662827465*ser);
57	}
58
59	/*********************************************************/
60	/* END OF NUMERICAL RECIPES ROUTINES */
61	/*********************************************************/
62
63	double Pbinom(int N, int ni, double p)
64	{
65	return bico(N,ni)pow(p,ni)pow(1-p,N-ni);
66	}
67
68	/*********************************************************/
69
70	void Plim_Binom(double pH0, int N, double pinf, double psup)
71	{
72	pinf = pH0 - 1.64sqrt(pH0*(1-pH0)/(double)N);
73	if(pinf < 0) pinf = .0;
74	psup = pH0 + 1.64sqrt(pH0*(1-pH0)/(double)N);
75	}
76
77	/*********************************************************/
78
79	double LnGamma (double alpha)
80	{
81	/* returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.
82	Stirling's formula is used for the central polynomial part of the procedure.
83	Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
84	Communications of the Association for Computing Machinery, 9:684
85	*/
86	double x=alpha, f=0, z;
87
88	if (x<7) {
89	f=1; z=x-1;
90	while (++z<7) f*=z;
91	x=z; f=-log(f);
92	}
93	z = 1/(x*x);
94	return f + (x-0.5)*log(x) - x + .918938533204673
95	+ (((-.000595238095238z+.000793650793651)z-.002777777777778)*z
96	+.083333333333333)/x;
97	}
98
99	/*********************************************************/
100
101	double IncompleteGamma (double x, double alpha, double ln_gamma_alpha)
102	{
103	/* returns the incomplete gamma ratio I(x,alpha) where x is the upper
104	limit of the integration and alpha is the shape parameter.
105	returns (-1) if in error
106	ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
107	(1) series expansion if (alpha>x \|\| x<=1)
108	(2) continued fraction otherwise
109	RATNEST FORTRAN by
110	Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics,
111	19: 285-287 (AS32)
112	*/
113	int i;
114	double p=alpha, g=ln_gamma_alpha;
115	double accurate=1e-8, overflow=1e30;
116	double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, pn[6];
117
118	if (x==0) return (0);
119	if (x<0 \|\| p<=0) return (-1);
120
121	factor=exp(p*log(x)-x-g);
122	if (x>1 && x>=p) goto l30;
123	/* (1) series expansion */
124	gin=1; term=1; rn=p;
125	l20:
126	rn++;
127	term*=x/rn; gin+=term;
128
129	if (term > accurate) goto l20;
130	gin*=factor/p;
131	goto l50;
132	l30:
133	/* (2) continued fraction */
134	a=1-p; b=a+x+1; term=0;
135	pn[0]=1; pn[1]=x; pn[2]=x+1; pn[3]=x*b;
136	gin=pn[2]/pn[3];
137	l32:
138	a++; b+=2; term++; an=a*term;
139	for (i=0; i<2; i++) pn[i+4]=bpn[i+2]-anpn[i];
140	if (pn[5] == 0) goto l35;
141	rn=pn[4]/pn[5]; dif=fabs(gin-rn);
142	if (dif>accurate) goto l34;
143	if (dif<=accurate*rn) goto l42;
144	l34:
145	gin=rn;
146	l35:
147	for (i=0; i<4; i++) pn[i]=pn[i+2];
148	if (fabs(pn[4]) < overflow) goto l32;
149	for (i=0; i<4; i++) pn[i]/=overflow;
150	goto l32;
151	l42:
152	gin=1-factor*gin;
153
154	l50:
155	return (gin);
156	}
157
158
159	/*********************************************************/
160
161	double PointChi2 (double prob, double v)
162	{
163	/* returns z so that Prob{x<z}=prob where x is Chi2 distributed with df=v
164	returns -1 if in error. 0.000002<prob<0.999998
165	RATNEST FORTRAN by
166	Best DJ & Roberts DE (1975) The percentage points of the
167	Chi2 distribution. Applied Statistics 24: 385-388. (AS91)
168	Converted into C by Ziheng Yang, Oct. 1993.
169	*/
170	double e=.5e-6, aa=.6931471805, p=prob, g;
171	double xx, c, ch, a=0,q=0,p1=0,p2=0,t=0,x=0,b=0,s1,s2,s3,s4,s5,s6;
172
173	if (p<.000002 \|\| p>.999998 \|\| v<=0) return (-1);
174
175	g = LnGamma (v/2);
176	xx=v/2; c=xx-1;
177	if (v >= -1.24*log(p)) goto l1;
178
179	ch=pow((pxxexp(g+xx*aa)), 1/xx);
180	if (ch-e<0) return (ch);
181	goto l4;
182	l1:
183	if (v>.32) goto l3;
184	ch=0.4; a=log(1-p);
185	l2:
186	q=ch; p1=1+ch(4.67+ch); p2=ch(6.73+ch*(6.66+ch));
187	t=-0.5+(4.67+2ch)/p1 - (6.73+ch(13.32+3*ch))/p2;
188	ch-=(1-exp(a+g+.5ch+caa)*p2/p1)/t;
189	if (fabs(q/ch-1)-.01 <= 0) goto l4;
190	else goto l2;
191
192	l3:
193	x=PointNormal (p);
194	p1=0.222222/v; ch=vpow((xsqrt(p1)+1-p1), 3.0);
195	if (ch>2.2v+6) ch=-2(log(1-p)-clog(.5ch)+g);
196	l4:
197	q=ch; p1=.5*ch;
198	if ((t=IncompleteGamma (p1, xx, g))<0) {
199	printf ("\nerr IncompleteGamma");
200	return (-1);
201	}
202	p2=p-t;
203	t=p2exp(xxaa+g+p1-c*log(ch));
204	b=t/ch; a=0.5t-bc;
205
206	s1=(210+a(140+a(105+a(84+a(70+60*a))))) / 420;
207	s2=(420+a(735+a(966+a(1141+1278a))))/2520;
208	s3=(210+a(462+a(707+932*a)))/2520;
209	s4=(252+a(672+1182a)+c(294+a(889+1740*a)))/5040;
210	s5=(84+264a+c(175+606*a))/2520;
211	s6=(120+c(346+127c))/5040;
212	ch+=t(1+0.5ts1-bc(s1-b(s2-b(s3-b(s4-b(s5-bs6))))));
213	if (fabs(q/ch-1) > e) goto l4;
214
215	return (ch);
216	}
217
218	/*********************************************************/
219
220	double PointNormal (double prob)
221	{
222	/* returns z so that Prob{x<z}=prob where x ~ N(0,1) and (1e-12)<prob<1-(1e-12)
223	returns (-9999) if in error
224	Odeh RE & Evans JO (1974) The percentage points of the normal distribution.
225	Applied Statistics 22: 96-97 (AS70)
226
227	Newer methods:
228	Wichura MJ (1988) Algorithm AS 241: the percentage points of the
229	normal distribution. 37: 477-484.
230	Beasley JD & Springer SG (1977). Algorithm AS 111: the percentage
231	points of the normal distribution. 26: 118-121.
232
233	*/
234	double a0=-.322232431088, a1=-1, a2=-.342242088547, a3=-.0204231210245;
235	double a4=-.453642210148e-4, b0=.0993484626060, b1=.588581570495;
236	double b2=.531103462366, b3=.103537752850, b4=.0038560700634;
237	double y, z=0, p=prob, p1;
238
239	p1 = (p<0.5 ? p : 1-p);
240	if (p1<1e-20) return (-9999);
241
242	y = sqrt (log(1/(p1*p1)));
243	z = y + ((((ya4+a3)y+a2)y+a1)y+a0) / ((((yb4+b3)y+b2)y+b1)y+b0);
244	return (p<0.5 ? -z : z);
245	}
246	/*********************************************************/
247
248	int DiscreteGamma (double freqK[], double rK[],
249	double alfa, double beta, int K, int median)
250	{
251	/* discretization of gamma distribution with equal proportions in each
252	category
253	*/
254	int i;
255	double gap05=1.0/(2.0K), t, factor=alfa/betaK, lnga1;
256
257	if(K==1)
258	{
259	rK[0] = 1.0;
260	return 0;
261	}
262
263	if (median) {
264	for (i=0; i<K; i++) rK[i]=PointGamma((i2.0+1)gap05, alfa, beta);
265	for (i=0,t=0; i<K; i++) t+=rK[i];
266	for (i=0; i<K; i++) rK[i]*=factor/t;
267	}
268	else {
269	lnga1=LnGamma(alfa+1);
270	for (i=0; i<K-1; i++)
271	freqK[i]=PointGamma((i+1.0)/K, alfa, beta);
272	for (i=0; i<K-1; i++)
273	freqK[i]=IncompleteGamma(freqK[i]*beta, alfa+1, lnga1);
274	rK[0] = freqK[0]*factor;
275	rK[K-1] = (1-freqK[K-2])*factor;
276	for (i=1; i<K-1; i++) rK[i] = (freqK[i]-freqK[i-1])*factor;
277	}
278	for (i=0; i<K; i++) freqK[i]=1.0/K;
279
280	return (0);
281	}
282
283	/*********************************************************/
284
285	arbre Read_Tree(char s_tree)
286	{
287	char **subs;
288	int i,n_ext,n_int,n_otu;
289	char *sub_tp;
290	arbre *tree;
291	int degree;
292
293
294	tree=(arbre *)mCalloc(1,sizeof(arbre));
295
296
297	n_otu=0;
298	For(i,(int)strlen(s_tree)) if(s_tree[i] == ',') n_otu++;
299	n_otu+=1;
300	tree->n_otu = n_otu;
301
302
303	Init_Tree(tree);
304
305
306	tree->noeud[n_otu]=(node *)mCalloc(1,sizeof(node));
307	tree->noeud[n_otu]->v = NULL;
308	Make_Node_Light(tree->noeud[n_otu]);
309	tree->noeud[n_otu]->num=n_otu;
310	tree->noeud[n_otu]->tax=0;
311
312
313
314	subs = Sub_Trees(s_tree,&degree);
315	Clean_Multifurcation(subs,degree,3);
316
317
318	if(!strlen(subs[2])) Exit("\n. Err: unrooted tree is needed\n");
319	sub_tp=(char *)mCalloc((int)strlen(s_tree),sizeof(char));
320
321
322	tree->has_branch_lengths = 1;
323
324	For(i,degree)
325	{
326	strcpy(sub_tp,subs[i]);
327	strcat(sub_tp,":");
328	if(strstr(s_tree,sub_tp))
329	{
330	tree->noeud[n_otu]->l[i] =
331	atof((char *)strstr(s_tree,sub_tp)+strlen(subs[i])+1);
332	}
333	else
334	tree->has_branch_lengths = 0;
335	}
336	Free(sub_tp);
337
338	n_int = n_ext = 0;
339	For(i,degree)
340	R_rtree(subs[i],tree->noeud[n_otu],tree,&n_int,&n_ext);
341
342	Make_All_Edges_Light(tree->noeud[0],tree->noeud[0]->v[0]);
343	i = 0;
344	Init_Tree_Edges(tree->noeud[0],tree->noeud[0]->v[0],tree,&i);
345
346	For(i,NODE_DEG_MAX)
347	Free(subs[i]);
348	Free(subs);
349
350	return tree;
351	}
352
353	/*********************************************************/
354
355	void Make_All_Edges_Light(node a, node d)
356	{
357	int i;
358
359	Make_Edge_Light(a,d);
360	if(d->tax) return;
361	else
362	{
363	For(i,3)
364	{
365	if(d->v[i] != a)
366	Make_All_Edges_Light(d,d->v[i]);
367	}
368	}
369	}
370
371	/*********************************************************/
372
373	void Make_All_Edges_Lk(node a, node d, arbre *tree)
374	{
375	int i;
376
377	Make_Edge_Lk(a,d,tree);
378	if(d->tax) return;
379	else
380	{
381	For(i,3)
382	{
383	if(d->v[i] != a)
384	Make_All_Edges_Lk(d,d->v[i],tree);
385	}
386	}
387	}
388
389	/*********************************************************/
390
391	void R_rtree(char s_tree, node pere, arbre tree, int n_int, int *n_ext)
392	{
393	int i;
394	char *sub_tp;
395	node *fils;
396	int n_otu = tree->n_otu;
397
398	if(strstr(s_tree," ")) Exit("\n Err : tree must not contain a ' ' character\n");
399
400	fils=(node *)mCalloc(1,sizeof(node));
401	fils->v = NULL;
402	Make_Node_Light(fils);
403
404	if(s_tree[0] == '(')
405	{
406	char **subs;
407	int degree;
408
409
410	(*n_int)+=1;
411	tree->noeud[n_otu+*n_int]=fils;
412	fils->num=n_otu+*n_int;
413	fils->tax=0;
414
415	if(s_tree[(int)strlen(s_tree)-1] == '*')
416	{
417	fils->check_branch = 1;
418	s_tree[(int)strlen(s_tree)-1] = '\0';
419	}
420
421
422	For(i,3)
423	{
424	if(!pere->v[i])
425	{
426	pere->v[i]=fils;
427	fils->l[0]=pere->l[i];
428	break;
429	}
430	}
431
432	fils->v[0]=pere;
433	subs=Sub_Trees(s_tree,&degree);
434
435	Clean_Multifurcation(subs,degree,2);
436
437	sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char));
438	strcpy(sub_tp,subs[0]);
439	strcat(sub_tp,":");
440	if(strstr(s_tree,sub_tp))
441	{
442	fils->l[1] = atof((char *)strstr(s_tree,sub_tp)
443	+(int)strlen(subs[0])+1);
444	}
445
446	strcpy(sub_tp,subs[1]);
447	strcat(sub_tp,":");
448	if(strstr(s_tree,sub_tp))
449	{
450	fils->l[2] = atof((char *)strstr(s_tree,sub_tp)
451	+(int)strlen(subs[1])+1);
452	}
453
454	Free(sub_tp);
455	R_rtree(subs[0],fils,tree,n_int,n_ext);
456	R_rtree(subs[1],fils,tree,n_int,n_ext);
457	For(i,NODE_DEG_MAX) Free(subs[i]);
458	Free(subs);
459	}
460
461	else
462	{
463	tree->noeud[*n_ext]=fils;
464	fils->tax=1;
465	For(i,3)
466	{
467	if(!pere->v[i])
468	{
469	pere->v[i]=fils;
470	fils->l[0]=pere->l[i];
471	break;
472	}
473	}
474
475	if(s_tree[(int)strlen(s_tree)-1] == '*')
476	{
477	fils->check_branch = 1;
478	s_tree[(int)strlen(s_tree)-1] = '\0';
479	}
480
481
482	fils->v[0]=pere;
483	strcpy(fils->name,s_tree);
484	fils->num=*n_ext;
485	(*n_ext)+=1;
486	}
487	}
488
489	/*********************************************************/
490
491	void Clean_Multifurcation(char **subtrees, int current_deg, int end_deg)
492	{
493
494	if(current_deg <= end_deg) return;
495	else
496	{
497	char *s_tmp;
498	int i;
499
500	s_tmp = (char *)mCalloc(T_MAX_LINE,sizeof(char));
501
502	strcat(s_tmp,"(\0");
503	strcat(s_tmp,subtrees[0]);
504	strcat(s_tmp,",\0");
505	strcat(s_tmp,subtrees[1]);
506	strcat(s_tmp,")\0");
507	Free(subtrees[0]);
508	subtrees[0] = s_tmp;
509
510
511	for(i=1;i<current_deg-1;i++)
512	{
513	strcpy(subtrees[i],subtrees[i+1]);
514	}
515
516	Clean_Multifurcation(subtrees,current_deg-1,end_deg);
517	}
518	}
519
520	/*********************************************************/
521
522	char *Sub_Trees(char tree, int *degree)
523	{
524	char **subs;
525	int posbeg,posend;
526	int i;
527
528	subs=(char *)mCalloc(NODE_DEG_MAX,sizeof(char ));
529
530	For(i,NODE_DEG_MAX) subs[i]=(char *)mCalloc(strlen(tree)+1,sizeof(char));
531
532	posbeg=posend=1;
533	(*degree)=0;
534	do
535	{
536	posbeg = posend;
537	if(tree[posend] != '(')
538	{
539	while((tree[posend] != ',' ) &&
540	(tree[posend] != ':' ) &&
541	(tree[posend] != ')' ))
542	posend += 1;
543	posend -= 1;
544	}
545	else posend=Next_Par(tree,posend);
546
547	while((tree[posend+1] != ',') &&
548	(tree[posend+1] != ':') &&
549	(tree[posend+1] != ')')) {posend++;}
550
551
552	strncpy(subs[(*degree)],tree+posbeg,posend-posbeg+1);
553	strcat(subs[(*degree)],"\0");
554
555	posend += 1;
556	while((tree[posend] != ',') &&
557	(tree[posend] != ')')) {posend++;}
558	posend+=1;
559
560
561	(*degree)++;
562	if((*degree) == NODE_DEG_MAX)
563	{
564	For(i,(*degree))
565	printf("\n. Subtree %d : %s\n",i+1,subs[i]);
566
567	printf("\n. The degree of a node cannot be greater than %d\n",NODE_DEG_MAX);
568	Exit("\n");
569	}
570	}
571	while(tree[posend-1] != ')');
572
573	return subs;
574	}
575
576
577	/*********************************************************/
578
579	int Next_Par(char *s, int pos)
580	{
581	int curr;
582
583	curr=pos+1;
584
585	while(*(s+curr) != ')')
586	{
587	if(*(s+curr) == '(') curr=Next_Par(s,curr);
588	curr++;
589	}
590
591	return curr;
592	}
593
594	/*********************************************************/
595
596	char Write_Tree(arbre tree)
597	{
598
599	char *s;
600	int i;
601
602	s=(char *)mCalloc(T_MAX_LINE,sizeof(char));
603
604	s[0]='(';
605
606	i = 0;
607	while((!tree->noeud[tree->n_otu+i]->v[0]) \|\|
608	(!tree->noeud[tree->n_otu+i]->v[1]) \|\|
609	(!tree->noeud[tree->n_otu+i]->v[2])) i++;
610
611	R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[0],s,tree);
612	R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[1],s,tree);
613	R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[2],s,tree);
614
615	s[(int)strlen(s)-1]=')';
616	s[(int)strlen(s)]=';';
617
618
619	return s;
620	}
621
622	/*********************************************************/
623
624	void R_wtree(node pere, node fils, char s_tree, arbre tree)
625	{
626	int i,p;
627
628	p = -1;
629	if(fils->tax)
630	{
631
632	strcat(s_tree,fils->name);
633	if((fils->b[0]) && (fils->b[0]->l != -1))
634	{
635	strcat(s_tree,":");
636	sprintf(s_tree+(int)strlen(s_tree),"%f",fils->b[0]->l);
637	fflush(stdout);
638	}
639	sprintf(s_tree+(int)strlen(s_tree),",");
640	}
641
642	else
643	{
644	s_tree[(int)strlen(s_tree)]='(';
645	For(i,3)
646	{
647	if(fils->v[i] != pere)
648	R_wtree(fils,fils->v[i],s_tree,tree);
649	else p=i;
650	}
651	s_tree[(int)strlen(s_tree)-1]=')';
652	if(fils->b[0]->l != -1)
653	{
654	if(tree->print_boot_val)
655	sprintf(s_tree+(int)strlen(s_tree),"%d",fils->b[p]->bip_score);
656	strcat(s_tree,":");
657	sprintf(s_tree+(int)strlen(s_tree),"%f,",fils->b[p]->l);
658	fflush(stdout);
659	}
660	}
661	}
662
663	/*********************************************************/
664
665	void Init_Tree(arbre *tree)
666	{
667	int i;
668
669	tree->noeud = (node *)mCalloc(2tree->n_otu-2,sizeof(node *));
670	tree->t_edges = (edge *)mCalloc(2tree->n_otu-3,sizeof(edge *));
671	For(i,2*tree->n_otu-3)
672	tree->t_edges[i] = NULL;
673
674	tree->best_tree = NULL;
675	tree->old_tree = NULL;
676
677	tree->has_bip = 0;
678
679	tree->best_loglk = UNLIKELY;
680	tree->tot_loglk = UNLIKELY;
681	tree->n_swap = 0;
682
683	tree->tot_dloglk = NULL;
684	tree->tot_d2loglk = NULL;
685
686	tree->n_pattern = -1;
687
688	tree->print_boot_val = 0;
689	}
690
691	/*********************************************************/
692
693	void Make_Edge_Light(node a, node d)
694	{
695	edge *b;
696
697	b = (edge *)mCalloc(1,sizeof(edge));
698
699	Init_Edge_Light(b);
700
701
702	b->left = a; b->rght = d;
703	if(a->tax) {b->rght = a; b->left = d;} /* root */
704	/* a tip is necessary on the right side of the edge */
705
706	Make_Edge_Dirs(b,a,d);
707
708	b->l = a->l[b->l_r];
709	if(a->tax) b->l = a->l[b->r_l];
710	if(b->l < BL_MIN) b->l = BL_MIN;
711	b->l_old = b->l;
712	}
713
714	/*********************************************************/
715
716	void Init_Edge_Light(edge *b)
717	{
718	b->bip_score = 0;
719	b->nj_score = .0;
720
721	b->site_dlk = -1.;
722	b->site_d2lk = -1.;
723	b->site_dlk_rr = NULL;
724	b->site_d2lk_rr = NULL;
725	b->p_lk_left = NULL;
726	b->p_lk_rght = NULL;
727	b->Pij_rr = NULL;
728	b->dPij_rr = NULL;
729	b->d2Pij_rr = NULL;
730
731	}
732
733	/*********************************************************/
734
735	void Make_Edge_Dirs(edge b, node a, node *d)
736	{
737	int i;
738
739	b->l_r = b->r_l = -1;
740	For(i,3)
741	{
742	if((a->v[i]) && (a->v[i] == d))
743	{
744	b->l_r = i;
745	a->b[i] = b;
746	}
747	if((d->v[i]) && (d->v[i] == a))
748	{
749	b->r_l = i;
750	d->b[i] = b;
751	}
752	}
753
754	if(a->tax) {b->r_l = 0; For(i,3) if(d->v[i]==a) {b->l_r = i; break;}}
755
756
757	b->l_v1 = b->l_v2 = b->r_v1 = b->r_v2 = -1;
758	For(i,3)
759	{
760	if(b->left->v[i] != b->rght)
761	{
762	if(b->l_v1 < 0) b->l_v1 = i;
763	else b->l_v2 = i;
764	}
765
766	if(b->rght->v[i] != b->left)
767	{
768	if(b->r_v1 < 0) b->r_v1 = i;
769	else b->r_v2 = i;
770	}
771	}
772	}
773
774	/*********************************************************/
775
776	void Make_Edge_Lk(node a, node d, arbre *tree)
777	{
778	int i,j;
779	int len;
780	edge *b;
781
782	b = NULL;
783
784	For(i,3) if((a->v[i]) && (a->v[i] == d)) {b = a->b[i]; break;}
785
786	len = (int)tree->data->crunch_len;
787
788	b->diff_lk = 0.0;
789	b->l_old = b->l;
790
791	if(!b->Pij_rr)
792	{
793	b->best_conf = 1;
794	b->ql = (double *)mCalloc(3,sizeof(double));
795
796
797	b->Pij_rr = (double *)mCalloc(tree->mod->n_catg,sizeof(double ));
798	b->dPij_rr = (double *)mCalloc(tree->mod->n_catg,sizeof(double ));
799	b->d2Pij_rr = (double *)mCalloc(tree->mod->n_catg,sizeof(double ));
800
801	For(i,tree->mod->n_catg)
802	{
803	b->Pij_rr[i] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
804	b->dPij_rr[i] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
805	b->d2Pij_rr[i] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
806
807	For(j,tree->mod->ns)
808	{
809	b->Pij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
810	b->dPij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
811	b->d2Pij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double ));
812	}
813	}
814
815	b->site_p_lk_left = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double ));
816	For(i,tree->mod->n_catg) b->site_p_lk_left[i] = (double *)mCalloc(tree->mod->ns,sizeof(double));
817	b->site_p_lk_rght = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double ));
818	For(i,tree->mod->n_catg) b->site_p_lk_rght[i] = (double *)mCalloc(tree->mod->ns,sizeof(double));
819
820	b->site_dlk_rr = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double));
821	b->site_d2lk_rr = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double));
822
823	b->scale_left = b->scale_rght = 0;
824
825	b->p_lk_left = NULL;
826	b->p_lk_rght = NULL;
827
828	b->sum_scale_f_rght = NULL;
829	b->sum_scale_f_left = NULL;
830
831	b->get_p_lk_left = 0;
832	b->get_p_lk_rght = 0;
833	b->ud_p_lk_left = 0;
834	b->ud_p_lk_rght = 0;
835	}
836	}
837
838	/*********************************************************/
839
840	void Make_Node_Light(node *n)
841	{
842
843	if(n->v) return;
844
845	n->v = (node *)mCalloc(3,sizeof(node ));
846	n->l = (double *)mCalloc(3,sizeof(double));
847	n->b = (edge *)mCalloc(3,sizeof(edge ));
848	n->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
849	n->score = (double *)mCalloc(3,sizeof(double));
850
851	Init_Node_Light(n);
852
853	}
854
855	/*********************************************************/
856
857	void Init_Node_Light(node *n)
858	{
859	int i;
860
861	n->check_branch = 0;
862	/* n->is_attach = 0; */
863	/* n->is_free = 0; */
864
865	For(i,3)
866	{
867	n->v[i]=NULL;
868	n->b[i]=NULL;
869	n->l[i]=-1;
870	}
871	}
872
873	/*********************************************************/
874
875	void Make_Node_Lk(node *n)
876	{
877	/* n->n_ex_nodes = (int )mCalloc(2,sizeof(int)); /
878	return;
879	}
880
881	/*********************************************************/
882
883	seq *Get_Seq(option input, int rw)
884	{
885	seq **data;
886	int i,j;
887	char **buff;
888	int n_unkn,n_removed,pos;
889	int *remove;
890
891
892	/* rewind(fp_seq); */
893
894	if(input->interleaved) data = Read_Seq_Interleaved(input->fp_seq,&(input->mod->n_otu));
895	else data = Read_Seq_Sequential(input->fp_seq,&(input->mod->n_otu));
896
897	if(data)
898	{
899	buff = (char *)mCalloc(input->mod->n_otu,sizeof(char ));
900	For(i,input->mod->n_otu) buff[i] = (char *)mCalloc(data[0]->len,sizeof(char));
901	remove = (int *)mCalloc(data[0]->len,sizeof(int));
902
903	n_removed = 0;
904
905	For(i,data[0]->len)
906	{
907	For(j,input->mod->n_otu)
908	{
909	if((data[j]->state[i] == '?') \|\|
910	(data[j]->state[i] == '-')) data[j]->state[i] = 'X';
911
912	if(data[j]->state[i] == 'U') data[j]->state[i] = 'T';
913
914	if((input->mod->datatype == NT) && (data[j]->state[i] == 'N')) data[j]->state[i] = 'X';
915
916	}
917
918	n_unkn = 0;
919	For(j,input->mod->n_otu) if(data[j]->state[i] == 'X') n_unkn++;
920
921	if(n_unkn == input->mod->n_otu)
922	{
923	remove[i] = 1;
924	n_removed++;
925	}
926
927	For(j,input->mod->n_otu) buff[j][i] = data[j]->state[i];
928	}
929
930	if(n_removed > 0)
931	{
932	if(input->mod->datatype == NT)
933	printf("\n. %d sites are made from completely undetermined states ('X', '-', '?' or 'N')...\n",n_removed);
934	else
935	printf("\n. %d sites are made from completely undetermined states ('X', '-', '?')...\n",n_removed);
936	}
937
938	pos = 0;
939	For(i,data[0]->len)
940	{
941	/* if(!remove[i]) */
942	/* { */
943	For(j,input->mod->n_otu) data[j]->state[pos] = buff[j][i];
944	pos++;
945	/* } */
946	}
947
948	For(i,input->mod->n_otu) data[i]->len = pos;
949	For(i,input->mod->n_otu) Free(buff[i]);
950	Free(buff);
951	Free(remove);
952	}
953	return data;
954	}
955
956	/*********************************************************/
957
958	seq *Read_Seq_Sequential(FILE in, int *n_otu)
959	{
960	int i;
961	char *line;
962	int len,readok;
963	seq **data;
964	char c;
965	char format = (char )mCalloc(20, sizeof(char));
966
967	line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
968
969	readok = len = 0;
970	do
971	{
972	if(fscanf(in,"%s",line) == EOF)
973	{
974	Free(line); return NULL;
975	}
976	else
977	{
978	if(strcmp(line,"\n") && strcmp(line,"\n") && strcmp(line,"\t"))
979	{
980	*n_otu = atoi(line);
981	data = (seq *)mCalloc(n_otu,sizeof(seq *));
982	if(*n_otu <= 0) Exit("\n. Problem with sequence format\n");
983	fscanf(in,"%s",line);
984	len = atoi(line);
985	if(len <= 0) Exit("\n. Problem with sequence format\n");
986	else readok = 1;
987	}
988	}
989	}while(!readok);
990
991
992	/* while((c=fgetc(in))!='\n'); */
993	while(((c=fgetc(in))!='\n') && (c != ' ') && (c != '\r') && (c != '\t'));
994
995	For(i,*n_otu)
996	{
997	data[i] = (seq *)mCalloc(1,sizeof(seq));
998	data[i]->len = 0;
999	data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1000	data[i]->state = (char *)mCalloc(T_MAX_SEQ,sizeof(char));
1001	sprintf(format, "%%%ds", T_MAX_NAME);
1002	fscanf(in, format, data[i]->name);
1003
1004	while(data[i]->len < len)
1005	Read_One_Line_Seq(&data,i,in);
1006
1007	if(data[i]->len != len)
1008	{
1009	printf("\n. Err: Problem with species %s's sequence (check the format)\n",
1010	data[i]->name);
1011	Exit("");
1012	}
1013	}
1014
1015	/* fgets(line,T_MAX_LINE,in); */
1016	/* inter data sets */
1017
1018	Free(format);
1019	Free(line);
1020	return data;
1021	}
1022
1023	/*********************************************************/
1024
1025	seq *Read_Seq_Interleaved(FILE in, int *n_otu)
1026	{
1027	int i,end,num_block;
1028	char *line;
1029	int len,readok;
1030	seq **data;
1031	char c;
1032	char *format;
1033
1034	line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1035	format = (char *)mCalloc(T_MAX_NAME, sizeof(char));
1036
1037	readok = len = 0;
1038	do
1039	{
1040	if(fscanf(in,"%s",line) == EOF)
1041	{
1042	Free(format);
1043	Free(line); return NULL;
1044	}
1045	else
1046	{
1047	if(strcmp(line,"\n") && strcmp(line,"\r") && strcmp(line,"\t"))
1048	{
1049	*n_otu = atoi(line);
1050	data = (seq *)mCalloc(n_otu,sizeof(seq *));
1051	if(*n_otu <= 0) Exit("\n. Problem with sequence format\n");
1052	fscanf(in,"%s",line);
1053	len = atoi(line);
1054	if(len <= 0) Exit("\n. Problem with sequence format\n");
1055	else readok = 1;
1056	}
1057	}
1058	}while(!readok);
1059
1060
1061	while(((c=fgetc(in))!='\n') && (c != ' ') && (c != '\r') && (c != '\t'));
1062
1063	end = 0;
1064	For(i,*n_otu)
1065	{
1066	data[i] = (seq *)mCalloc(1,sizeof(seq));
1067	data[i]->len = 0;
1068	data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1069	data[i]->state = (char *)mCalloc(T_MAX_SEQ,sizeof(char));
1070	sprintf(format, "%%%ds", T_MAX_NAME);
1071	fscanf(in, format, data[i]->name);
1072	if(!Read_One_Line_Seq(&data,i,in))
1073	{
1074	end = 1;
1075	if((i != n_otu) && (i != n_otu-1))
1076	{
1077	printf("\n. Err: Problem with species %s's sequence\n",data[i]->name);
1078	Exit("");
1079	}
1080	break;
1081	}
1082	}
1083
1084	if(data[0]->len == len) end = 1;
1085
1086	if(!end)
1087	{
1088	end = 0;
1089
1090	num_block = 1;
1091	do
1092	{
1093	num_block++;
1094
1095	/* interblock */
1096	if(!fgets(line,T_MAX_LINE,in)) break;
1097
1098	if(line[0] != 13 && line[0] != 10)
1099	{
1100	printf("\n. One or more missing sequences in block %d\n",num_block-1);
1101	Exit("");
1102	}
1103
1104	For(i,*n_otu)
1105	if(data[i]->len != len)
1106	break;
1107
1108	if(i == *n_otu) break;
1109
1110
1111	For(i,*n_otu)
1112	{
1113	if(data[i]->len > len)
1114	{
1115	printf("\n. Err: Problem with species %s's sequence\n",data[i]->name);
1116	Exit("");
1117	}
1118	else if(!Read_One_Line_Seq(&data,i,in))
1119	{
1120	end = 1;
1121	if((i != n_otu) && (i != n_otu-1))
1122	{
1123	printf("\n. Err: Problem with species %s's sequence\n",data[i]->name);
1124	Exit("");
1125	}
1126	break;
1127	}
1128	}
1129	}while(!end);
1130	}
1131
1132	For(i,*n_otu)
1133	{
1134	if(data[i]->len != len)
1135	{
1136	printf("\n. Check sequence '%s' length...\n",data[i]->name);
1137	Exit("");
1138	}
1139	}
1140
1141	Free(format);
1142	Free(line);
1143	return data;
1144	}
1145
1146	/*********************************************************/
1147
1148	int Read_One_Line_Seq(seq **data, int num_otu, FILE in)
1149	{
1150	char c;
1151
1152	c=' ';
1153	while(1)
1154	{
1155	/* if((c == EOF) \|\| (c == '\n') \|\| (c == '\r')) break; */
1156	if((c == EOF) \|\| (c == 13) \|\| (c == 10)) break;
1157	else if((c==' ') \|\| (c=='\t')) {c=(char)fgetc(in); continue;}
1158	Uppercase(&c);
1159
1160	/if(strchr("ACGTUMRWSYKBDHVNXO?-.",c) == NULL)/
1161	if (strchr("ABCDEFGHIKLMNOPQRSTUVWXYZ?-.", c) == NULL)
1162	{
1163	printf("\n. Err: bad symbol: \"%c\" at position %d of species %s\n",
1164	c,(data)[num_otu]->len,(data)[num_otu]->name);
1165	Exit("");
1166	}
1167
1168	if(c == '.')
1169	{
1170	c = (data)[0]->state[(data)[num_otu]->len];
1171	if(!num_otu)
1172	Exit("\n. Err: Symbol \".\" should not appear in the first sequence\n");
1173	}
1174	(data)[num_otu]->state[(data)[num_otu]->len]=c;
1175	(*data)[num_otu]->len++;
1176	/* if(c=='U') c='T'; */
1177	c = (char)fgetc(in);
1178	}
1179	if(c == EOF) return 0;
1180	else return 1;
1181	}
1182
1183	/*********************************************************/
1184
1185	void Uppercase(char *ch)
1186	{
1187	/* convert ch to upper case -- either ASCII or EBCDIC */
1188	ch = isupper((int)ch) ? ch : toupper((int)ch);
1189	}
1190
1191	/*********************************************************/
1192
1193	allseq Compact_Seq(seq data, option input)
1194	{
1195	allseq *alldata;
1196	int i,j,k/,diff/,site;
1197	int n_patt,which_patt,n_invar;
1198	char **sp_names;
1199	int n_otu;
1200
1201
1202	n_otu = input->mod->n_otu;
1203
1204	sp_names = (char *)mCalloc(n_otu,sizeof(char ));
1205	For(i,n_otu)
1206	{
1207	sp_names[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1208	strcpy(sp_names[i],data[i]->name);
1209	}
1210
1211	alldata = Make_Seq(n_otu,data[0]->len,sp_names);
1212
1213	For(i,n_otu) Free(sp_names[i]);
1214	Free(sp_names);
1215
1216	n_patt = which_patt = 0;
1217	/* diff = -1;*/
1218
1219
1220	if(data[0]->len%input->mod->stepsize)
1221	{
1222	printf("\n. Sequence length is not a multiple of %d\n",input->mod->stepsize);
1223	Exit("");
1224	}
1225
1226	Fors(site,data[0]->len,input->mod->stepsize)
1227	{
1228	Fors(k,n_patt,input->mod->stepsize)
1229	{
1230	For(j,n_otu)
1231	{
1232	if(/!Compare_Two_States/strncmp(alldata->c_seq[j]->state+k,
1233	data[j]->state+site,
1234	input->mod->stepsize))
1235	break;
1236	}
1237
1238	if(j == n_otu)
1239	{
1240	which_patt = k;
1241	break;
1242	}
1243	}
1244
1245	/* k = n_patt; */
1246
1247	if(k == n_patt)
1248	{
1249	For(j,n_otu)
1250	Copy_One_State(data[j]->state+site,
1251	alldata->c_seq[j]->state+n_patt,
1252	input->mod->stepsize);
1253
1254
1255	for(j=0;j<n_otu;j++)
1256	{
1257	/* if((Is_Ambigu(alldata->c_seq[j]->state+n_patt,input->mod->datatype,input->mod->stepsize) \|\| */
1258	/* (/\!Compare_Two_States\/strncmp(alldata->c_seq[j]->state+n_patt, */
1259	/* alldata->c_seq[0]->state+n_patt, */
1260	/* input->mod->stepsize)))) */
1261	if(!(Are_Compatible(alldata->c_seq[j]->state+n_patt,
1262	alldata->c_seq[0]->state+n_patt,
1263	input->mod->stepsize,
1264	input->mod->datatype)))
1265	break;
1266	}
1267
1268	if(j==n_otu)
1269	{
1270	For(j,n_otu)
1271	{
1272	alldata->invar[n_patt] = Assign_State(alldata->c_seq[j]->state+n_patt,
1273	input->mod->datatype,
1274	input->mod->stepsize);
1275	break;
1276	}
1277	}
1278	else
1279	alldata->invar[n_patt] = -1;
1280
1281	/* Print_Site(alldata,k,n_otu,"\n",input->mod->stepsize); */
1282
1283	alldata->sitepatt[site] = n_patt;
1284	alldata->wght[n_patt] += 1.;
1285	n_patt+=input->mod->stepsize;
1286	}
1287	else
1288	{
1289	alldata->sitepatt[site] = which_patt;
1290	alldata->wght[which_patt] += 1.;
1291	}
1292	}
1293
1294
1295	alldata->init_len = data[0]->len;
1296	alldata->crunch_len = n_patt;
1297	For(i,n_otu) alldata->c_seq[i]->len = n_patt;
1298
1299	/* fprintf(stderr,"%d patterns found\n",n_patt); */
1300
1301	/* For(site,alldata->crunch_len) printf("%1.0f",alldata->wght[site]); */
1302
1303	n_invar=0;
1304	For(i,alldata->crunch_len) if(alldata->invar[i]>-1) n_invar+=(int)alldata->wght[i];
1305
1306	if(input->mod->datatype == NT)
1307	Get_Base_Freqs(alldata);
1308	else
1309	Get_AA_Freqs(alldata);
1310
1311	/* fprintf(stderr,"Average nucleotides frequencies : \n"); */
1312	/* fprintf(stderr,"%f %f %f %f\n", */
1313	/* alldata->b_frq[0], */
1314	/* alldata->b_frq[1], */
1315	/* alldata->b_frq[2], */
1316	/* alldata->b_frq[3]); */
1317	return alldata;
1318	}
1319
1320	/*********************************************************/
1321
1322	allseq Compact_CSeq(allseq data, model *mod)
1323	{
1324	allseq *alldata;
1325	int i,j,k,site;
1326	int n_patt,which_patt;
1327	int n_otu;
1328
1329	n_otu = data->n_otu;
1330
1331	alldata = (allseq *)mCalloc(1,sizeof(allseq));
1332	alldata->n_otu=n_otu;
1333	alldata->c_seq = (seq *)mCalloc(n_otu,sizeof(seq ));
1334	alldata->wght = (double *)mCalloc(data->crunch_len,sizeof(double));
1335	alldata->b_frq = (double *)mCalloc(mod->ns,sizeof(double));
1336	alldata->ambigu = (int *)mCalloc(data->crunch_len,sizeof(int));
1337	alldata->invar = (int *)mCalloc(data->crunch_len,sizeof(int));
1338
1339	alldata->crunch_len = alldata->init_len = -1;
1340	For(j,n_otu)
1341	{
1342	alldata->c_seq[j] = (seq *)mCalloc(1,sizeof(seq));
1343	alldata->c_seq[j]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1344	strcpy(alldata->c_seq[j]->name,data->c_seq[j]->name);
1345	alldata->c_seq[j]->state = (char *)mCalloc(data->crunch_len,sizeof(char));
1346	alldata->c_seq[j]->state[0] = data->c_seq[j]->state[0];
1347	}
1348
1349	n_patt = which_patt = 0;
1350
1351
1352	Fors(site,data->crunch_len,mod->stepsize)
1353	{
1354	Fors(k,n_patt,mod->stepsize)
1355	{
1356	For(j,n_otu)
1357	{
1358	if(/!Compare_Two_States/strncmp(alldata->c_seq[j]->state+k,
1359	data->c_seq[j]->state+site,
1360	mod->stepsize))
1361	break;
1362	}
1363
1364	if(j == n_otu)
1365	{
1366	which_patt = k;
1367	break;
1368	}
1369	}
1370
1371	if(k == n_patt)
1372	{
1373	For(j,n_otu) Copy_One_State(data->c_seq[j]->state+site,
1374	alldata->c_seq[j]->state+n_patt,
1375	mod->stepsize);
1376
1377	for(j=1;j<n_otu;j++)
1378	if(/!Compare_Two_States/strncmp(alldata->c_seq[j]->state+n_patt,
1379	alldata->c_seq[j-1]->state+n_patt,
1380	mod->stepsize)) break;
1381
1382	if(j==n_otu) alldata->invar[n_patt] = 1;
1383	alldata->wght[n_patt] += data->wght[site];
1384	n_patt+=mod->stepsize;
1385	}
1386	else alldata->wght[which_patt] += data->wght[site];
1387
1388	/* Print_Site(alldata,k,n_otu,"\n",mod->stepsize); */
1389
1390	}
1391
1392	alldata->init_len = data->crunch_len;
1393	alldata->crunch_len = n_patt;
1394	For(i,n_otu) alldata->c_seq[i]->len = n_patt;
1395
1396
1397	(mod->datatype == NT)?
1398	(Get_Base_Freqs(alldata)):
1399	(Get_AA_Freqs(alldata));
1400
1401	return alldata;
1402	}
1403
1404
1405	/*********************************************************/
1406
1407	void Get_Base_Freqs(allseq *data)
1408	{
1409	int i,j,k;
1410	double A,C,G,T;
1411	double fA,fC,fG,fT;
1412	double w;
1413
1414	fA = fC = fG = fT = .25;
1415
1416	For(k,8)
1417	{
1418	A = C = G = T = .0;
1419	For(i,data->n_otu)
1420	{
1421	For(j,data->crunch_len)
1422	{
1423	w = data->wght[j];
1424	if(w)
1425	{
1426	switch(data->c_seq[i]->state[j]){
1427	case 'A' : A+=w;
1428	break;
1429	case 'C' : C+=w;
1430	break;
1431	case 'G' : G+=w;
1432	break;
1433	case 'T' : T+=w;
1434	break;
1435	case 'U' : T+=w;
1436	break;
1437	case 'M' : C+=wfC/(fC+fA); A+=wfA/(fA+fC);
1438	break;
1439	case 'R' : G+=wfG/(fA+fG); A+=wfA/(fA+fG);
1440	break;
1441	case 'W' : T+=wfT/(fA+fT); A+=wfA/(fA+fT);
1442	break;
1443	case 'S' : C+=wfC/(fC+fG); G+=wfG/(fC+fG);
1444	break;
1445	case 'Y' : C+=wfC/(fC+fT); T+=wfT/(fT+fC);
1446	break;
1447	case 'K' : G+=wfG/(fG+fT); T+=wfT/(fT+fG);
1448	break;
1449	case 'B' : C+=wfC/(fC+fG+fT); G+=wfG/(fC+fG+fT); T+=w*fT/(fC+fG+fT);
1450	break;
1451	case 'D' : A+=wfA/(fA+fG+fT); G+=wfG/(fA+fG+fT); T+=w*fT/(fA+fG+fT);
1452	break;
1453	case 'H' : A+=wfA/(fA+fC+fT); C+=wfC/(fA+fC+fT); T+=w*fT/(fA+fC+fT);
1454	break;
1455	case 'V' : A+=wfA/(fA+fC+fG); C+=wfC/(fA+fC+fG); G+=w*fG/(fA+fC+fG);
1456	break;
1457	case 'N' : case 'X' : case '?' : case 'O' : case '-' :
1458	A+=wfA; C+=wfC; G+=wfG; T+=wfT; break;
1459	default : break;
1460	}
1461	}
1462	}
1463	}
1464	fA = A/(A+C+G+T);
1465	fC = C/(A+C+G+T);
1466	fG = G/(A+C+G+T);
1467	fT = T/(A+C+G+T);
1468	}
1469
1470	data->b_frq[0] = fA;
1471	data->b_frq[1] = fC;
1472	data->b_frq[2] = fG;
1473	data->b_frq[3] = fT;
1474	}
1475
1476	/*********************************************************/
1477
1478	void Get_AA_Freqs(allseq *data)
1479	{
1480	int i,j,k;
1481	double A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y;
1482	double fA,fC,fD,fE,fF,fG,fH,fI,fK,fL,fM,fN,fP,fQ,fR,fS,fT,fV,fW,fY;
1483	double w;
1484	double sum;
1485
1486	fA = fC = fD = fE = fF = fG = fH = fI = fK = fL =
1487	fM = fN = fP = fQ = fR = fS = fT = fV = fW = fY = 1./20.;
1488
1489	For(k,8)
1490	{
1491
1492	A = C = D = E = F = G = H = I = K = L =
1493	M = N = P = Q = R = S = T = V = W = Y = .0;
1494
1495	For(i,data->n_otu)
1496	{
1497	For(j,data->crunch_len)
1498	{
1499	w = data->wght[j];
1500	if(w)
1501	{
1502	switch(data->c_seq[i]->state[j]){
1503	case 'A' : A+=w; break;
1504	case 'C' : C+=w; break;
1505	case 'D' : D+=w; break;
1506	case 'E' : E+=w; break;
1507	case 'F' : F+=w; break;
1508	case 'G' : G+=w; break;
1509	case 'H' : H+=w; break;
1510	case 'I' : I+=w; break;
1511	case 'K' : K+=w; break;
1512	case 'L' : L+=w; break;
1513	case 'M' : M+=w; break;
1514	case 'N' : N+=w; break;
1515	case 'P' : P+=w; break;
1516	case 'Q' : Q+=w; break;
1517	case 'R' : R+=w; break;
1518	case 'S' : S+=w; break;
1519	case 'T' : T+=w; break;
1520	case 'V' : V+=w; break;
1521	case 'W' : W+=w; break;
1522	case 'Y' : Y+=w; break;
1523	case 'Z' : Q+=w; break;
1524	case 'X' : case '?' : case 'O' : case '-' :
1525	A+=w*fA;
1526	C+=w*fC;
1527	D+=w*fD;
1528	E+=w*fE;
1529	F+=w*fF;
1530	G+=w*fG;
1531	H+=w*fH;
1532	I+=w*fI;
1533	K+=w*fK;
1534	L+=w*fL;
1535	M+=w*fM;
1536	N+=w*fN;
1537	P+=w*fP;
1538	Q+=w*fQ;
1539	R+=w*fR;
1540	S+=w*fS;
1541	T+=w*fT;
1542	V+=w*fV;
1543	W+=w*fW;
1544	Y+=w*fY;
1545	break;
1546	default : break;
1547	}
1548	}
1549	}
1550	}
1551	sum = (A+C+D+E+F+G+H+I+K+L+M+N+P+Q+R+S+T+V+W+Y);
1552	fA = A/sum; fC = C/sum; fD = D/sum; fE = E/sum;
1553	fF = F/sum; fG = G/sum; fH = H/sum; fI = I/sum;
1554	fK = K/sum; fL = L/sum; fM = M/sum; fN = N/sum;
1555	fP = P/sum; fQ = Q/sum; fR = R/sum; fS = S/sum;
1556	fT = T/sum; fV = V/sum; fW = W/sum; fY = Y/sum;
1557	}
1558
1559	data->b_frq[0] = fA; data->b_frq[1] = fR; data->b_frq[2] = fN; data->b_frq[3] = fD;
1560	data->b_frq[4] = fC; data->b_frq[5] = fQ; data->b_frq[6] = fE; data->b_frq[7] = fG;
1561	data->b_frq[8] = fH; data->b_frq[9] = fI; data->b_frq[10] = fL; data->b_frq[11] = fK;
1562	data->b_frq[12] = fM; data->b_frq[13] = fF; data->b_frq[14] = fP; data->b_frq[15] = fS;
1563	data->b_frq[16] = fT; data->b_frq[17] = fW; data->b_frq[18] = fY; data->b_frq[19] = fV;
1564	}
1565
1566	/*********************************************************/
1567
1568	arbre Read_Tree_File(FILE fp_input_tree)
1569	{
1570	char *line;
1571	arbre *tree;
1572	int i;
1573	char c;
1574
1575	line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1576
1577	do
1578	c=fgetc(fp_input_tree);
1579	while((c != '(') && (c != EOF));
1580
1581	if(c==EOF)
1582	{
1583	Free(line);
1584	return NULL;
1585	}
1586
1587	i=0;
1588	for(;;)
1589	{
1590	if((c == ' ') \|\| (c == '\n'))
1591	{
1592	c=fgetc(fp_input_tree);
1593	if(c==EOF) break;
1594	else continue;
1595	}
1596
1597	line[i]=c;
1598	i++;
1599	c=fgetc(fp_input_tree);
1600	if(c==EOF \|\| c==';') break;
1601	}
1602
1603	tree = Read_Tree(line);
1604	Free(line);
1605	return tree;
1606	}
1607
1608	/*********************************************************/
1609
1610	void Init_Tree_Edges(node a, node d, arbre tree, int cur)
1611	{
1612	int i,dir_a_d;
1613
1614	dir_a_d = -1;
1615	For(i,3) if(a->v[i] == d) {dir_a_d = i; break;}
1616
1617
1618	tree->t_edges[*cur] = a->b[dir_a_d];
1619	tree->t_edges[cur]->num = cur;
1620	cur = cur + 1;
1621
1622	if(d->tax) return;
1623	else
1624	{
1625	For(i,3)
1626	{
1627	if(d->v[i] != a)
1628	Init_Tree_Edges(d,d->v[i],tree,cur);
1629	}
1630	}
1631	}
1632
1633	/*********************************************************/
1634
1635	void Exit(char *message)
1636	{
1637	fprintf(stderr,"%s",message);
1638	exit(1);
1639	}
1640
1641	/*********************************************************/
1642
1643	void *mCalloc(int nb, size_t size)
1644	{
1645	void *allocated;
1646
1647	if((allocated = calloc((size_t)nb,(size_t)size)) != NULL)
1648	{
1649	return allocated;
1650	}
1651	else
1652	Exit("\n. Err: low memory\n");
1653
1654	return NULL;
1655	}
1656
1657	/*********************************************************/
1658
1659	void mRealloc(void p,int nb, size_t size)
1660	{
1661	if((p = realloc(p,(size_t)nb*size)) != NULL)
1662	return p;
1663	else
1664	Exit("\n. Err: low memory\n");
1665
1666	return NULL;
1667	}
1668
1669	/*********************************************************/
1670
1671	arbre *Make_Light_Tree_Struct(int n_otu)
1672	{
1673	arbre *tree;
1674	int i;
1675
1676	tree = (arbre *)mCalloc(1,sizeof(arbre ));
1677	tree->t_edges = (edge *)mCalloc(2n_otu-3,sizeof(edge *));
1678	tree->noeud = (node *)mCalloc(2n_otu-2,sizeof(node *));
1679	tree->n_otu = n_otu;
1680
1681	For(i,2*n_otu-3)
1682	{
1683	tree->t_edges[i] = (edge *)mCalloc(1,sizeof(edge));
1684	Init_Edge_Light(tree->t_edges[i]);
1685	}
1686
1687	For(i,2*n_otu-2)
1688	{
1689	tree->noeud[i] = (node *)mCalloc(1,sizeof(node));
1690	tree->noeud[i]->v = NULL;
1691	Make_Node_Light(tree->noeud[i]);
1692	}
1693	return tree;
1694	}
1695
1696	/*********************************************************/
1697
1698	int Sort_Double_Decrease(const void a, const void b)
1699	{
1700	if(((double )(a)) >= ((double )(b))) return -1;
1701	else return 1;
1702	}
1703
1704	/*********************************************************/
1705
1706	void qksort(double* A, int ilo, int ihi)
1707	{
1708	double pivot; // pivot value for partitioning array
1709	int ulo, uhi; // indices at ends of unpartitioned region
1710	int ieq; // least index of array entry with value equal to pivot
1711	double tempEntry; // temporary entry used for swapping
1712
1713	if (ilo >= ihi) {
1714	return;
1715	}
1716	// Select a pivot value.
1717	pivot = A[(ilo + ihi)/2];
1718	// Initialize ends of unpartitioned region and least index of entry
1719	// with value equal to pivot.
1720	ieq = ulo = ilo;
1721	uhi = ihi;
1722	// While the unpartitioned region is not empty, try to reduce its size.
1723	while (ulo <= uhi) {
1724	if (A[uhi] > pivot) {
1725	// Here, we can reduce the size of the unpartitioned region and
1726	// try again.
1727	uhi--;
1728	} else {
1729	// Here, A[uhi] <= pivot, so swap entries at indices ulo and
1730	// uhi.
1731	tempEntry = A[ulo];
1732	A[ulo] = A[uhi];
1733	A[uhi] = tempEntry;
1734	// After the swap, A[ulo] <= pivot.
1735	if (A[ulo] < pivot) {
1736	// Swap entries at indices ieq and ulo.
1737	tempEntry = A[ieq];
1738	A[ieq] = A[ulo];
1739	A[ulo] = tempEntry;
1740	// After the swap, A[ieq] < pivot, so we need to change
1741	// ieq.
1742	ieq++;
1743	// We also need to change ulo, but we also need to do
1744	// that when A[ulo] = pivot, so we do it after this if
1745	// statement.
1746	}
1747	// Once again, we can reduce the size of the unpartitioned
1748	// region and try again.
1749	ulo++;
1750	}
1751	}
1752	// Now, all entries from index ilo to ieq - 1 are less than the pivot
1753	// and all entries from index uhi to ihi + 1 are greater than the
1754	// pivot. So we have two regions of the array that can be sorted
1755	// recursively to put all of the entries in order.
1756	qksort(A, ilo, ieq - 1);
1757	qksort(A, uhi + 1, ihi);
1758	}
1759
1760	/********************************************************/
1761
1762	void Print_Site(allseq alldata, int num, int n_otu, char sep, int stepsize)
1763	{
1764	int i,j;
1765	For(i,n_otu)
1766	{
1767	printf("%s ",alldata->c_seq[i]->name);
1768	For(j,stepsize)
1769	printf("%c",alldata->c_seq[i]->state[num+j]);
1770	printf("%s",sep);
1771	}
1772	fprintf(stderr,"%s",sep);
1773	}
1774
1775	/*********************************************************/
1776
1777	void Print_Site_Lk(arbre tree, FILE fp)
1778	{
1779	int site;
1780	int catg;
1781	char *s;
1782
1783	s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
1784
1785	fprintf(fp,"Note : P(D\|M) is the probability of site D given the model M (i.e., the site likelihood)\n");
1786	if(tree->mod->n_catg > 1 \|\| tree->mod->invar)
1787	fprintf(fp,"P(D\|M,rr[x]) is the probability of site D given the model M and the relative rate\nof evolution rr[x], where x is the class of rate to be considered.\nWe have P(D\|M) = \\sum_x P(x) x P(D\|M,rr[x]).\n");
1788	fprintf(fp,"\n\n");
1789
1790	sprintf(s,"Site");
1791	fprintf(fp, "%-7s",s);
1792
1793	sprintf(s,"P(D\|M)");
1794	fprintf(fp,"%-16s",s);
1795
1796	if(tree->mod->n_catg > 1)
1797	{
1798	For(catg,tree->mod->n_catg)
1799	{
1800	sprintf(s,"P(D\|M,rr[%d]=%5.4f)",catg+1,tree->mod->rr[catg]);
1801	fprintf(fp,"%-22s",s);
1802	}
1803	}
1804
1805	if(tree->mod->invar)
1806	{
1807	sprintf(s,"P(D\|M,rr[0]=0)");
1808	fprintf(fp,"%-16s",s);
1809	}
1810
1811	fprintf(fp,"\n");
1812
1813	For(site,tree->data->init_len)
1814	{
1815	fprintf(fp,"%-7d",site+1);
1816	fprintf(fp,"%-16g",exp(tree->site_lk[tree->data->sitepatt[site]]));
1817	if(tree->mod->n_catg > 1)
1818	{
1819	For(catg,tree->mod->n_catg)
1820	fprintf(fp,"%-22g",exp(tree->log_site_lk_cat[catg][tree->data->sitepatt[site]]));
1821	}
1822	if(tree->mod->invar)
1823	{
1824	if ((double)tree->data->invar[site] > -0.5)
1825	fprintf(fp,"%-16g",tree->mod->pi[tree->data->invar[tree->data->sitepatt[site]]]);
1826	else
1827	fprintf(fp,"%-16g",0.0);
1828	}
1829	fprintf(fp,"\n");
1830	}
1831
1832	Free(s);
1833	}
1834
1835
1836	/*********************************************************/
1837
1838	void Print_Seq(seq **data, int n_otu)
1839	{
1840	int i,j;
1841
1842	printf("%d\t%d\n",n_otu,data[0]->len);
1843	For(i,n_otu)
1844	{
1845	For(j,23)
1846	{
1847	if(j<(int)strlen(data[i]->name))
1848	putchar(data[i]->name[j]);
1849	else putchar(' ');
1850	}
1851	For(j,data[i]->len) /FLT uncommented/
1852	/* For(j,2000)//FLT commented*/
1853	{
1854	printf("%c",data[i]->state[j]);
1855	}
1856	printf("\n");
1857	}
1858	}
1859
1860	/*********************************************************/
1861
1862	void Print_CSeq(FILE fp, allseq alldata)
1863	{
1864	int i,j,k;
1865	int n_otu;
1866
1867	n_otu = alldata->n_otu;
1868	fprintf(fp,"%d\t%d\n",n_otu,alldata->init_len);
1869	For(i,n_otu)
1870	{
1871	For(j,23)
1872	{
1873	if(j<(int)strlen(alldata->c_seq[i]->name))
1874	fputc(alldata->c_seq[i]->name[j],fp);
1875	else fputc(' ',fp);
1876	}
1877
1878	For(j,alldata->crunch_len)
1879	{
1880	For(k,alldata->wght[j])
1881	fprintf(fp,"%c",alldata->c_seq[i]->state[j]);
1882	}
1883	fprintf(fp,"\n");
1884	}
1885	fprintf(fp,"\n");
1886
1887	/* printf("\t"); */
1888	/* For(j,alldata->crunch_len) */
1889	/* printf("%.0f ",alldata->wght[j]); */
1890	/* printf("\n"); */
1891	}
1892
1893	/*********************************************************/
1894
1895	void Order_Tree_Seq(arbre tree, seq *data)
1896	{
1897	int i,j,n_otu;
1898	seq *buff;
1899
1900	n_otu = tree->n_otu;
1901
1902	For(i,n_otu)
1903	{
1904	For(j,n_otu)
1905	{
1906	if(!strcmp(tree->noeud[i]->name,data[j]->name))
1907	break;
1908	}
1909	buff = data[j];
1910	data[j] = data[i];
1911	data[i] = buff;
1912	}
1913	}
1914
1915	/*********************************************************/
1916
1917	void Order_Tree_CSeq(arbre tree, allseq data)
1918	{
1919	int i,j,n_otu_tree,n_otu_seq;
1920	seq *buff;
1921
1922
1923	n_otu_tree = tree->n_otu;
1924	n_otu_seq = data->n_otu;
1925
1926
1927	if(n_otu_tree != n_otu_seq)
1928	{
1929	/* printf("%d(tree) != %d(seq) \n",n_otu_tree,n_otu_seq); */
1930	Exit("\n. The number of tips in the tree is not the same as the number of sequences\n");
1931	}
1932	For(i,MAX(n_otu_tree,n_otu_seq))
1933	{
1934	For(j,MIN(n_otu_tree,n_otu_seq))
1935	{
1936	if(!strcmp(tree->noeud[i]->name,data->c_seq[j]->name))
1937	break;
1938	}
1939
1940	if(j==MIN(n_otu_tree,n_otu_seq))
1941	{
1942	printf("\n. Err: %s is not found in sequences data set\n",
1943	tree->noeud[i]->name);
1944	Exit("");
1945	}
1946	buff = data->c_seq[j];
1947	data->c_seq[j] = data->c_seq[i];
1948	data->c_seq[i] = buff;
1949	}
1950	}
1951
1952	/*********************************************************/
1953
1954	matrix *Make_Mat(int n_otu)
1955	{
1956	matrix *mat;
1957	int i;
1958
1959	mat = (matrix *)mCalloc(1,sizeof(matrix));
1960
1961	mat->n_otu = n_otu;
1962
1963	mat->P = (double *)mCalloc(n_otu,sizeof(double ));
1964	mat->Q = (double *)mCalloc(n_otu,sizeof(double ));
1965	mat->dist = (double *)mCalloc(n_otu,sizeof(double ));
1966	mat->on_off = (int *)mCalloc(n_otu,sizeof(int));
1967	mat->name = (char *)mCalloc(n_otu,sizeof(char ));
1968	mat->tip_node = (node *)mCalloc(n_otu,sizeof(node ));
1969
1970
1971	For(i,n_otu)
1972	{
1973	mat->P[i] = (double *)mCalloc(n_otu,sizeof(double));
1974	mat->Q[i] = (double *)mCalloc(n_otu,sizeof(double));
1975	mat->dist[i] = (double *)mCalloc(n_otu,sizeof(double));
1976	mat->name[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
1977	}
1978
1979	return mat;
1980	}
1981
1982	/*********************************************************/
1983
1984	void Init_Mat(matrix mat, allseq data)
1985	{
1986	int i;
1987
1988	mat->n_otu = data->n_otu;
1989	mat->r = mat->n_otu;
1990	mat->curr_int = mat->n_otu;
1991	mat->method = 1;
1992
1993	For(i,data->n_otu)
1994	{
1995	strcpy(mat->name[i],data->c_seq[i]->name);
1996	mat->on_off[i] = 1;
1997	}
1998	}
1999
2000	/*********************************************************/
2001
2002	arbre Make_Tree(allseq data)
2003	{
2004	arbre *tree;
2005	int i;
2006
2007	tree = (arbre *)mCalloc(1,sizeof(arbre ));
2008	tree->n_otu = data->n_otu;
2009	Init_Tree(tree);
2010
2011	For(i,2*tree->n_otu-2)
2012	{
2013	tree->noeud[i] = (node *)mCalloc(1,sizeof(node));
2014	tree->noeud[i]->v = NULL;
2015	Make_Node_Light(tree->noeud[i]);
2016	}
2017
2018	For(i,tree->n_otu)
2019	{
2020	strcpy(tree->noeud[i]->name,data->c_seq[i]->name);
2021	tree->noeud[i]->tax = 1;
2022	tree->noeud[i]->num = i;
2023	}
2024
2025	return tree;
2026	}
2027
2028	/*********************************************************/
2029
2030	void Print_Dist(matrix *mat)
2031	{
2032	int i,j;
2033
2034	For(i,mat->n_otu)
2035	{
2036	printf("%s ",mat->name[i]);
2037
2038	For(j,mat->n_otu)
2039	printf("%9.6f ",mat->dist[i][j]);
2040	printf("\n");
2041	}
2042	}
2043
2044	/*********************************************************/
2045
2046	void Print_Node(node a, node d, arbre *tree)
2047	{
2048	int i;
2049
2050	printf("N %2d %2d ",a->num,d->num);
2051	For(i,3) if(a->v[i] == d) {printf("%2d %f\n",
2052	a->b[i]->num,
2053	/* a->b[i]->check_this_one, */
2054	a->b[i]->nj_score); break;}
2055	if(d->tax) return;
2056	else
2057	For(i,3)
2058	if(d->v[i] != a) Print_Node(d,d->v[i],tree);
2059	}
2060
2061	/*********************************************************/
2062
2063	void Share_Lk_Struct(arbre t_full, arbre t_empt)
2064	{
2065	int i,n_otu;
2066	edge b_e,b_f;
2067
2068
2069	n_otu = t_full->n_otu;
2070
2071	t_empt->root = t_empt->noeud[0];
2072
2073	t_empt->tot_loglk_sorted = t_full->tot_loglk_sorted;
2074	t_empt->log_site_lk_cat = t_full->log_site_lk_cat;
2075	t_empt->site_lk = t_full->site_lk;
2076	t_empt->tot_dloglk = t_full->tot_dloglk;
2077	t_empt->tot_d2loglk = t_full->tot_d2loglk;
2078	/* t_empt->mod = t_full->mod; */
2079	/* t_empt->data = t_full->data; */
2080	/* t_empt->mod->s_opt = t_full->mod->s_opt; */
2081
2082	/* For(i,2n_otu-2) /
2083	/* t_empt->noeud[i]->n_ex_nodes = t_full->noeud[i]->n_ex_nodes; */
2084
2085	For(i,2*n_otu-3)
2086	{
2087	b_f = t_full->t_edges[i];
2088	b_e = t_empt->t_edges[i];
2089
2090	b_e->Pij_rr = b_f->Pij_rr;
2091	b_e->dPij_rr = b_f->dPij_rr;
2092	b_e->d2Pij_rr = b_f->d2Pij_rr;
2093
2094	b_e->p_lk_left = b_f->p_lk_left;
2095	b_e->p_lk_rght = b_f->p_lk_rght;
2096
2097	b_e->sum_scale_f_left = b_f->sum_scale_f_left;
2098	b_e->sum_scale_f_rght = b_f->sum_scale_f_rght;
2099
2100	b_e->site_p_lk_left = b_f->site_p_lk_left;
2101	b_e->site_p_lk_rght = b_f->site_p_lk_rght;
2102
2103	b_e->site_dlk_rr = b_f->site_dlk_rr;
2104	b_e->site_d2lk_rr = b_f->site_d2lk_rr;
2105
2106	b_e->ql = b_f->ql;
2107	}
2108	}
2109
2110	/*********************************************************/
2111
2112	/* void Init_Constant() */
2113	/* { */
2114
2115	/* NODE_DEG_MAX = 50; */
2116	/* BRENT_ITMAX = 100; */
2117	/* BRENT_CGOLD = 0.3819660; */
2118	/* BRENT_ZEPS = 1.e-10; */
2119	/* MNBRAK_GOLD = 1.618034; */
2120	/* MNBRAK_GLIMIT = 100.0; */
2121	/* MNBRAK_TINY = 1.e-20; */
2122	/* ALPHA_MIN = 0.04; */
2123	/* ALPHA_MAX = 100; */
2124	/* BL_MIN = 1.e-10; */
2125	/* BL_START = 1.e-03; */
2126	/* BL_MAX = 1.e+05; */
2127	/* MIN_DIFF_LK = 1.e-06; */
2128	/* GOLDEN_R = 0.61803399; */
2129	/* GOLDEN_C = (1.0-GOLDEN_R); */
2130	/* T_MAX_FILE = 200; */
2131	/* T_MAX_LINE = 100000; */
2132	/* T_MAX_NAME = 100; */
2133	/* T_MAX_SEQ = 1000000; */
2134	/* N_MAX_INSERT = 20; */
2135	/* N_MAX_OTU = 4000; */
2136	/* UNLIKELY = -1.e10; */
2137	/* NJ_SEUIL = 0.1; */
2138	/* ROUND_MAX = 100; */
2139	/* DIST_MAX = 2.00; */
2140	/* AROUND_LK = 50.0; */
2141	/* PROP_STEP = 1.0; */
2142	/* T_MAX_ALPHABET = 100; */
2143	/* MDBL_MIN = 2.225074E-308; */
2144	/* MDBL_MAX = 1.797693E+308; */
2145	/* POWELL_ITMAX = 200; */
2146	/* LINMIN_TOL = 2.0E-04; */
2147	/* LIM_SCALE = 3; */
2148	/* LIM_SCALE_VAL = 1.E-50; */
2149	/* /\* LIM_SCALE = 300; \/ /
2150	/* /\* LIM_SCALE_VAL = 1.E-500; \/ /
2151	/* } */
2152
2153	/*********************************************************/
2154
2155	void Print_Mat(matrix *mat)
2156	{
2157	int i,j;
2158
2159	printf("%d",mat->n_otu);
2160	printf("\n");
2161
2162	For(i,mat->n_otu)
2163	{
2164	For(j,13)
2165	{
2166	if(j>=(int)strlen(mat->name[i])) putchar(' ');
2167	else putchar(mat->name[i][j]);
2168	}
2169
2170	For(j,mat->n_otu)
2171	{
2172	if(mat->dist[i][j] == -1)
2173	printf(" - ");
2174	else
2175	printf("%7.8f ",mat->dist[i][j]);
2176	}
2177	printf("\n");
2178	}
2179	}
2180
2181	/*********************************************************/
2182
2183	int Sort_Edges_Diff_Lk(arbre tree, edge *sorted_edges, int n_elem)
2184	{
2185	int i,j;
2186	edge *buff;
2187
2188	For(i,n_elem-1)
2189	{
2190	for(j=i+1;j<n_elem;j++)
2191	{
2192	if(sorted_edges[j]->diff_lk < sorted_edges[i]->diff_lk)
2193	{
2194	buff = sorted_edges[j];
2195	sorted_edges[j] = sorted_edges[i];
2196	sorted_edges[i] = buff;
2197	}
2198	}
2199	}
2200	return 1;
2201	}
2202	/*********************************************************/
2203
2204	void NNI(arbre tree, edge b_fcus, int do_swap)
2205	{
2206	int l_r, r_l, l_v1, l_v2, r_v3, r_v4;
2207	node v1,v2,v3,v4;
2208	double lk1, lk2, lk3;
2209	double lk1_init, lk2_init, lk3_init;
2210	double bl_init;
2211	double l1,l2,l3;
2212	double l_infa, l_infb, l_max;
2213	/* double lk_infa, lk_infb, lk_max; */
2214	double lk_init;
2215
2216	bl_init = b_fcus->l;
2217	lk_init = tree->tot_loglk;
2218
2219	b_fcus->best_conf = 1;
2220	b_fcus->diff_lk = .0;
2221
2222	lk1 = lk2 = lk3 = UNLIKELY;
2223	v1 = v2 = v3 = v4 = NULL;
2224
2225	l_r = r_l = l_v1 = l_v2 = r_v3 = r_v4 = -1;
2226
2227	l_r = b_fcus->l_r;
2228	r_l = b_fcus->r_l;
2229
2230	v1 = b_fcus->left->v[b_fcus->l_v1];
2231	v2 = b_fcus->left->v[b_fcus->l_v2];
2232	v3 = b_fcus->rght->v[b_fcus->r_v1];
2233	v4 = b_fcus->rght->v[b_fcus->r_v2];
2234
2235
2236	l1 = l2 = l3 = -1.;
2237
2238
2239
2240	/***********/
2241	Swap(v2,b_fcus->left,b_fcus->rght,v3,tree);
2242	tree->mod->s_opt->opt_bl = 0;
2243	tree->both_sides = 1;
2244	lk2_init = Update_Lk_At_Given_Edge(b_fcus,tree);
2245
2246
2247	l_infa = 10.*b_fcus->l;
2248	l_max = b_fcus->l;
2249	l_infb = BL_MIN;
2250
2251	lk2 = Br_Len_Brent(l_infa,l_max,l_infb,
2252	1.e-6,
2253	&(b_fcus->l),
2254	b_fcus,tree,1000);
2255
2256	if(lk2 < lk2_init - MIN_DIFF_LK)
2257	{
2258	printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
2259	printf("%f -- %f \n",lk2_init,lk2);
2260	printf("\n. Err. in Optimize_Br_Len_Serie\n");
2261	}
2262
2263	l2 = b_fcus->l;
2264	Swap(v3,b_fcus->left,b_fcus->rght,v2,tree);
2265	/***********/
2266
2267
2268	/***********/
2269	Swap(v2,b_fcus->left,b_fcus->rght,v4,tree);
2270	b_fcus->l = bl_init;
2271	tree->mod->s_opt->opt_bl = 0;
2272	tree->both_sides = 1;
2273	lk3_init = Update_Lk_At_Given_Edge(b_fcus,tree);
2274
2275
2276	l_infa = 10.*b_fcus->l;
2277	l_max = b_fcus->l;
2278	l_infb = BL_MIN;
2279
2280	lk3 = Br_Len_Brent(l_infa,l_max,l_infb,
2281	1.e-6,
2282	&(b_fcus->l),
2283	b_fcus,tree,1000);
2284
2285	if(lk3 < lk3_init - MIN_DIFF_LK)
2286	{
2287	printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
2288	printf("%f -- %f \n",lk3_init,lk3);
2289	printf("\n. Err. in NNI (2)\n");
2290	}
2291
2292
2293	l3 = b_fcus->l;
2294	Swap(v4,b_fcus->left,b_fcus->rght,v2,tree);
2295	/***********/
2296
2297
2298
2299	/***********/
2300	b_fcus->l = bl_init;
2301	tree->mod->s_opt->opt_bl = 0;
2302	tree->both_sides = 1;
2303
2304	lk1_init = Update_Lk_At_Given_Edge(b_fcus,tree);
2305
2306	if((lk1_init < lk_init - MIN_DIFF_LK) \|\|
2307	(lk1_init > lk_init + MIN_DIFF_LK))
2308	{
2309	printf("\n\n. lk_init = %E; lk = %E\n",
2310	lk_init,
2311	lk1_init);
2312	Exit("\n. Err. in NNI (3)\n");
2313	}
2314
2315	l_infa = 10.*b_fcus->l;
2316	l_max = b_fcus->l;
2317	l_infb = BL_MIN;
2318
2319	lk1 = Br_Len_Brent(l_infa,l_max,l_infb,
2320	1.e-6,
2321	&(b_fcus->l),
2322	b_fcus,tree,1000);
2323
2324	if(lk1 < lk_init - MIN_DIFF_LK)
2325	{
2326	printf("\n\n%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l);
2327	printf("%f -- %f \n",lk1_init,lk1);
2328	printf("\n. Err. in NNI (3)\n");
2329	}
2330
2331	l1 = b_fcus->l;
2332	/***********/
2333
2334
2335
2336	b_fcus->ql[0] = l1;
2337	b_fcus->ql[1] = l2;
2338	b_fcus->ql[2] = l3;
2339
2340
2341	b_fcus->diff_lk = lk1 - MAX(lk2,lk3);
2342
2343
2344
2345	if(lk2 > lk3) b_fcus->best_conf = 2;
2346	else b_fcus->best_conf = 3;
2347
2348
2349	if((do_swap) && ((lk2 > lk1+MDBL_MIN) \|\| (lk3 > lk1+MDBL_MIN)))
2350	{
2351	tree->n_swap++;
2352	printf("Swap edge %d -> %f\n",b_fcus->num,MAX(lk2,lk3));
2353	fflush(stdout);
2354
2355	if(lk2 > lk3)
2356	{
2357	tree->best_loglk = lk2;
2358	Swap(v2,b_fcus->left,b_fcus->rght,v3,tree);
2359	b_fcus->l = l2;
2360	tree->both_sides = 1;
2361	Lk(tree,tree->data);
2362	}
2363	else
2364	{
2365	tree->best_loglk = lk3;
2366	Swap(v2,b_fcus->left,b_fcus->rght,v4,tree);
2367	b_fcus->l = l3;
2368	tree->both_sides = 1;
2369	Lk(tree,tree->data);
2370	}
2371	}
2372	else
2373	{
2374	b_fcus->l = bl_init;
2375	Update_PMat_At_Given_Edge(b_fcus,tree);
2376	tree->tot_loglk = lk_init;
2377	}
2378	}
2379
2380	/*********************************************************/
2381
2382	void Swap(node a, node b, node c, node d, arbre *tree)
2383	{
2384	int ab, ba, cd, dc;
2385	int i;
2386
2387
2388	/* \ /d \ /a
2389	\ / \ /
2390	\b__...__c/ -> \b__...__c/
2391	/ \ / \
2392	/ \ / \
2393	/a \ /d \ */
2394
2395
2396
2397	ab = ba = cd = dc = -1;
2398
2399	For(i,3) if(a->v[i] == b) { ab = i; break; }
2400	For(i,3) if(b->v[i] == a) { ba = i; break; }
2401	For(i,3) if(c->v[i] == d) { cd = i; break; }
2402	For(i,3) if(d->v[i] == c) { dc = i; break; }
2403
2404	a->v[ab] = c;
2405	d->v[dc] = b;
2406	b->v[ba] = d;
2407	c->v[cd] = a;
2408	b->b[ba] = d->b[dc];
2409	c->b[cd] = a->b[ab];
2410
2411
2412	(a->b[ab]->left == b)?
2413	(a->b[ab]->left = c):
2414	(a->b[ab]->rght = c);
2415
2416	(d->b[dc]->left == c)?
2417	(d->b[dc]->left = b):
2418	(d->b[dc]->rght = b);
2419
2420	For(i,3)
2421	{
2422	if(a->b[ab]->left->v[i] == a->b[ab]->rght) a->b[ab]->l_r = i;
2423	if(a->b[ab]->rght->v[i] == a->b[ab]->left) a->b[ab]->r_l = i;
2424	if(d->b[dc]->left->v[i] == d->b[dc]->rght) d->b[dc]->l_r = i;
2425	if(d->b[dc]->rght->v[i] == d->b[dc]->left) d->b[dc]->r_l = i;
2426	}
2427
2428
2429	a->b[ab]->l_v1 = a->b[ab]->l_v2 =
2430	a->b[ab]->r_v1 = a->b[ab]->r_v2 =
2431	d->b[dc]->l_v1 = d->b[dc]->l_v2 =
2432	d->b[dc]->r_v1 = d->b[dc]->r_v2 = -1;
2433
2434	For(i,3)
2435	{
2436	if(i != a->b[ab]->l_r)
2437	{
2438	if(a->b[ab]->l_v1 < 0) a->b[ab]->l_v1 = i;
2439	else a->b[ab]->l_v2 = i;
2440	}
2441	if(i != a->b[ab]->r_l)
2442	{
2443	if(a->b[ab]->r_v1 < 0) a->b[ab]->r_v1 = i;
2444	else a->b[ab]->r_v2 = i;
2445	}
2446	if(i != d->b[dc]->l_r)
2447	{
2448	if(d->b[dc]->l_v1 < 0) d->b[dc]->l_v1 = i;
2449	else d->b[dc]->l_v2 = i;
2450	}
2451	if(i != d->b[dc]->r_l)
2452	{
2453	if(d->b[dc]->r_v1 < 0) d->b[dc]->r_v1 = i;
2454	else d->b[dc]->r_v2 = i;
2455	}
2456	}
2457	}
2458
2459	/*********************************************************/
2460
2461	void Update_All_Partial_Lk(edge b_fcus, arbre tree)
2462	{
2463
2464	Update_SubTree_Partial_Lk(b_fcus->left->b[b_fcus->l_v1],
2465	b_fcus->left,
2466	b_fcus->left->v[b_fcus->l_v1],
2467	tree);
2468
2469	Update_SubTree_Partial_Lk(b_fcus->left->b[b_fcus->l_v2],
2470	b_fcus->left,
2471	b_fcus->left->v[b_fcus->l_v2],
2472	tree);
2473
2474	Update_SubTree_Partial_Lk(b_fcus->rght->b[b_fcus->r_v1],
2475	b_fcus->rght,
2476	b_fcus->rght->v[b_fcus->r_v1],
2477	tree);
2478
2479	Update_SubTree_Partial_Lk(b_fcus->rght->b[b_fcus->r_v2],
2480	b_fcus->rght,
2481	b_fcus->rght->v[b_fcus->r_v2],
2482	tree);
2483
2484	tree->tot_loglk = Lk_At_Given_Edge(tree,b_fcus);
2485	}
2486
2487	/*********************************************************/
2488
2489	void Update_SubTree_Partial_Lk(edge b_fcus, node a, node d, arbre tree)
2490	{
2491	int i;
2492
2493	Update_P_Lk(tree,b_fcus,a);
2494	if(d->tax) return;
2495	else For(i,3) if(d->v[i] != a)
2496	Update_SubTree_Partial_Lk(d->b[i],d,d->v[i],tree);
2497	}
2498
2499	/*********************************************************/
2500
2501	double Update_Lk_At_Given_Edge(edge b_fcus, arbre tree)
2502	{
2503
2504	/* if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; */
2505	/* For(i,tree->mod->n_catg) */
2506	/* { */
2507	/* PMat(b_fcus->ltree->mod->rr[i], /
2508	/* tree->mod, */
2509	/* &b_fcus->Pij_rr[i]); */
2510	/* } */
2511
2512
2513	/* Updating partial likelihood after branch swapping */
2514	Update_P_Lk(tree,b_fcus,b_fcus->left);
2515	Update_P_Lk(tree,b_fcus,b_fcus->rght);
2516
2517	tree->tot_loglk = Lk_At_Given_Edge(tree,b_fcus);
2518	return tree->tot_loglk;
2519	}
2520
2521	/*********************************************************/
2522
2523	void Update_PMat_At_Given_Edge(edge b_fcus, arbre tree)
2524	{
2525
2526	int i;
2527	double len;
2528
2529	len = -1.0;
2530
2531	For(i,tree->mod->n_catg)
2532	{
2533	len = b_fcus->l*tree->mod->rr[i];
2534	if(len < BL_MIN) len = BL_MIN;
2535	PMat(len,tree->mod,&b_fcus->Pij_rr[i]);
2536	}
2537
2538	}
2539
2540	/*********************************************************/
2541
2542	allseq Make_Seq(int n_otu, int len, char *sp_names)
2543	{
2544	allseq *alldata;
2545	int j;
2546
2547	alldata = (allseq *)mCalloc(1,sizeof(allseq));
2548	alldata->n_otu = n_otu;
2549	alldata->c_seq = (seq *)mCalloc(n_otu,sizeof(seq ));
2550	alldata->wght = (double *)mCalloc(len,sizeof(double));
2551	alldata->b_frq = (double *)mCalloc(T_MAX_ALPHABET,sizeof(double));
2552	alldata->ambigu = (int *)mCalloc(len,sizeof(int));
2553	alldata->sitepatt = (int *)mCalloc(len,sizeof(int ));
2554	alldata->invar = (int *)mCalloc(len,sizeof(int));
2555
2556	alldata->crunch_len = alldata->init_len = -1;
2557
2558	For(j,n_otu)
2559	{
2560	alldata->c_seq[j] = (seq *)mCalloc(1,sizeof(seq));
2561	alldata->c_seq[j]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char));
2562	strcpy(alldata->c_seq[j]->name,sp_names[j]);
2563	alldata->c_seq[j]->state = (char *)mCalloc(len+1,sizeof(char));
2564	}
2565	return alldata;
2566	}
2567
2568	/*********************************************************/
2569
2570	allseq Copy_CData(allseq ori, model *mod)
2571	{
2572	allseq *new;
2573	int i,j,n_otu;
2574	char **sp_names;
2575
2576	n_otu = ori->n_otu;
2577
2578	sp_names = (char *)mCalloc(n_otu,sizeof(char ));
2579	For(i,n_otu)
2580	{
2581	sp_names[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char));
2582	strcpy(sp_names[i],ori->c_seq[i]->name);
2583	}
2584
2585	new = Make_Seq(n_otu,ori->init_len,sp_names);
2586
2587	For(i,n_otu) Free(sp_names[i]);
2588	Free(sp_names);
2589
2590	For(i,ori->init_len)
2591	new->sitepatt[i] = ori->sitepatt[i];
2592
2593	For(j,ori->crunch_len)
2594	{
2595	For(i,ori->n_otu)
2596	new->c_seq[i]->state[j] = ori->c_seq[i]->state[j];
2597	new->wght[j] = ori->wght[j];
2598	new->ambigu[j] = ori->ambigu[j];
2599	new->invar[j] = ori->invar[j];
2600	}
2601
2602	For(i,ori->n_otu)
2603	{
2604	new->c_seq[i]->len = ori->c_seq[i]->len;
2605	strcpy(new->c_seq[i]->name,ori->c_seq[i]->name);
2606	}
2607
2608	new->init_len = ori->init_len;
2609	new->clean_len = ori->clean_len;
2610	new->crunch_len = ori->crunch_len;
2611	For(i,mod->ns) new->b_frq[i] = ori->b_frq[i];
2612	new->n_otu = ori->n_otu;
2613	return new;
2614	}
2615
2616	/*********************************************************/
2617
2618	optimiz *Alloc_Optimiz()
2619	{
2620	optimiz *s_opt;
2621	s_opt = (optimiz *)mCalloc(1,sizeof(optimiz));
2622	return s_opt;
2623	}
2624
2625	/*********************************************************/
2626
2627	void Init_Optimiz(optimiz *s_opt)
2628	{
2629	s_opt->print = 1;
2630	s_opt->last_opt = 1;
2631	s_opt->opt_alpha = 0;
2632	s_opt->opt_kappa = 0;
2633	s_opt->opt_bl = 0;
2634	s_opt->opt_pinvar = 0;
2635	s_opt->init_lk = UNLIKELY;
2636	s_opt->n_it_max = 1000;
2637	}
2638
2639	/*********************************************************/
2640
2641	int Filexists(char *filename)
2642	{
2643	FILE *fp;
2644	fp =fopen(filename,"r");
2645	if (fp) {
2646	fclose(fp);
2647	return 1;
2648	} else
2649	return 0;
2650	}
2651
2652	/*********************************************************/
2653
2654	FILE Openfile(char filename, int mode)
2655	{
2656	/* mode = 0 -> read */
2657	/* mode = 1 -> write */
2658	/* mode = 2 -> append */
2659
2660	FILE *fp;
2661	char *s;
2662	int open_test=0;
2663
2664	/* s = (char )mCalloc(T_MAX_FILE,sizeof(char)); /
2665
2666	/* strcpy(s,filename); */
2667
2668	s = filename;
2669
2670	fp = NULL;
2671
2672	switch(mode)
2673	{
2674	case 0 :
2675	{
2676	while(!(fp = (FILE *)fopen(s,"r")) && ++open_test<10)
2677	{
2678	printf("\nCan't open file %s, enter a new name : ",s);
2679	Getstring_Stdin(s);
2680	fflush(stdout);
2681	}
2682	break;
2683	}
2684	case 1 :
2685	{
2686	fp = (FILE *)fopen(s,"w");
2687	break;
2688	}
2689	case 2 :
2690	{
2691	fp = (FILE *)fopen(s,"a");
2692	break;
2693	}
2694
2695	default : break;
2696
2697	}
2698
2699	/* Free(s); */
2700
2701	return fp;
2702	}
2703
2704	/*********************************************************/
2705
2706	void Print_Fp_Out(FILE fp_out, time_t t_beg, time_t t_end, arbre tree, option *input, int n_data_set)
2707	{
2708	char *s;
2709	div_t hour,min;
2710
2711	fprintf(fp_out,". Sequence file : [%s]\n\n", input->seqfile);
2712
2713	/* fprintf(fp_out,". Data set [#%d]\n",n_data_set); FLT*/
2714
2715	(tree->mod->datatype == NT)?
2716	(fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",input->modelname)):
2717	(fprintf(fp_out,". Model of amino acids substitution : %s\n\n",input->modelname));
2718
2719	/was after Sequence file ; moved here FLT/
2720	s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2721	fprintf(fp_out,". Initial tree : [%s]\n\n",
2722	(!input->inputtree)?("BIONJ"):
2723	(strcat(strcat(strcat(s,"user tree ("),input->inputtreefile),")")));
2724	Free(s);
2725
2726	fprintf(fp_out,". Number of taxa : %d\n\n",tree->n_otu);/added FLT/
2727
2728	fprintf(fp_out,"\n");
2729
2730	fprintf(fp_out,". Likelihood : loglk = %.5f\n\n",tree->tot_loglk);/was last ; moved here FLT/
2731
2732	fprintf(fp_out,". Discrete gamma model : %s\n",
2733	(tree->mod->n_catg>1)?("Yes"):("No\n"));
2734	if(tree->mod->n_catg > 1)
2735	{
2736	fprintf(fp_out," - Number of categories : %d\n",tree->mod->n_catg);
2737	fprintf(fp_out," - Gamma shape parameter : %.3f\n\n",tree->mod->alpha);
2738	}
2739
2740	if(tree->mod->invar)
2741	fprintf(fp_out,". Proportion of invariant : %.3f\n\n",tree->mod->pinvar);
2742
2743	/was before Discrete gamma model ; moved here FLT/
2744	if(tree->mod->whichmodel <= 5)
2745	{
2746	fprintf(fp_out,". Transition/transversion ratio : %.3f\n\n",tree->mod->kappa);
2747	}
2748	else if(tree->mod->whichmodel == 6)
2749	{
2750	fprintf(fp_out,". Transition/transversion ratio for purines : %.3f\n",
2751	tree->mod->kappa2.tree->mod->lambda/(1.+tree->mod->lambda));
2752	fprintf(fp_out,". Transition/transversion ratio for pyrimidines : %.3f\n\n",
2753	tree->mod->kappa*2./(1.+tree->mod->lambda));
2754	}
2755
2756	if(tree->mod->datatype == NT)
2757	{
2758	fprintf(fp_out,". Nucleotides frequencies :\n\n");
2759	fprintf(fp_out," - f(A)=%8.5f\n",tree->mod->pi[0]);
2760	fprintf(fp_out," - f(C)=%8.5f\n",tree->mod->pi[1]);
2761	fprintf(fp_out," - f(G)=%8.5f\n",tree->mod->pi[2]);
2762	fprintf(fp_out," - f(T)=%8.5f\n\n",tree->mod->pi[3]);
2763	}
2764
2765
2766
2767	/*****************************************/
2768	if((tree->mod->whichmodel == 7) \|\|
2769	(tree->mod->whichmodel == 8))
2770	{
2771	int i,j;
2772
2773	printf("\n");
2774	fprintf(fp_out,". GTR relative rate parameters : \n\n");
2775	fprintf(fp_out,"A <-> C %8.5f\n",*(tree->mod->rr_param[0]));
2776	fprintf(fp_out,"A <-> G %8.5f\n",*(tree->mod->rr_param[1]));
2777	fprintf(fp_out,"A <-> T %8.5f\n",*(tree->mod->rr_param[2]));
2778	fprintf(fp_out,"C <-> G %8.5f\n",*(tree->mod->rr_param[3]));
2779	fprintf(fp_out,"C <-> T %8.5f\n",*(tree->mod->rr_param[4]));
2780	fprintf(fp_out,"G <-> T 1.0 (fixed)\n\n");
2781
2782
2783	fprintf(fp_out,"\n. Instantaneous rate matrix : \n");
2784	fprintf(fp_out,"\n[A---------C---------G---------T------]\n");
2785	For(i,4)
2786	{
2787	For(j,4)
2788	fprintf(fp_out,"%8.5f ",tree->mod->mat_Q[i*4+j]);
2789	fprintf(fp_out,"\n");
2790	}
2791	fprintf(fp_out,"\n");
2792	fprintf(fp_out,"eg., the instantaneous rate of change from 'C' to 'A' is %8.5f x %8.5f = %8.5f\n\n",
2793	tree->mod->pi[0],
2794	*(tree->mod->rr_param[0]),
2795	tree->mod->mat_Q[1*4+0]);
2796	}
2797	/*****************************************/
2798
2799
2800	hour = div(t_end-t_beg,3600);
2801	min = div(t_end-t_beg,60 );
2802
2803	min.quot -= hour.quot*60;
2804
2805	fprintf(fp_out,". Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60);
2806	if(t_end-t_beg > 60)
2807	fprintf(fp_out,". -> %d seconds\n",(int)(t_end-t_beg));
2808
2809	fprintf(fp_out,"\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n");
2810
2811	fflush(fp_out);
2812
2813	}
2814
2815	/*********************************************************/
2816	/FLT wrote this function/
2817	void Print_Fp_Out_Lines(FILE fp_out, time_t t_beg, time_t t_end, arbre tree, option *input, int n_data_set)
2818	{
2819	char *s;
2820	/div_t hour,min;/
2821
2822	if (n_data_set==1)
2823	{
2824
2825	fprintf(fp_out,". Sequence file : [%s]\n\n", input->seqfile);
2826
2827	(tree->mod->datatype == NT)?
2828	(fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",input->modelname)):
2829	(fprintf(fp_out,". Model of amino acids substitution : %s\n\n",input->modelname));
2830
2831	s = (char *)mCalloc(T_MAX_LINE,sizeof(char));
2832	fprintf(fp_out,". Initial tree : [%s]\n\n",
2833	(!input->inputtree)?("BIONJ"):
2834	(strcat(strcat(strcat(s,"user tree ("),input->inputtreefile),")")));
2835	Free(s);
2836
2837	fprintf(fp_out,"\n");
2838
2839	/headline 1/
2840	fprintf(fp_out, ". Data\t");
2841
2842	fprintf(fp_out,"Nb of \t");
2843
2844	fprintf(fp_out,"Likelihood\t");
2845
2846	fprintf(fp_out, "Discrete \t");
2847
2848	if(tree->mod->n_catg > 1)
2849	fprintf(fp_out, "Number of \tGamma shape\t");
2850
2851	fprintf(fp_out,"Proportion of\t");
2852
2853	if(tree->mod->whichmodel <= 6)
2854	fprintf(fp_out,"Transition/ \t");
2855
2856	fprintf(fp_out,"Nucleotides frequencies \t");
2857
2858	if(tree->mod->whichmodel == 7)
2859	fprintf(fp_out,"Instantaneous rate matrix \t");
2860
2861	/* fprintf(fp_out,"Time\t");*/
2862
2863	fprintf(fp_out, "\n");
2864
2865
2866	/headline 2/
2867	fprintf(fp_out, " set\t");
2868
2869	fprintf(fp_out,"taxa\t");
2870
2871	fprintf(fp_out,"loglk \t");
2872
2873	fprintf(fp_out, "gamma model\t");
2874
2875	if(tree->mod->n_catg > 1)
2876	fprintf(fp_out, "categories\tparameter \t");
2877
2878	fprintf(fp_out,"invariant \t");
2879
2880	if(tree->mod->whichmodel <= 6)
2881	fprintf(fp_out,"transversion\t");
2882
2883	fprintf(fp_out,"f(A) f(C) f(G) f(T) \t");
2884
2885	if(tree->mod->whichmodel == 7)
2886	fprintf(fp_out,"[A---------C---------G---------T------]\t");
2887
2888	/* fprintf(fp_out,"used\t");*/
2889
2890	fprintf(fp_out, "\n");
2891
2892
2893	/headline 3/
2894	if(tree->mod->whichmodel == 6) {
2895	fprintf(fp_out," \t \t \t \t");
2896	if(tree->mod->n_catg > 1) fprintf(fp_out," \t \t");
2897	fprintf(fp_out," \t");
2898	fprintf(fp_out,"purines pyrimid.\t");
2899
2900	fprintf(fp_out, "\n");
2901	}
2902
2903	fprintf(fp_out, "\n");
2904	}
2905
2906
2907	/line items/
2908
2909	fprintf(fp_out," #%d\t",n_data_set);
2910
2911	fprintf(fp_out,"%d \t",tree->n_otu);
2912
2913	fprintf(fp_out,"%.5f\t",tree->tot_loglk);
2914
2915	fprintf(fp_out,"%s \t",
2916	(tree->mod->n_catg>1)?("Yes"):("No "));
2917	if(tree->mod->n_catg > 1)
2918	{
2919	fprintf(fp_out,"%d \t",tree->mod->n_catg);
2920	fprintf(fp_out,"%.3f \t",tree->mod->alpha);
2921	}
2922
2923	/if(tree->mod->invar)/
2924	fprintf(fp_out,"%.3f \t",tree->mod->pinvar);
2925
2926	if(tree->mod->whichmodel <= 5)
2927	{
2928	fprintf(fp_out,"%.3f \t",tree->mod->kappa);
2929	}
2930	else if(tree->mod->whichmodel == 6)
2931	{
2932	fprintf(fp_out,"%.3f ",
2933	tree->mod->kappa2.tree->mod->lambda/(1.+tree->mod->lambda));
2934	fprintf(fp_out,"%.3f\t",
2935	tree->mod->kappa*2./(1.+tree->mod->lambda));
2936	}
2937
2938
2939	if(tree->mod->datatype == NT)
2940	{
2941	fprintf(fp_out,"%8.5f ",tree->mod->pi[0]);
2942	fprintf(fp_out,"%8.5f ",tree->mod->pi[1]);
2943	fprintf(fp_out,"%8.5f ",tree->mod->pi[2]);
2944	fprintf(fp_out,"%8.5f\t",tree->mod->pi[3]);
2945	}
2946	/*
2947	hour = div(t_end-t_beg,3600);
2948	min = div(t_end-t_beg,60 );
2949
2950	min.quot -= hour.quot*60;
2951
2952	fprintf(fp_out,"%dh%dm%ds\t", hour.quot,min.quot,(int)(t_end-t_beg)%60);
2953	if(t_end-t_beg > 60)
2954	fprintf(fp_out,". -> %d seconds\t",(int)(t_end-t_beg));
2955	*/
2956
2957	/*****************************************/
2958	if((tree->mod->whichmodel == 7) \|\| (tree->mod->whichmodel == 8))
2959	{
2960	int i,j;
2961
2962	For(i,4)
2963	{
2964	if (i!=0) {
2965	/format/
2966	fprintf(fp_out," \t \t \t \t");
2967	if(tree->mod->n_catg > 1) fprintf(fp_out," \t \t");
2968	fprintf(fp_out," \t \t");
2969	}
2970	For(j,4)
2971	fprintf(fp_out,"%8.5f ",tree->mod->mat_Q[i*4+j]);
2972	if (i<3) fprintf(fp_out,"\n");
2973	}
2974	}
2975	/*****************************************/
2976
2977	fprintf(fp_out, "\n\n");
2978
2979
2980	fflush(fp_out);
2981
2982	}
2983
2984	/*********************************************************/
2985
2986	void Alloc_All_P_Lk(arbre *tree)
2987	{
2988	int i,j,k;
2989	int nbytes;
2990
2991	nbytes = 0;
2992
2993	For(i,2*tree->n_otu-3)
2994	{
2995	tree->t_edges[i]->get_p_lk_left = 1;
2996	tree->t_edges[i]->get_p_lk_rght = 1;
2997
2998	tree->t_edges[i]->p_lk_left =
2999	(double *)mCalloc(tree->data->crunch_len,sizeof(double ));
3000
3001	nbytes += tree->data->crunch_len * sizeof(double **);
3002
3003	tree->t_edges[i]->p_lk_rght =
3004	(double *)mCalloc(tree->data->crunch_len,sizeof(double ));
3005
3006	nbytes += tree->data->clean_len * sizeof(double **);
3007
3008	For(j,tree->data->crunch_len)
3009	{
3010	tree->t_edges[i]->p_lk_left[j] =
3011	(double *)mCalloc(tree->mod->n_catg,sizeof(double ));
3012
3013	nbytes += tree->mod->n_catg * sizeof(double *);
3014
3015	For(k,tree->mod->n_catg)
3016	{
3017	tree->t_edges[i]->p_lk_left[j][k] =
3018	(double *)mCalloc(tree->mod->ns,sizeof(double ));
3019	nbytes += tree->mod->ns * sizeof(double);
3020	}
3021
3022
3023	tree->t_edges[i]->p_lk_rght[j] =
3024	(double *)mCalloc(tree->mod->n_catg,sizeof(double ));
3025
3026	nbytes += tree->mod->n_catg * sizeof(double *);
3027
3028	For(k,tree->mod->n_catg)
3029	{
3030	tree->t_edges[i]->p_lk_rght[j][k] =
3031	(double *)mCalloc(tree->mod->ns,sizeof(double ));
3032	nbytes += tree->mod->ns * sizeof(double);
3033	}
3034	}
3035
3036	tree->t_edges[i]->sum_scale_f_left =
3037	(double *)mCalloc(tree->data->crunch_len,sizeof(double ));
3038
3039	tree->t_edges[i]->sum_scale_f_rght =
3040	(double *)mCalloc(tree->data->crunch_len,sizeof(double ));
3041
3042	}
3043
3044	printf("NBYTES = %d\n",nbytes);
3045	}
3046
3047	/*********************************************************/
3048
3049	matrix K2P_dist(allseq data, double g_shape)
3050	{
3051	int i,j,k;
3052	int diff;
3053	double unc_len;
3054	matrix *mat;
3055	double **len;
3056
3057	len = (double *)mCalloc(data->n_otu,sizeof(double ));
3058	For(i,data->n_otu)
3059	len[i] = (double *)mCalloc(data->n_otu,sizeof(double));
3060
3061	unc_len = .0;
3062
3063	mat = Make_Mat(data->n_otu);
3064	Init_Mat(mat,data);
3065
3066
3067	For(i,data->c_seq[0]->len)
3068	{
3069	For(j,data->n_otu-1)
3070	{
3071	for(k=j+1;k<data->n_otu;k++)
3072	{
3073	if(((data->c_seq[j]->state[i] == 'A' \|\| data->c_seq[j]->state[i] == 'G') &&
3074	(data->c_seq[k]->state[i] == 'C' \|\| data->c_seq[k]->state[i] == 'T'))\|\|
3075	((data->c_seq[j]->state[i] == 'C' \|\| data->c_seq[j]->state[i] == 'T') &&
3076	(data->c_seq[k]->state[i] == 'A' \|\| data->c_seq[k]->state[i] == 'G')))
3077	{
3078	diff++;
3079	mat->Q[j][k]+=data->wght[i];
3080	len[j][k]+=data->wght[i];
3081	len[k][j]=len[j][k];
3082	}
3083
3084	else
3085	if(((data->c_seq[j]->state[i] == 'A' && data->c_seq[k]->state[i] == 'G') \|\|
3086	(data->c_seq[j]->state[i] == 'G' && data->c_seq[k]->state[i] == 'A'))\|\|
3087	((data->c_seq[j]->state[i] == 'C' && data->c_seq[k]->state[i] == 'T') \|\|
3088	(data->c_seq[j]->state[i] == 'T' && data->c_seq[k]->state[i] == 'C')))
3089	{
3090	diff++;
3091	mat->P[j][k]+=data->wght[i];
3092	len[j][k]+=data->wght[i];
3093	len[k][j]=len[j][k];
3094	}
3095	else
3096	if((data->c_seq[j]->state[i] == 'A' \|\|
3097	data->c_seq[j]->state[i] == 'C' \|\|
3098	data->c_seq[j]->state[i] == 'G' \|\|
3099	data->c_seq[j]->state[i] == 'T')&&
3100	(data->c_seq[k]->state[i] == 'A' \|\|
3101	data->c_seq[k]->state[i] == 'C' \|\|
3102	data->c_seq[k]->state[i] == 'G' \|\|
3103	data->c_seq[k]->state[i] == 'T'))
3104	{
3105	len[j][k]+=data->wght[i];
3106	len[k][j]=len[j][k];
3107	}
3108	}
3109	}
3110	}
3111
3112
3113	For(i,data->n_otu-1)
3114	for(j=i+1;j<data->n_otu;j++)
3115	{
3116	if(len[i][j])
3117	{
3118	mat->P[i][j] /= len[i][j];
3119	mat->Q[i][j] /= len[i][j];
3120	}
3121	else
3122	{
3123	mat->P[i][j] = .5;
3124	mat->Q[i][j] = .5;
3125	}
3126
3127	mat->P[j][i] = mat->P[i][j];
3128	mat->Q[j][i] = mat->Q[i][j];
3129
3130
3131	if((1-2mat->P[i][j]-mat->Q[i][j] <= .0) \|\| (1-2mat->Q[i][j] <= .0))
3132	{
3133	mat->dist[i][j] = -1.;
3134	mat->dist[j][i] = -1.;
3135	continue;
3136	}
3137
3138	mat->dist[i][j] = (g_shape/2)*
3139	(pow(1-2*mat->P[i][j]-mat->Q[i][j],-1./g_shape) +
3140	0.5pow(1-2mat->Q[i][j],-1./g_shape) - 1.5);
3141
3142
3143	if(mat->dist[i][j] > DIST_MAX)
3144	{
3145	mat->dist[i][j] = DIST_MAX;
3146	}
3147	mat->dist[j][i] = mat->dist[i][j];
3148	}
3149
3150	For(i,data->n_otu) free(len[i]);
3151	free(len);
3152	return mat;
3153	}
3154
3155	/*********************************************************/
3156
3157	matrix JC69_Dist(allseq data, model *mod)
3158	{
3159	int site,i,j,k;
3160	double unc_len;
3161	matrix *mat;
3162	double **len;
3163
3164
3165	len = (double *)mCalloc(data->n_otu,sizeof(double ));
3166	For(i,data->n_otu)
3167	len[i] = (double *)mCalloc(data->n_otu,sizeof(double));
3168
3169	unc_len = .0;
3170
3171	mat = Make_Mat(data->n_otu);
3172	Init_Mat(mat,data);
3173
3174	Fors(site,data->c_seq[0]->len,mod->stepsize)
3175	{
3176	For(j,data->n_otu-1)
3177	{
3178	for(k=j+1;k<data->n_otu;k++)
3179	{
3180	if((!Is_Ambigu(data->c_seq[j]->state+site,mod->datatype,mod->stepsize)) &&
3181	(!Is_Ambigu(data->c_seq[k]->state+site,mod->datatype,mod->stepsize)))
3182	{
3183	len[j][k]+=data->wght[site];
3184	len[k][j]=len[j][k];
3185	if(/!Compare_Two_States/strncmp(data->c_seq[j]->state+site,
3186	data->c_seq[k]->state+site,
3187	mod->stepsize))
3188	mat->P[j][k]+=data->wght[site];
3189	}
3190	}
3191	}
3192	}
3193
3194
3195	For(i,data->n_otu-1)
3196	for(j=i+1;j<data->n_otu;j++)
3197	{
3198	if(len[i][j])
3199	{
3200	mat->P[i][j] /= len[i][j];
3201	}
3202	else
3203	{
3204	mat->P[i][j] = 1.;
3205	}
3206
3207	mat->P[j][i] = mat->P[i][j];
3208
3209	if((1.-(mod->ns)/(mod->ns-1.)*mat->P[i][j]) < .0)
3210	{
3211	mat->dist[i][j] = DIST_MAX;
3212	}
3213	else
3214	mat->dist[i][j] = -(mod->ns-1.)/(mod->ns)log(1.-(mod->ns)/(mod->ns-1.)mat->P[i][j]);
3215
3216
3217	if(mat->dist[i][j] > DIST_MAX)
3218	{
3219	mat->dist[i][j] = DIST_MAX;
3220	}
3221	mat->dist[j][i] = mat->dist[i][j];
3222	}
3223
3224	For(i,data->n_otu) free(len[i]);
3225	free(len);
3226
3227	return mat;
3228	}
3229
3230	/*********************************************************/
3231
3232	matrix Hamming_Dist(allseq data, model *mod)
3233	{
3234	int i,j,k;
3235	double unc_len;
3236	matrix *mat;
3237	double **len;
3238
3239
3240	len = (double *)mCalloc(data->n_otu,sizeof(double ));
3241	For(i,data->n_otu)
3242	len[i] = (double *)mCalloc(data->n_otu,sizeof(double));
3243
3244	unc_len = .0;
3245
3246	mat = Make_Mat(data->n_otu);
3247	Init_Mat(mat,data);
3248
3249	For(i,data->c_seq[0]->len)
3250	{
3251	For(j,data->n_otu-1)
3252	{
3253	for(k=j+1;k<data->n_otu;k++)
3254	{
3255	if((!Is_Ambigu(data->c_seq[j]->state+i,mod->datatype,mod->stepsize)) &&
3256	(!Is_Ambigu(data->c_seq[k]->state+i,mod->datatype,mod->stepsize)))
3257	{
3258	len[j][k]+=data->wght[i];
3259	len[k][j]=len[j][k];
3260	if(data->c_seq[j]->state[i] != data->c_seq[k]->state[i])
3261	mat->P[j][k]+=data->wght[i];
3262	}
3263	}
3264	}
3265	}
3266
3267
3268	For(i,data->n_otu-1)
3269	for(j=i+1;j<data->n_otu;j++)
3270	{
3271	if(len[i][j])
3272	{
3273	mat->P[i][j] /= len[i][j];
3274	}
3275	else
3276	{
3277	mat->P[i][j] = 1.;
3278	}
3279
3280	mat->P[j][i] = mat->P[i][j];
3281
3282	mat->dist[i][j] = mat->P[i][j];
3283
3284
3285	if(mat->dist[i][j] > DIST_MAX)
3286	{
3287	mat->dist[i][j] = DIST_MAX;
3288	}
3289	mat->dist[j][i] = mat->dist[i][j];
3290	}
3291
3292	For(i,data->n_otu) free(len[i]);
3293	free(len);
3294
3295	return mat;
3296	}
3297
3298	/*********************************************************/
3299
3300	int Is_Ambigu(char *state, int datatype, int stepsize)
3301	{
3302	int i;
3303
3304	if(datatype == NT)
3305	{
3306	For(i,stepsize)
3307	{
3308	if(strchr("MRWSYKBDHVNXO?-.",state[i]))
3309	return 1;
3310	}
3311	}
3312	else
3313	{
3314	if(strchr("X?-.",state[0])) return 1;
3315	}
3316
3317	return 0;
3318	}
3319
3320	/*********************************************************/
3321
3322	void Check_Ambiguities(allseq *data, int datatype, int stepsize)
3323	{
3324	int i,j;
3325
3326	Fors(j,data->crunch_len,stepsize) For(i,data->n_otu)
3327	{
3328	if(Is_Ambigu(data->c_seq[i]->state+j,
3329	datatype,
3330	stepsize))
3331	{
3332	data->ambigu[j] = 1;
3333	break;
3334	}
3335	}
3336	}
3337
3338	/*********************************************************/
3339
3340	int Assign_State(char *c, int datatype, int stepsize)
3341	{
3342	int state[3];
3343	int i;
3344
3345	state[0] = -1;
3346	if(datatype == NT)
3347	{
3348	For(i,stepsize)
3349	{
3350	switch(c[i])
3351	{
3352	case 'A' : state[i]=0; break;
3353	case 'C' : state[i]=1; break;
3354	case 'G' : state[i]=2; break;
3355	case 'T' : state[i]=3; break;
3356	case 'U' : state[i]=3; break;
3357
3358	default : state[i]=-1;
3359	break;
3360	}
3361	}
3362	return (stepsize>1)?(state[0]16+state[1]4+state[2]):(state[0]);
3363	}
3364	else
3365	{
3366	switch(c[0]){
3367	case 'A' : state[0]=0; break;
3368	case 'R' : state[0]=1; break;
3369	case 'N' : state[0]=2; break;
3370	case 'D' : state[0]=3; break;
3371	case 'C' : state[0]=4; break;
3372	case 'Q' : state[0]=5; break;
3373	case 'E' : state[0]=6; break;
3374	case 'G' : state[0]=7; break;
3375	case 'H' : state[0]=8; break;
3376	case 'I' : state[0]=9; break;
3377	case 'L' : state[0]=10; break;
3378	case 'K' : state[0]=11; break;
3379	case 'M' : state[0]=12; break;
3380	case 'F' : state[0]=13; break;
3381	case 'P' : state[0]=14; break;
3382	case 'S' : state[0]=15; break;
3383	case 'T' : state[0]=16; break;
3384	case 'W' : state[0]=17; break;
3385	case 'Y' : state[0]=18; break;
3386	case 'V' : state[0]=19; break;
3387
3388	case 'B' : state[0] = 2; break;
3389	case 'Z' : state[0] = 5; break;
3390
3391	default : state[0]=-1;
3392	break;
3393	}
3394	return state[0];
3395	}
3396	return -1;
3397	}
3398
3399	/*********************************************************/
3400
3401	void Bootstrap(arbre *tree)
3402	{
3403	int *site_num, n_site;
3404	int replicate,j,k;
3405	int position,init_len;
3406	double buff;
3407	allseq *boot_data;
3408	arbre *boot_tree;
3409	model *boot_mod;
3410	matrix *boot_mat;
3411	char *s;
3412	/* double rf; */
3413
3414
3415
3416	tree->mod->s_opt->last_opt = 0;
3417	tree->print_boot_val = 1;
3418
3419	Alloc_Bip(tree);
3420
3421	Get_Bip(tree->noeud[0],
3422	tree->noeud[0]->v[0],
3423	tree);
3424
3425	site_num = (int *)mCalloc(tree->data->init_len,sizeof(int));
3426
3427	n_site = 0;
3428	For(j,tree->data->crunch_len) For(k,tree->data->wght[j])
3429	{
3430	site_num[n_site] = j;
3431	n_site++;
3432	}
3433
3434	boot_data = Copy_CData(tree->data,tree->mod);
3435
3436	boot_tree = NULL;
3437
3438	printf("\n\n. Non parametric bootstrap analysis \n\n");
3439	printf(" ["); fflush(NULL);
3440
3441	For(replicate,tree->mod->bootstrap)
3442	{
3443	For(j,boot_data->crunch_len) boot_data->wght[j] = .0;
3444
3445	init_len = 0;
3446	For(j,boot_data->init_len)
3447	{
3448	buff = rand();
3449	buff /= (RAND_MAX+1.);
3450	buff *= tree->data->init_len;
3451	position = (int)floor(buff);
3452	boot_data->wght[site_num[position]] += 1.;
3453	init_len++;
3454	}
3455
3456	if(init_len != tree->data->init_len) Exit("\n. Pb in copying sequences\n");
3457
3458
3459	(!tree->mod->datatype)?
3460	(Get_Base_Freqs(boot_data)):
3461	(Get_AA_Freqs(boot_data));
3462
3463	boot_mod = Copy_Model(tree->mod);
3464
3465	Init_Model(boot_data,boot_mod);
3466
3467
3468	if(tree->input->inputtree)
3469	{
3470	rewind(tree->input->fp_input_tree);
3471
3472	boot_tree = Read_Tree_File(tree->input->fp_input_tree);
3473	}
3474	else
3475	{
3476	boot_mat = ML_Dist(boot_data,boot_mod);
3477
3478	boot_mat->tree = Make_Tree(boot_data);
3479
3480	Bionj(boot_mat);
3481
3482	boot_tree = boot_mat->tree;
3483
3484	Free_Mat(boot_mat);
3485	}
3486
3487
3488	boot_tree->mod = boot_mod;
3489	boot_tree->input = tree->input;
3490	boot_tree->data = boot_data;
3491	boot_tree->both_sides = 1;
3492	boot_tree->n_pattern = boot_tree->data->crunch_len/
3493	boot_tree->mod->stepsize;
3494
3495	boot_tree->mod->s_opt->print = 0;
3496
3497	Order_Tree_CSeq(boot_tree,boot_data);
3498
3499	Share_Lk_Struct(tree,boot_tree);
3500
3501	Init_P_Lk_Tips(boot_tree);
3502
3503	if(boot_tree->mod->s_opt->opt_topo)
3504	Simu(boot_tree,1000);
3505	else
3506	{
3507	if(boot_tree->mod->s_opt->opt_free_param)
3508	Round_Optimize(boot_tree,boot_tree->data);
3509	else
3510	Lk(boot_tree,boot_data);
3511	}
3512
3513
3514	Alloc_Bip(boot_tree);
3515
3516	Get_Bip(boot_tree->noeud[0],
3517	boot_tree->noeud[0]->v[0],
3518	boot_tree);
3519
3520	Compare_Bip(tree,boot_tree);
3521
3522
3523	if(tree->input->print_boot_trees)
3524	{
3525	s = Write_Tree(boot_tree);
3526	fprintf(tree->input->fp_boot_tree,"%s\n",s);
3527	Free(s);
3528	Print_Fp_Out_Lines(tree->input->fp_boot_stats,0,0,boot_tree,tree->input,replicate+1);
3529	}
3530
3531
3532	/* rf = .0; */
3533	/* For(j,2tree->n_otu-3) /
3534	/* rf += tree->t_edges[j]->bip_score; */
3535
3536
3537	printf(".");
3538	if(!((replicate+1)%10))
3539	{
3540	printf("] %d/%d\n ",replicate+1,tree->mod->bootstrap);
3541	if(replicate != tree->mod->bootstrap-1) printf("[");
3542	}
3543	fflush(NULL);
3544
3545
3546	Free_Tree(boot_tree);
3547
3548	Free_Model(boot_mod);
3549
3550	}
3551	if(((replicate)%10)) printf("] %d/%d\n ",replicate,tree->mod->bootstrap);
3552
3553	if(tree->input->print_boot_trees)
3554	{
3555	fclose(tree->input->fp_boot_tree);
3556	fclose(tree->input->fp_boot_stats);
3557	}
3558
3559	Free_Cseq(boot_data);
3560
3561	Free(site_num);
3562	}
3563
3564	/*********************************************************/
3565
3566	void Update_BrLen_Invar(arbre *tree)
3567	{
3568	int i;
3569	For(i,2tree->n_otu-3) tree->t_edges[i]->l=(1.0-tree->mod->pinvar);
3570	}
3571
3572	/*********************************************************/
3573
3574	void Getstring_Stdin(char *file_name)
3575	{
3576	fgets(file_name,T_MAX_LINE,stdin);
3577	if (strchr(file_name, '\n') != NULL)
3578	*strchr(file_name, '\n') = '\0';
3579	}
3580
3581	/*********************************************************/
3582
3583	void Print_Freq(arbre *tree)
3584	{
3585
3586	switch(tree->mod->datatype)
3587	{
3588	case NT:
3589	{
3590	printf("A : %f\n",tree->mod->pi[0]);
3591	printf("C : %f\n",tree->mod->pi[1]);
3592	printf("G : %f\n",tree->mod->pi[2]);
3593	printf("T : %f\n",tree->mod->pi[3]);
3594
3595	printf("U : %f\n",tree->mod->pi[4]);
3596	printf("M : %f\n",tree->mod->pi[5]);
3597	printf("R : %f\n",tree->mod->pi[6]);
3598	printf("W : %f\n",tree->mod->pi[7]);
3599	printf("S : %f\n",tree->mod->pi[8]);
3600	printf("Y : %f\n",tree->mod->pi[9]);
3601	printf("K : %f\n",tree->mod->pi[10]);
3602	printf("B : %f\n",tree->mod->pi[11]);
3603	printf("D : %f\n",tree->mod->pi[12]);
3604	printf("H : %f\n",tree->mod->pi[13]);
3605	printf("V : %f\n",tree->mod->pi[14]);
3606	printf("N : %f\n",tree->mod->pi[15]);
3607	break;
3608	}
3609	case AA:
3610	{
3611	printf("A : %f\n",tree->mod->pi[0]);
3612	printf("R : %f\n",tree->mod->pi[1]);
3613	printf("N : %f\n",tree->mod->pi[2]);
3614	printf("D : %f\n",tree->mod->pi[3]);
3615	printf("C : %f\n",tree->mod->pi[4]);
3616	printf("Q : %f\n",tree->mod->pi[5]);
3617	printf("E : %f\n",tree->mod->pi[6]);
3618	printf("G : %f\n",tree->mod->pi[7]);
3619	printf("H : %f\n",tree->mod->pi[8]);
3620	printf("I : %f\n",tree->mod->pi[9]);
3621	printf("L : %f\n",tree->mod->pi[10]);
3622	printf("K : %f\n",tree->mod->pi[11]);
3623	printf("M : %f\n",tree->mod->pi[12]);
3624	printf("F : %f\n",tree->mod->pi[13]);
3625	printf("P : %f\n",tree->mod->pi[14]);
3626	printf("S : %f\n",tree->mod->pi[15]);
3627	printf("T : %f\n",tree->mod->pi[16]);
3628	printf("W : %f\n",tree->mod->pi[17]);
3629	printf("Y : %f\n",tree->mod->pi[18]);
3630	printf("V : %f\n",tree->mod->pi[19]);
3631
3632	printf("N : %f\n",tree->mod->pi[20]);
3633	break;
3634	}
3635	default : {break;}
3636	}
3637	}
3638
3639	/*********************************************************/
3640
3641	double Num_Derivatives_One_Param(double (func)(arbre tree), arbre *tree,
3642	double f0, double *param, double stepsize,
3643	double *err, int precise)
3644	{
3645	int i,j;
3646	double errt,fac,hh,**a,ans;
3647	int n_iter;
3648	a = (double *)mCalloc(11,sizeof(double ));
3649	For(i,11) a[i] = (double *)mCalloc(11,sizeof(double));
3650
3651
3652	n_iter = 10; /* */
3653
3654	ans = .0;
3655
3656	if (stepsize == 0.0) Exit("\n. h must be nonzero in Dfridr.");
3657
3658	hh=stepsize;
3659
3660	if(!precise)
3661	{
3662
3663	param = param+hh;
3664	a[0][0] = (*func)(tree);
3665	a[0][0] -= f0;
3666	a[0][0] /= hh;
3667	param = param-hh;
3668
3669	ans = a[0][0];
3670	}
3671	else
3672	{
3673	param = param+hh;
3674	a[0][0] = (*func)(tree);
3675	/* param = param-2hh; /
3676	/* a[0][0] -= (func)(tree); /
3677	/* a[0][0] /= (2.0hh); /
3678	/* param = param+hh; */
3679	a[0][0] -= f0;
3680	a[0][0] /= hh;
3681	param = param-hh;
3682
3683
3684	*err=1e30;
3685	for(i=1;i<n_iter;i++)
3686	{
3687	hh /= 1.4;
3688
3689	/* param = param+hh; */
3690	/* a[0][i] = (func)(tree); /
3691	/* param = param-2hh; /
3692	/* a[0][i] -= (func)(tree); /
3693	/* a[0][i] /= (2.0hh); /
3694	/* param = param+hh; */
3695
3696
3697	param = param+hh;
3698	a[0][i] = (*func)(tree);
3699	/* param = param-2hh; /
3700	/* a[0][i] -= (func)(tree); /
3701	/* a[0][i] /= (2.0hh); /
3702	/* param = param+hh; */
3703	a[0][i] -= f0;
3704	a[0][i] /= hh;
3705	param = param-hh;
3706
3707
3708	fac=1.4*1.4;
3709	for (j=1;j<=i;j++)
3710	{
3711	a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
3712	fac=1.41.4fac;
3713
3714	errt=MAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));
3715
3716	if (errt <= *err)
3717	{
3718	*err=errt;
3719	ans=a[j][i];
3720	}
3721	}
3722
3723	if(fabs(a[i][i]-a[i-1][i-1]) >= 2.0(err))
3724	break;
3725	}
3726	}
3727	For(i,11) Free(a[i]);
3728	Free(a);
3729
3730	return ans;
3731	}
3732
3733	/*********************************************************/
3734
3735	void Num_Derivative_Several_Param(arbre tree, double param, int n_param, double stepsize,
3736	double (func)(arbre tree), double *derivatives)
3737	{
3738	int i;
3739	double err,f0;
3740
3741	f0 = (*func)(tree);
3742
3743	For(i,n_param)
3744	{
3745	derivatives[i] = Num_Derivatives_One_Param(func,
3746	tree,
3747	f0,
3748	param+i,
3749	stepsize,
3750	&err,
3751	0
3752	);
3753	}
3754	}
3755
3756	/*********************************************************/
3757
3758	int Compare_Two_States(char state1, char state2, int state_size)
3759	{
3760	/* 1 the two states are identical */
3761	/* 0 the two states are different */
3762	int i;
3763
3764	For(i,state_size) if(state1[i] != state2[i]) break;
3765
3766	return (i==state_size)?(1):(0);
3767	}
3768
3769	/*********************************************************/
3770
3771	void Copy_One_State(char from, char to, int state_size)
3772	{
3773	int i;
3774	For(i,state_size) to[i] = from[i];
3775	}
3776
3777	/*********************************************************/
3778
3779	model *Make_Model_Basic()
3780	{
3781	model *mod;
3782	int i;
3783
3784	mod = (model *)mCalloc(1,sizeof(model));
3785
3786	mod->custom_mod_string = (char *)mCalloc(T_MAX_LINE,sizeof(char));
3787	mod->user_b_freq = (double *)mCalloc(100,sizeof(double));
3788
3789	mod->rr_param = (double *)mCalloc(6,sizeof(double ));
3790	mod->rr_param_values = (double *)mCalloc(6,sizeof(double));
3791	mod->rr_param_num = (int *)mCalloc(6,sizeof(int ));
3792	mod->n_rr_param_per_cat = (int *)mCalloc(6,sizeof(int));
3793	mod->s_opt = (optimiz *)Alloc_Optimiz();
3794
3795	For(i,6)
3796	mod->rr_param_num[i] = (int *)mCalloc(6,sizeof(int));
3797
3798	return mod;
3799	}
3800
3801	/*********************************************************/
3802
3803	void Make_Model_Complete(model *mod)
3804	{
3805	int i,j;
3806
3807	mod->pi = (double *)mCalloc(mod->ns,sizeof(double));
3808	mod->r_proba = (double *)mCalloc(mod->n_catg,sizeof(double));
3809	mod->rr = (double *)mCalloc(mod->n_catg,sizeof(double));
3810
3811	mod->Pij_rr = (double *)mCalloc(mod->n_catg,sizeof(double ));
3812	mod->dPij_rr = (double *)mCalloc(mod->n_catg,sizeof(double ));
3813	mod->d2Pij_rr = (double *)mCalloc(mod->n_catg,sizeof(double ));
3814
3815	For(i,mod->n_catg)
3816	{
3817	mod->Pij_rr[i] = (double *)mCalloc(mod->ns,sizeof(double ));
3818	mod->dPij_rr[i] = (double *)mCalloc(mod->ns,sizeof(double ));
3819	mod->d2Pij_rr[i] = (double *)mCalloc(mod->ns,sizeof(double ));
3820
3821	For(j,mod->ns)
3822	{
3823	mod->Pij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double));
3824	mod->dPij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double));
3825	mod->d2Pij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double));
3826	}
3827	}
3828
3829	mod->mat_Q = (double )mCalloc(mod->nsmod->ns,sizeof(double));
3830	mod->mat_Vr = (double )mCalloc(mod->nsmod->ns,sizeof(double));
3831	mod->mat_Vi = (double )mCalloc(mod->nsmod->ns,sizeof(double));
3832	mod->vct_eDmr = (double *)mCalloc(mod->ns ,sizeof(double));
3833	mod->vct_ev = (double *)mCalloc(mod->ns ,sizeof(double));
3834	}
3835
3836	/*********************************************************/
3837
3838	model Copy_Model(model ori)
3839	{
3840	int i,j;
3841	model *cpy;
3842
3843
3844	cpy = Make_Model_Basic();
3845
3846	Copy_Optimiz(ori->s_opt,cpy->s_opt);
3847	/* cpy->c_code = ori->c_code; */
3848
3849
3850	cpy->ns = ori->ns;
3851	cpy->n_catg = ori->n_catg;
3852
3853
3854	Make_Model_Complete(cpy);
3855
3856
3857	cpy->datatype = ori->datatype;
3858	cpy->n_otu = ori->n_otu;
3859	cpy->alpha_old = ori->alpha_old;
3860	cpy->kappa_old = ori->alpha_old;
3861	cpy->lambda_old = ori->lambda_old;
3862	cpy->pinvar_old = ori->pinvar_old;
3863	cpy->whichmodel = ori->whichmodel;
3864	cpy->seq_len = ori->seq_len;
3865	cpy->update_eigen = ori->update_eigen;
3866	cpy->kappa = ori->kappa;
3867	cpy->alpha = ori->alpha;
3868	cpy->lambda = ori->lambda;
3869	cpy->bootstrap = ori->bootstrap;
3870	cpy->invar = ori->invar;
3871	cpy->pinvar = ori->pinvar;
3872	cpy->stepsize = ori->stepsize;
3873	cpy->n_diff_rr_param = ori->n_diff_rr_param;
3874
3875
3876
3877	For(i,cpy->n_diff_rr_param)
3878	{
3879	cpy->n_rr_param_per_cat[i] = ori->n_rr_param_per_cat[i];
3880	For(j,cpy->n_rr_param_per_cat[i])
3881	{
3882	cpy->rr_param_num[i][j] = ori->rr_param_num[i][j];
3883	}
3884	}
3885
3886	For(i,6)
3887	{
3888	cpy->rr_param_values[i] = ori->rr_param_values[i];
3889	cpy->rr_param[i] = cpy->rr_param_values+i;
3890	}
3891
3892	For(i,cpy->ns)
3893	{
3894	cpy->pi[i] = ori->pi[i];
3895	cpy->user_b_freq[i] = ori->user_b_freq[i];
3896	}
3897
3898	For(i,cpy->n_catg)
3899	{
3900	cpy->r_proba[i] = ori->r_proba[i];
3901	cpy->rr[i] = ori->rr[i];
3902	}
3903
3904	return cpy;
3905	}
3906
3907	/*********************************************************/
3908
3909	void Set_Defaults_Input(option* input)
3910	{
3911
3912	input->mod->datatype = 0;
3913	strcpy(input->modelname,"HKY");
3914	strcpy(input->nt_or_cd,"nucleotides");
3915	input->n_data_sets = 1;
3916	input->interleaved = 1;
3917	input->inputtree = 0;
3918	input->tree = NULL;
3919	input->phyml_tree_file_open_mode = 1;
3920	input->phyml_stat_file_open_mode = 1;
3921	input->seq_len = -1;
3922	input->n_data_set_asked = -1;
3923	input->print_boot_trees = 1;
3924
3925
3926	}
3927
3928	/*********************************************************/
3929
3930	void Set_Defaults_Model(model *mod)
3931	{
3932	int i;
3933
3934	strcpy(mod->custom_mod_string,"000000");
3935	mod->whichmodel = 4;
3936	mod->n_catg = 1;
3937	mod->kappa = 4.0;
3938	mod->alpha = 2.0;
3939	mod->lambda = 1.0;
3940	mod->bootstrap = 0;
3941	mod->invar = 0;
3942	mod->pinvar = 0.0;
3943	mod->stepsize = 1;
3944	mod->ns = 4;
3945	mod->n_diff_rr_param = 0;
3946
3947	For(i,6)
3948	{
3949	mod->rr_param_values[i] = 1.0;
3950	mod->rr_param[i] = mod->rr_param_values+i;
3951	}
3952
3953	For(i,4) mod->user_b_freq[i] = -1.;
3954
3955	}
3956
3957	/*********************************************************/
3958
3959	void Set_Defaults_Optimiz(optimiz *s_opt)
3960	{
3961	s_opt->opt_alpha = 0;
3962	s_opt->opt_kappa = 0;
3963	s_opt->opt_lambda = 0;
3964	s_opt->opt_bl = 0;
3965	s_opt->opt_pinvar = 0;
3966	s_opt->opt_rr_param = 0;
3967	s_opt->opt_topo = 1;
3968	s_opt->opt_free_param = 1;
3969	}
3970
3971	/*********************************************************/
3972
3973	void Copy_Optimiz(optimiz ori, optimiz cpy)
3974	{
3975	cpy->print = ori->print;
3976	cpy->opt_alpha = ori->opt_alpha;
3977	cpy->opt_kappa = ori->opt_kappa;
3978	cpy->opt_lambda = ori->opt_lambda;
3979	cpy->opt_pinvar = ori->opt_pinvar;
3980	cpy->opt_rr_param = ori->opt_rr_param;
3981	cpy->opt_free_param = ori->opt_free_param;
3982	cpy->opt_bl = ori->opt_bl;
3983	cpy->init_lk = ori->init_lk;
3984	cpy->n_it_max = ori->n_it_max;
3985	cpy->opt_topo = ori->opt_topo;
3986	}
3987
3988	/*********************************************************/
3989
3990	void Get_Bip(node a, node d, arbre *tree)
3991	{
3992	if(d->tax)
3993	{
3994	d->bip_node[0][0] = d;
3995	d->bip_size[0] = 1;
3996	return;
3997	}
3998	else
3999	{
4000	int i,j,k;
4001	int d_a;
4002
4003
4004	d_a = -1;
4005
4006	For(i,3)
4007	{
4008	if(d->v[i] != a)
4009	Get_Bip(d,d->v[i],tree);
4010	else d_a = i;
4011	}
4012
4013	d->bip_size[d_a] = 0;
4014	For(i,3)
4015	if(d->v[i] !=a )
4016	{
4017	For(j,3)
4018	{
4019	if(d->v[i]->v[j] == d)
4020	{
4021	For(k,d->v[i]->bip_size[j])
4022	{
4023	d->bip_node[d_a][d->bip_size[d_a]] = d->v[i]->bip_node[j][k];
4024	strcpy(d->bip_name[d_a][d->bip_size[d_a]],d->v[i]->bip_node[j][k]->name);
4025	d->bip_size[d_a]++;
4026	}
4027	break;
4028	}
4029	}
4030	}
4031
4032	qsort(d->bip_name[d_a],d->bip_size[d_a],sizeof(char *),Sort_String);
4033
4034	For(i,3)
4035	if(a->v[i] == d)
4036	{
4037	a->bip_size[i] = 0;
4038	For(j,tree->n_otu)
4039	{
4040	For(k,d->bip_size[d_a])
4041	{
4042	if(d->bip_node[d_a][k] == tree->noeud[j])
4043	break;
4044	}
4045
4046	if(k == d->bip_size[d_a])
4047	{
4048	a->bip_node[i][a->bip_size[i]] = tree->noeud[j];
4049	strcpy(a->bip_name[i][a->bip_size[i]],tree->noeud[j]->name);
4050	a->bip_size[i]++;
4051	}
4052	}
4053
4054	qsort(a->bip_name[i],a->bip_size[i],sizeof(char *),Sort_String);
4055
4056	if(a->bip_size[i] != tree->n_otu - d->bip_size[d_a])
4057	{
4058	printf("%d %d \n",a->bip_size[i],tree->n_otu - d->bip_size[d_a]);
4059	Exit("\n. Problem in counting bipartitions \n");
4060	}
4061	break;
4062	}
4063	}
4064	}
4065
4066	/*********************************************************/
4067
4068	void Alloc_Bip(arbre *tree)
4069	{
4070	int i,j,k;
4071
4072	tree->has_bip = 1;
4073
4074	For(i,2*tree->n_otu-2)
4075	{
4076	tree->noeud[i]->bip_size = (int *)mCalloc(3,sizeof(int));
4077	tree->noeud[i]->bip_node = (node *)mCalloc(3,sizeof(node ));
4078	tree->noeud[i]->bip_name = (char *)mCalloc(3,sizeof(char ));
4079	For(j,3)
4080	{
4081	tree->noeud[i]->bip_node[j] =
4082	(node *)mCalloc(tree->n_otu,sizeof(node ));
4083
4084	tree->noeud[i]->bip_name[j] =
4085	(char *)mCalloc(tree->n_otu,sizeof(char ));
4086
4087	For(k,tree->n_otu)
4088	tree->noeud[i]->bip_name[j][k] =
4089	(char *)mCalloc(T_MAX_NAME,sizeof(char ));
4090	}
4091	}
4092	}
4093
4094	/*********************************************************/
4095
4096	int Sort_Double_Increase(const void a, const void b)
4097	{
4098	if(((double )(a)) <= ((double )(b))) return -1;
4099	else return 1;
4100	}
4101
4102	/*********************************************************/
4103
4104	int Sort_String(const void a, const void b)
4105	{
4106	return(strcmp(((const char )(a)), ((const char **)(b))));
4107	}
4108
4109	/*********************************************************/
4110
4111	void Compare_Bip(arbre tree1, arbre tree2)
4112	{
4113	int i,j,k;
4114	edge b1,b2;
4115	char bip1,bip2;
4116	int bip_size;
4117
4118
4119	For(i,2*tree1->n_otu-3)
4120	{
4121	if((!tree1->t_edges[i]->left->tax) &&
4122	(!tree1->t_edges[i]->rght->tax))
4123	{
4124
4125	b1 = tree1->t_edges[i];
4126
4127	For(j,2*tree2->n_otu-3)
4128	{
4129	if((!tree2->t_edges[j]->left->tax) &&
4130	(!tree2->t_edges[j]->rght->tax))
4131	{
4132
4133	b2 = tree2->t_edges[j];
4134
4135	if(MIN(b1->left->bip_size[b1->l_r],b1->rght->bip_size[b1->r_l]) ==
4136	MIN(b2->left->bip_size[b2->l_r],b2->rght->bip_size[b2->r_l]))
4137	{
4138	bip_size = MIN(b1->left->bip_size[b1->l_r],b1->rght->bip_size[b1->r_l]);
4139
4140	if(b1->left->bip_size[b1->l_r] == b1->rght->bip_size[b1->r_l])
4141	{
4142	if(b1->left->bip_name[b1->l_r][0][0] < b1->rght->bip_name[b1->r_l][0][0])
4143	{
4144	bip1 = b1->left->bip_name[b1->l_r];
4145	}
4146	else
4147	{
4148	bip1 = b1->rght->bip_name[b1->r_l];
4149	}
4150	}
4151	else if(b1->left->bip_size[b1->l_r] < b1->rght->bip_size[b1->r_l])
4152	{
4153	bip1 = b1->left->bip_name[b1->l_r];
4154	}
4155	else
4156	{
4157	bip1 = b1->rght->bip_name[b1->r_l];
4158	}
4159
4160
4161	if(b2->left->bip_size[b2->l_r] == b2->rght->bip_size[b2->r_l])
4162	{
4163	if(b2->left->bip_name[b2->l_r][0][0] < b2->rght->bip_name[b2->r_l][0][0])
4164	{
4165	bip2 = b2->left->bip_name[b2->l_r];
4166	}
4167	else
4168	{
4169	bip2 = b2->rght->bip_name[b2->r_l];
4170	}
4171	}
4172	else if(b2->left->bip_size[b2->l_r] < b2->rght->bip_size[b2->r_l])
4173	{
4174	bip2 = b2->left->bip_name[b2->l_r];
4175	}
4176	else
4177	{
4178	bip2 = b2->rght->bip_name[b2->r_l];
4179	}
4180
4181	if(bip_size == 1) Exit("\n. Problem in Compare_Bip\n");
4182
4183
4184	For(k,bip_size)
4185	{
4186	if(strcmp(bip1[k],bip2[k])) break;
4187	}
4188
4189	if(k == bip_size)
4190	{
4191	b1->bip_score++;
4192	b2->bip_score++;
4193	break;
4194	}
4195	}
4196	}
4197	}
4198	}
4199	}
4200	}
4201
4202	/*********************************************************/
4203
4204	void Test_Multiple_Data_Set_Format(option *input)
4205	{
4206	char *line;
4207
4208	line = (char *)mCalloc(T_MAX_LINE,sizeof(char));
4209
4210	input->n_trees = 0;
4211
4212	while(fgets(line,T_MAX_LINE,input->fp_input_tree)) if(strstr(line,";")) input->n_trees++;
4213
4214	Free(line);
4215
4216	/* if((input->n_trees != input->n_data_sets) && */
4217	/* (input->n_data_sets > 1)) */
4218	/* Exit("\n. The number of trees should be the same as\n the number of data sets\n\n"); */
4219
4220	if((input->mod->bootstrap > 1) && (input->n_trees > 1))
4221	Exit("\n. Bootstrap option is not allowed with multiple trees\n");
4222
4223
4224	rewind(input->fp_input_tree);
4225
4226	return;
4227	}
4228
4229	/*********************************************************/
4230
4231	int Are_Compatible(char statea, char stateb, int stepsize, int datatype)
4232	{
4233	int i,j;
4234	char a,b;
4235
4236
4237	if(datatype == NT)
4238	{
4239	For(i,stepsize)
4240	{
4241	a = statea[i];
4242	For(j,stepsize)
4243	{
4244	b = stateb[j];
4245
4246	switch(a)
4247	{
4248	case 'A':
4249	{
4250	switch(b)
4251	{
4252	case 'A' :
4253	case 'M' :
4254	case 'R' :
4255	case 'W' :
4256	case 'D' :
4257	case 'H' :
4258	case 'V' :
4259	case 'X' : {b=b; break;}
4260	default : return 0;
4261	}
4262	break;
4263	}
4264	case 'G':
4265	{
4266	switch(b)
4267	{
4268	case 'G' :
4269	case 'R' :
4270	case 'S' :
4271	case 'K' :
4272	case 'B' :
4273	case 'D' :
4274	case 'V' :
4275	case 'X' : {b=b; break;}
4276	default : return 0;
4277	}
4278	break;
4279	}
4280	case 'C':
4281	{
4282	switch(b)
4283	{
4284	case 'C' :
4285	case 'M' :
4286	case 'S' :
4287	case 'Y' :
4288	case 'B' :
4289	case 'H' :
4290	case 'V' :
4291	case 'X' : {b=b; break;}
4292	default : return 0;
4293	}
4294	break;
4295	}
4296	case 'T':
4297	{
4298	switch(b)
4299	{
4300	case 'T' :
4301	case 'W' :
4302	case 'Y' :
4303	case 'K' :
4304	case 'B' :
4305	case 'D' :
4306	case 'H' :
4307	case 'X' :
4308	{b=b; break;}
4309	default : return 0;
4310	}
4311	break;
4312	}
4313	case 'M' :
4314	{
4315	switch(b)
4316	{
4317	case 'M' :
4318	case 'A' :
4319	case 'C' :
4320	case 'R' :
4321	case 'W' :
4322	case 'S' :
4323	case 'Y' :
4324	case 'B' :
4325	case 'D' :
4326	case 'H' :
4327	case 'V' :
4328	case 'X' :
4329	{b=b; break;}
4330	default : return 0;
4331	}
4332	break;
4333	}
4334	case 'R' :
4335	{
4336	switch(b)
4337	{
4338	case 'R' :
4339	case 'A' :
4340	case 'G' :
4341	case 'M' :
4342	case 'W' :
4343	case 'S' :
4344	case 'K' :
4345	case 'B' :
4346	case 'D' :
4347	case 'H' :
4348	case 'V' :
4349	case 'X' : {b=b; break;}
4350	default : return 0;
4351	}
4352	break;
4353	}
4354
4355	case 'W' :
4356	{
4357	switch(b)
4358	{
4359	case 'W' :
4360	case 'A' :
4361	case 'T' :
4362	case 'M' :
4363	case 'R' :
4364	case 'Y' :
4365	case 'K' :
4366	case 'B' :
4367	case 'D' :
4368	case 'H' :
4369	case 'V' :
4370	case 'X' : {b=b; break;}
4371	default : return 0;
4372	}
4373	break;
4374	}
4375
4376	case 'S' :
4377	{
4378	switch(b)
4379	{
4380	case 'S' :
4381	case 'C' :
4382	case 'G' :
4383	case 'M' :
4384	case 'R' :
4385	case 'Y' :
4386	case 'K' :
4387	case 'B' :
4388	case 'D' :
4389	case 'H' :
4390	case 'V' :
4391	case 'X' : {b=b; break;}
4392	default : return 0;
4393	}
4394	break;
4395	}
4396
4397	case 'Y' :
4398	{
4399	switch(b)
4400	{
4401	case 'Y' :
4402	case 'C' :
4403	case 'T' :
4404	case 'M' :
4405	case 'W' :
4406	case 'S' :
4407	case 'K' :
4408	case 'B' :
4409	case 'D' :
4410	case 'H' :
4411	case 'V' :
4412	case 'X' : {b=b; break;}
4413	default : return 0;
4414	}
4415	break;
4416	}
4417
4418	case 'K' :
4419	{
4420	switch(b)
4421	{
4422	case 'K' :
4423	case 'G' :
4424	case 'T' :
4425	case 'R' :
4426	case 'W' :
4427	case 'S' :
4428	case 'Y' :
4429	case 'B' :
4430	case 'D' :
4431	case 'H' :
4432	case 'V' :
4433	case 'X' : {b=b; break;}
4434	default : return 0;
4435	}
4436	break;
4437	}
4438	case 'B' :
4439	{
4440	switch(b)
4441	{
4442	case 'B' :
4443	case 'C' :
4444	case 'G' :
4445	case 'T' :
4446	case 'M' :
4447	case 'R' :
4448	case 'W' :
4449	case 'S' :
4450	case 'Y' :
4451	case 'K' :
4452	case 'D' :
4453	case 'H' :
4454	case 'V' :
4455	case 'X' : {b=b; break;}
4456	default : return 0;
4457	}
4458	break;
4459	}
4460	case 'D' :
4461	{
4462	switch(b)
4463	{
4464	case 'D' :
4465	case 'A' :
4466	case 'G' :
4467	case 'T' :
4468	case 'M' :
4469	case 'R' :
4470	case 'W' :
4471	case 'S' :
4472	case 'Y' :
4473	case 'K' :
4474	case 'B' :
4475	case 'H' :
4476	case 'V' :
4477	case 'X' : {b=b; break;}
4478	default : return 0;
4479	}
4480	break;
4481	}
4482	case 'H' :
4483	{
4484	switch(b)
4485	{
4486	case 'H' :
4487	case 'A' :
4488	case 'C' :
4489	case 'T' :
4490	case 'M' :
4491	case 'R' :
4492	case 'W' :
4493	case 'S' :
4494	case 'Y' :
4495	case 'K' :
4496	case 'B' :
4497	case 'D' :
4498	case 'V' :
4499	case 'X' : {b=b; break;}
4500	default : return 0;
4501	}
4502	break;
4503	}
4504	case 'V' :
4505	{
4506	switch(b)
4507	{
4508	case 'V' :
4509	case 'A' :
4510	case 'C' :
4511	case 'G' :
4512	case 'M' :
4513	case 'R' :
4514	case 'W' :
4515	case 'S' :
4516	case 'Y' :
4517	case 'K' :
4518	case 'B' :
4519	case 'D' :
4520	case 'H' :
4521	case 'X' : {b=b; break;}
4522	default : return 0;
4523	}
4524	break;
4525	}
4526	case 'X' :
4527	{
4528	switch(b)
4529	{
4530	case 'X' :
4531	case 'A' :
4532	case 'C' :
4533	case 'G' :
4534	case 'T' :
4535	case 'M' :
4536	case 'R' :
4537	case 'W' :
4538	case 'S' :
4539	case 'Y' :
4540	case 'K' :
4541	case 'B' :
4542	case 'D' :
4543	case 'H' :
4544	case 'V' : {b=b; break;}
4545	default : return 0;
4546	}
4547	break;
4548	}
4549	default :
4550	{
4551	printf("\n. Err. in Are_Compatible\n");
4552	printf("\n. Please check that characters `%c` and `%c`\n",a,b);
4553	printf(" correspond to existing amino-acids.\n");
4554	Exit("\n");
4555	return 0;
4556	}
4557	}
4558	}
4559	}
4560	}
4561	else
4562	{
4563	a = statea[0]; b = stateb[0];
4564	switch(a)
4565	{
4566	case 'A' :
4567	{
4568	switch(b)
4569	{
4570	case 'A' :
4571	case 'X' : {b=b; break;}
4572	default : return 0;
4573	}
4574	break;
4575	}
4576	case 'R' :
4577	{
4578	switch(b)
4579	{
4580	case 'R' :
4581	case 'X' : {b=b; break;}
4582	default : return 0;
4583	}
4584	break;
4585	}
4586	case 'N' :
4587	{
4588	switch(b)
4589	{
4590	case 'N' :
4591	case 'B' :
4592	case 'X' : {b=b; break;}
4593	default : return 0;
4594	}
4595	break;
4596	}
4597	case 'B' :
4598	{
4599	switch(b)
4600	{
4601	case 'N' :
4602	case 'B' :
4603	case 'X' : {b=b; break;}
4604	default : return 0;
4605	}
4606	break;
4607	}
4608	case 'D' :
4609	{
4610	switch(b)
4611	{
4612	case 'D' :
4613	case 'X' : {b=b; break;}
4614	default : return 0;
4615	}
4616	break;
4617	}
4618	case 'C' :
4619	{
4620	switch(b)
4621	{
4622	case 'C' :
4623	case 'X' : {b=b; break;}
4624	default : return 0;
4625	}
4626	break;
4627	}
4628	case 'Q' :
4629	{
4630	switch(b)
4631	{
4632	case 'Q' :
4633	case 'Z' :
4634	case 'X' : {b=b; break;}
4635	default : return 0;
4636	}
4637	break;
4638	}
4639	case 'Z' :
4640	{
4641	switch(b)
4642	{
4643	case 'Q' :
4644	case 'Z' :
4645	case 'X' : {b=b; break;}
4646	default : return 0;
4647	}
4648	break;
4649	}
4650	case 'E' :
4651	{
4652	switch(b)
4653	{
4654	case 'E' :
4655	case 'X' : {b=b; break;}
4656	default : return 0;
4657	}
4658	break;
4659	}
4660	case 'G' :
4661	{
4662	switch(b)
4663	{
4664	case 'G' :
4665	case 'X' : {b=b; break;}
4666	default : return 0;
4667	}
4668	break;
4669	}
4670	case 'H' :
4671	{
4672	switch(b)
4673	{
4674	case 'H' :
4675	case 'X' : {b=b; break;}
4676	default : return 0;
4677	}
4678	break;
4679	}
4680	case 'I' :
4681	{
4682	switch(b)
4683	{
4684	case 'I' :
4685	case 'X' : {b=b; break;}
4686	default : return 0;
4687	}
4688	break;
4689	}
4690	case 'L' :
4691	{
4692	switch(b)
4693	{
4694	case 'L' :
4695	case 'X' : {b=b; break;}
4696	default : return 0;
4697	}
4698	break;
4699	}
4700	case 'K' :
4701	{
4702	switch(b)
4703	{
4704	case 'K' :
4705	case 'X' : {b=b; break;}
4706	default : return 0;
4707	}
4708	break;
4709	}
4710	case 'M' :
4711	{
4712	switch(b)
4713	{
4714	case 'M' :
4715	case 'X' : {b=b; break;}
4716	default : return 0;
4717	}
4718	break;
4719	}
4720	case 'F' :
4721	{
4722	switch(b)
4723	{
4724	case 'F' :
4725	case 'X' : {b=b; break;}
4726	default : return 0;
4727	}
4728	break;
4729	}
4730	case 'P' :
4731	{
4732	switch(b)
4733	{
4734	case 'P' :
4735	case 'X' : {b=b; break;}
4736	default : return 0;
4737	}
4738	break;
4739	}
4740	case 'S' :
4741	{
4742	switch(b)
4743	{
4744	case 'S' :
4745	case 'X' : {b=b; break;}
4746	default : return 0;
4747	}
4748	break;
4749	}
4750	case 'T' :
4751	{
4752	switch(b)
4753	{
4754	case 'T' :
4755	case 'X' : {b=b; break;}
4756	default : return 0;
4757	}
4758	break;
4759	}
4760	case 'W' :
4761	{
4762	switch(b)
4763	{
4764	case 'W' :
4765	case 'X' : {b=b; break;}
4766	default : return 0;
4767	}
4768	break;
4769	}
4770	case 'Y' :
4771	{
4772	switch(b)
4773	{
4774	case 'Y' :
4775	case 'X' : {b=b; break;}
4776	default : return 0;
4777	}
4778	break;
4779	}
4780	case 'V' :
4781	{
4782	switch(b)
4783	{
4784	case 'V' :
4785	case 'X' : {b=b; break;}
4786	default : return 0;
4787	}
4788	break;
4789	}
4790	case 'X' :
4791	{
4792	switch(b)
4793	{
4794	case 'A':case 'R':case 'N' :case 'B' :case 'D' :
4795	case 'C':case 'Q':case 'Z' :case 'E' :case 'G' :
4796	case 'H':case 'I':case 'L' :case 'K' :case 'M' :
4797	case 'F':case 'P':case 'S' :case 'T' :case 'W' :
4798	case 'Y':case 'V': case 'X' : {b=b; break;}
4799	default : return 0;
4800	}
4801	break;
4802	}
4803	default :
4804	{
4805	printf("\n. Err. in Are_Compatible\n");
4806	printf("\n. Please check that characters `%c` and `%c`\n",a,b);
4807	printf(" correspond to existing amino-acids.\n");
4808	Exit("\n");
4809	return 0;
4810	}
4811	}
4812	}
4813	return 1;
4814	}
4815
4816	/*********************************************************/
4817
4818	void Hide_Ambiguities(allseq *data)
4819	{
4820	int i;
4821
4822	For(i,data->crunch_len)
4823	{
4824	if(data->ambigu[i])
4825	{
4826	data->wght[i] = 0.0;
4827	}
4828	}
4829	}
4830
4831	/*********************************************************/
4832	/*********************************************************/
4833	/*********************************************************/
4834	/*********************************************************/
4835	/*********************************************************/

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

Download in other formats: