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

Last change on this file was 4073, checked in by westram, 13 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
3PHYML :  a program that  computes maximum likelihood  phylogenies from
4DNA or AA homologous sequences
5
6Copyright (C) Stephane Guindon. Oct 2003 onward
7
8All parts of  the source except where indicated  are distributed under
9the GNU public licence.  See http://www.opensource.org for details.
10
11*/
12
13#include "utilities.h"
14#include "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
26double bico(int n, int k)
27{
28   return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
29}
30
31double 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
41double 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
63double 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
70void Plim_Binom(double pH0, int N, double *pinf, double *psup)
71{
72  *pinf = pH0 - 1.64*sqrt(pH0*(1-pH0)/(double)N);
73  if(*pinf < 0) *pinf = .0;
74  *psup = pH0 + 1.64*sqrt(pH0*(1-pH0)/(double)N);
75}
76
77/*********************************************************/
78
79double 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          + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z
96               +.083333333333333)/x; 
97}
98
99/*********************************************************/
100
101double 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]=b*pn[i+2]-an*pn[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
161double 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((p*xx*exp(g+xx*aa)), 1/xx);
180   if (ch-e<0) return (ch);
181   goto l4;
182l1:
183   if (v>.32) goto l3;
184   ch=0.4;   a=log(1-p);
185l2:
186   q=ch;  p1=1+ch*(4.67+ch);  p2=ch*(6.73+ch*(6.66+ch));
187   t=-0.5+(4.67+2*ch)/p1 - (6.73+ch*(13.32+3*ch))/p2;
188   ch-=(1-exp(a+g+.5*ch+c*aa)*p2/p1)/t;
189   if (fabs(q/ch-1)-.01 <= 0) goto l4;
190   else                       goto l2;
191 
192l3: 
193   x=PointNormal (p);
194   p1=0.222222/v;   ch=v*pow((x*sqrt(p1)+1-p1), 3.0);
195   if (ch>2.2*v+6)  ch=-2*(log(1-p)-c*log(.5*ch)+g);
196l4:
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=p2*exp(xx*aa+g+p1-c*log(ch));   
204   b=t/ch;  a=0.5*t-b*c;
205
206   s1=(210+a*(140+a*(105+a*(84+a*(70+60*a))))) / 420;
207   s2=(420+a*(735+a*(966+a*(1141+1278*a))))/2520;
208   s3=(210+a*(462+a*(707+932*a)))/2520;
209   s4=(252+a*(672+1182*a)+c*(294+a*(889+1740*a)))/5040;
210   s5=(84+264*a+c*(175+606*a))/2520;
211   s6=(120+c*(346+127*c))/5040;
212   ch+=t*(1+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6))))));
213   if (fabs(q/ch-1) > e) goto l4;
214
215   return (ch);
216}
217
218/*********************************************************/
219
220double 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 + ((((y*a4+a3)*y+a2)*y+a1)*y+a0) / ((((y*b4+b3)*y+b2)*y+b1)*y+b0);
244   return (p<0.5 ? -z : z);
245}
246/*********************************************************/
247
248int 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.0*K), t, factor=alfa/beta*K, 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((i*2.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
285arbre *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
355void 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
373void 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
391void 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
491void 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
522char **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
579int 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
596char *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
624void 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
665void Init_Tree(arbre *tree)
666{
667  int i;
668
669  tree->noeud              = (node **)mCalloc(2*tree->n_otu-2,sizeof(node *));
670  tree->t_edges            = (edge **)mCalloc(2*tree->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
693void 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
716void 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
735void 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
776void 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
840void 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
857void 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
875void Make_Node_Lk(node *n)
876{
877/*   n->n_ex_nodes = (int *)mCalloc(2,sizeof(int)); */
878  return;
879}
880
881/*********************************************************/
882
883seq **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   
958seq **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
1025seq **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
1148int 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
1185void 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
1193allseq *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
1322allseq *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
1407void 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+=w*fC/(fC+fA); A+=w*fA/(fA+fC);
1438                    break;
1439                  case 'R' : G+=w*fG/(fA+fG); A+=w*fA/(fA+fG);
1440                    break;
1441                  case 'W' : T+=w*fT/(fA+fT); A+=w*fA/(fA+fT);
1442                    break;
1443                  case 'S' : C+=w*fC/(fC+fG); G+=w*fG/(fC+fG);
1444                    break;
1445                  case 'Y' : C+=w*fC/(fC+fT); T+=w*fT/(fT+fC);
1446                    break;
1447                  case 'K' : G+=w*fG/(fG+fT); T+=w*fT/(fT+fG);
1448                    break;
1449                  case 'B' : C+=w*fC/(fC+fG+fT); G+=w*fG/(fC+fG+fT); T+=w*fT/(fC+fG+fT);
1450                    break;
1451                  case 'D' : A+=w*fA/(fA+fG+fT); G+=w*fG/(fA+fG+fT); T+=w*fT/(fA+fG+fT);
1452                    break;
1453                  case 'H' : A+=w*fA/(fA+fC+fT); C+=w*fC/(fA+fC+fT); T+=w*fT/(fA+fC+fT);
1454                    break;
1455                  case 'V' : A+=w*fA/(fA+fC+fG); C+=w*fC/(fA+fC+fG); G+=w*fG/(fA+fC+fG);
1456                    break;
1457                  case 'N' : case 'X' : case '?' : case 'O' : case '-' : 
1458                    A+=w*fA; C+=w*fC; G+=w*fG; T+=w*fT; 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
1478void 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
1568arbre *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
1610void 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
1635void Exit(char *message)
1636{
1637  fprintf(stderr,"%s",message);
1638  exit(1);
1639}
1640
1641/*********************************************************/
1642
1643void *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
1659void *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
1671arbre *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(2*n_otu-3,sizeof(edge *));
1678  tree->noeud   = (node **)mCalloc(2*n_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
1698int 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
1706void 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
1762void 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
1777void 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
1838void 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
1862void 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
1895void 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
1917void 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
1954matrix *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
1984void 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
2002arbre *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
2030void 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
2046void 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
2063void 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,2*n_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
2155void 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
2183int 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
2204void 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
2382void 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
2461void 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
2489void 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
2501double 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->l*tree->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
2523void 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
2542allseq *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
2570allseq *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
2618optimiz *Alloc_Optimiz()
2619{
2620  optimiz *s_opt;
2621  s_opt = (optimiz *)mCalloc(1,sizeof(optimiz));
2622  return s_opt;
2623}
2624
2625/*********************************************************/
2626
2627void 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       
2641int 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
2654FILE *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
2706void 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->kappa*2.*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*/
2817void 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->kappa*2.*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
2986void 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
3049matrix *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-2*mat->P[i][j]-mat->Q[i][j] <= .0) || (1-2*mat->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.5*pow(1-2*mat->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
3157matrix *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
3232matrix *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
3300int 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
3322void 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
3340int 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
3401void 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,2*tree->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
3566void Update_BrLen_Invar(arbre *tree)
3567{
3568  int i;
3569  For(i,2*tree->n_otu-3) tree->t_edges[i]->l*=(1.0-tree->mod->pinvar);
3570}
3571
3572/*********************************************************/
3573
3574void 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
3583void 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
3641double 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-2*hh; */
3676      /*   a[0][0] -= (*func)(tree); */
3677      /*   a[0][0] /= (2.0*hh); */
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-2*hh; */
3692          /*       a[0][i] -= (*func)(tree); */
3693          /*       a[0][i] /= (2.0*hh); */
3694          /*       *param   = *param+hh; */
3695         
3696         
3697          *param   = *param+hh;
3698          a[0][i]  = (*func)(tree);
3699          /*   *param   = *param-2*hh; */
3700          /*   a[0][i] -= (*func)(tree); */
3701          /*   a[0][i] /= (2.0*hh); */
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.4*1.4*fac;
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
3735void 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
3758int 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
3771void 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
3779model *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
3803void 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->ns*mod->ns,sizeof(double));
3830  mod->mat_Vr   = (double *)mCalloc(mod->ns*mod->ns,sizeof(double)); 
3831  mod->mat_Vi   = (double *)mCalloc(mod->ns*mod->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
3838model *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
3909void 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
3930void 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
3959void 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
3973void 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
3990void 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
4068void 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
4096int 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
4104int Sort_String(const void *a, const void *b)
4105{
4106  return(strcmp((*(const char **)(a)), (*(const char **)(b))));
4107}
4108
4109/*********************************************************/
4110
4111void 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
4204void 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
4231int 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
4818void 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.