source: trunk/GDE/PHYML/utilities.c

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