source: trunk/GDE/PHYML/models.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: 82.1 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 "models.h"
15#include "eigen.h"
16#include "free.h"
17
18/*********************************************************/
19
20void PMat_K80(double l,double kappa, double ***Pij)
21{
22  double Ts,Tv,e1,e2,aux;
23
24  /*0 => A*/
25  /*1 => C*/
26  /*2 => G*/
27  /*3 => T*/
28
29  /* Ts -> transition*/
30  /* Tv -> transversion*/
31 
32 
33  aux = -2*l/(kappa+2);
34  e1 = exp(aux *2);
35  if (1.0!=kappa)
36    {
37      e2 = exp(aux *(kappa+1));
38      Tv = .25*(1-e1);
39      Ts = .25*(1+e1-2*e2);
40    }
41  else
42    {
43      Ts = Tv = .25*(1-e1);
44    }
45
46
47  (*Pij)[0][1] = (*Pij)[1][0] = Tv; 
48  (*Pij)[0][2] = (*Pij)[2][0] = Ts; 
49  (*Pij)[0][3] = (*Pij)[3][0] = Tv;
50
51  (*Pij)[1][2] = (*Pij)[2][1] = Tv; 
52  (*Pij)[1][3] = (*Pij)[3][1] = Ts; 
53 
54  (*Pij)[2][3] = (*Pij)[3][2] = Tv;
55
56  (*Pij)[0][0] = (*Pij)[1][1] = 
57  (*Pij)[2][2] = (*Pij)[3][3] = 1.-Ts-2.*Tv;
58}
59
60/*********************************************************/
61
62void dPMat_K80(double l, double ***dPij, double rr, double k)
63{
64
65  double aux,e1,e2;
66  double dTsl;
67  double dTvl;
68
69    /* Ts -> transition*/
70    /* Tv -> transversion*/
71
72 
73  aux = -2.*l*rr/(k+2.);
74  e1 = exp(aux *2);
75  if (1.0!=k)
76    e2 = exp(aux *(k+1));
77  else
78    e2 = e1;
79
80  dTsl = -rr/(k+2) * e1 + rr*(k+1)/(k+2) * e2;
81  dTvl = rr/(k+2) * e1;
82
83  /*First derivatives*/
84
85  /* branch lengths */
86
87  (*dPij)[0][0] = (*dPij)[1][1] = 
88  (*dPij)[2][2] = (*dPij)[3][3] = -dTsl-2*dTvl;
89
90  (*dPij)[0][1] = (*dPij)[1][0] = dTvl;
91  (*dPij)[0][2] = (*dPij)[2][0] = dTsl;
92  (*dPij)[0][3] = (*dPij)[3][0] = dTvl;
93
94  (*dPij)[1][2] = (*dPij)[2][1] = dTvl;
95  (*dPij)[1][3] = (*dPij)[3][1] = dTsl;
96
97  (*dPij)[2][3] = (*dPij)[3][2] = dTvl;
98}
99
100/*********************************************************/
101
102void d2PMat_K80(double l, double ***d2Pij, double rr, double k)
103{
104  double e1,e2,aux,aux2;
105  double d2Tsl;
106  double d2Tvl;
107
108    /* Ts -> transition*/
109    /* Tv -> transversion*/
110
111  aux = -2.*l*rr/(k+2.);
112  e1 = exp(aux *2);
113  if (1.0!=k)
114    e2 = exp(aux *(k+1));
115  else
116    e2 = e1;
117
118
119  aux2 = rr*rr/pow(k+2,2);
120  d2Tvl = -4.*aux2 * e1;
121  d2Tsl = -d2Tvl -2*aux2 *pow((k+1),2) * e2;
122
123  /*Scnd derivatives*/
124
125  /* branch lengths */
126  (*d2Pij)[0][0] = (*d2Pij)[1][1] = 
127  (*d2Pij)[2][2] = (*d2Pij)[3][3] = -d2Tsl-2*d2Tvl;
128
129  (*d2Pij)[0][1] = (*d2Pij)[1][0] = d2Tvl;
130  (*d2Pij)[0][2] = (*d2Pij)[2][0] = d2Tsl;
131  (*d2Pij)[0][3] = (*d2Pij)[3][0] = d2Tvl;
132
133  (*d2Pij)[1][2] = (*d2Pij)[2][1] = d2Tvl;
134  (*d2Pij)[1][3] = (*d2Pij)[3][1] = d2Tsl;
135
136  (*d2Pij)[2][3] = (*d2Pij)[3][2] = d2Tvl;
137
138}
139
140/*********************************************************/
141
142void PMat_TN93(double l, model *mod, double ***Pij)
143{
144  int i,j;
145  double e1,e2,e3;
146  double a1t,a2t,bt;
147  double A,C,G,T,R,Y;
148  double kappa1,kappa2;
149  int kappa_has_changed;
150
151
152  A = mod->pi[0]; C = mod->pi[1]; G = mod->pi[2]; T = mod->pi[3];
153  R = A+G;  Y = T+C;
154
155  kappa_has_changed = 0;
156  if(mod->kappa < .0) mod->kappa = 1.0e-5;
157
158  if(mod->whichmodel < 5) { mod->lambda = 1.; }
159  else if(mod->whichmodel == 5)
160    {
161      do
162        {
163          mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/(2.*mod->kappa));
164
165          if(mod->lambda < .0)
166            {
167              mod->kappa += mod->kappa/10.;
168              kappa_has_changed = 1;
169            }
170        }while(mod->lambda < .0);
171    }
172
173
174  if((!mod->s_opt->opt_kappa) && (kappa_has_changed))
175    {
176      printf("\n. WARNING: This transition/transversion ratio\n");
177      printf("  is impossible with these base frequencies!\n");
178      printf("  The ratio is now set to %.3f\n",mod->kappa);
179    }
180
181  kappa2 = mod->kappa*2./(1.+mod->lambda);
182  kappa1 = kappa2 * mod->lambda;
183
184 
185  bt = l/(2.*(A*G*kappa1+C*T*kappa2+R*Y));
186
187  a1t = kappa1;
188  a2t = kappa2;
189  a1t*=bt; a2t*=bt;
190
191  e1 = exp(-a1t*R-bt*Y);
192  e2 = exp(-a2t*Y-bt*R);
193  e3 = exp(-bt);
194
195
196  /*A->A*/(*Pij)[0][0] = A+Y*A/R*e3+G/R*e1; 
197  /*A->C*/(*Pij)[0][1] = C*(1-e3);
198  /*A->G*/(*Pij)[0][2] = G+Y*G/R*e3-G/R*e1;
199  /*A->T*/(*Pij)[0][3] = T*(1-e3);
200
201  /*C->A*/(*Pij)[1][0] = A*(1-e3);
202  /*C->C*/(*Pij)[1][1] = C+R*C/Y*e3+T/Y*e2;
203  /*C->G*/(*Pij)[1][2] = G*(1-e3);
204  /*C->T*/(*Pij)[1][3] = T+R*T/Y*e3-T/Y*e2;
205
206  /*G->A*/(*Pij)[2][0] = A+Y*A/R*e3-A/R*e1;
207  /*G->C*/(*Pij)[2][1] = C*(1-e3);
208  /*G->G*/(*Pij)[2][2] = G+Y*G/R*e3+A/R*e1;
209  /*G->T*/(*Pij)[2][3] = T*(1-e3);
210
211  /*T->A*/(*Pij)[3][0] = A*(1-e3);
212  /*T->C*/(*Pij)[3][1] = C+R*C/Y*e3-C/Y*e2;
213  /*T->G*/(*Pij)[3][2] = G*(1-e3);
214  /*T->T*/(*Pij)[3][3] = T+R*T/Y*e3+C/Y*e2;
215 
216  For(i,4) For(j,4)
217    if((*Pij)[i][j] < MDBL_MIN) (*Pij)[i][j] = MDBL_MIN;
218
219}
220
221/*********************************************************/
222
223void dPMat_TN93(double l, double ***dPij, model *mod, double rr)
224{
225  double kappa1,kappa2;
226  double de1dl,de2dl,de3dl;
227  double A,C,G,T,R,Y;
228  double Z;
229
230
231  A = mod->pi[0]; C = mod->pi[1]; G = mod->pi[2]; T = mod->pi[3];
232  R = A+G; Y = C+T;
233
234  if(mod->whichmodel < 5) { mod->lambda = 1.; }
235  else if(mod->whichmodel == 5)
236    {
237      mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/(2.*mod->kappa));
238    }
239
240  kappa2 = mod->kappa*2./(1.+mod->lambda);
241  kappa1 = kappa2 * mod->lambda;
242
243
244  Z = 2.*A*G*kappa1+2.*C*T*kappa2+2.*R*Y;
245 
246  de1dl = -rr*(kappa1*R/Z + Y/Z)*exp(-kappa1*R/Z*l*rr-Y/Z*l*rr);
247  de2dl = -rr*(kappa2*Y/Z + R/Z)*exp(-kappa2*Y/Z*l*rr-R/Z*l*rr);
248  de3dl = -rr/Z*exp(-l*rr/Z);
249
250  /*A->A*/(*dPij)[0][0] = Y*A/R*de3dl+G/R*de1dl;
251  /*A->C*/(*dPij)[0][1] = -C*de3dl;
252  /*A->G*/(*dPij)[0][2] = Y*G/R*de3dl-G/R*de1dl;
253  /*A->T*/(*dPij)[0][3] = -T*de3dl;
254
255  /*C->A*/(*dPij)[1][0] = -A*de3dl;
256  /*C->C*/(*dPij)[1][1] = R*C/Y*de3dl+T/Y*de2dl;
257  /*C->G*/(*dPij)[1][2] = -G*de3dl;
258  /*C->T*/(*dPij)[1][3] = R*T/Y*de3dl-T/Y*de2dl;
259 
260  /*G->A*/(*dPij)[2][0] = Y*A/R*de3dl-A/R*de1dl;
261  /*G->C*/(*dPij)[2][1] = -C*de3dl;
262  /*G->G*/(*dPij)[2][2] = Y*G/R*de3dl+A/R*de1dl;
263  /*G->T*/(*dPij)[2][3] = -T*de3dl;
264
265  /*T->A*/(*dPij)[3][0] = -A*de3dl;
266  /*T->C*/(*dPij)[3][1] = R*C/Y*de3dl-C/Y*de2dl;
267  /*T->G*/(*dPij)[3][2] = -G*de3dl;
268  /*T->T*/(*dPij)[3][3] = R*T/Y*de3dl+C/Y*de2dl;
269}
270
271/*********************************************************/
272
273void d2PMat_TN93(double l, double ***d2Pij, model *mod, double rr)
274{
275  double kappa1,kappa2;
276  double d2e1dl2,d2e2dl2,d2e3dl2;
277  double A,C,G,T,R,Y;
278  double Z;
279
280
281  A = mod->pi[0]; C = mod->pi[1]; G = mod->pi[2]; T = mod->pi[3];
282  R = A+G; Y = C+T;
283
284
285  if(mod->whichmodel < 5) { mod->lambda = 1.; }
286  else if(mod->whichmodel == 5)
287    {
288      mod->lambda = (Y+(R-Y)/(2.*mod->kappa))/(R-(R-Y)/(2.*mod->kappa));
289    }
290
291  kappa2 = mod->kappa*2./(1.+mod->lambda);
292  kappa1 = kappa2 * mod->lambda;
293
294
295  Z = 2.*A*G*kappa1+2.*C*T*kappa2+2.*R*Y;
296 
297  d2e1dl2 = (-rr*(kappa1*R/Z + Y/Z))*
298            (-rr*(kappa1*R/Z + Y/Z))*
299             exp(-kappa1*R/Z*l*rr-Y/Z*l*rr);
300  d2e2dl2 = (-rr*(kappa2*Y/Z + R/Z))*
301            (-rr*(kappa2*Y/Z + R/Z))*
302             exp(-kappa2*Y/Z*l*rr-R/Z*l*rr);
303  d2e3dl2 = (-rr/Z)*
304            (-rr/Z)*
305             exp(-l*rr/Z);
306
307  /*A->A*/(*d2Pij)[0][0] = Y*A/R*d2e3dl2+G/R*d2e1dl2;
308  /*A->C*/(*d2Pij)[0][1] = -C*d2e3dl2;
309  /*A->G*/(*d2Pij)[0][2] = Y*G/R*d2e3dl2-G/R*d2e1dl2;
310  /*A->T*/(*d2Pij)[0][3] = -T*d2e3dl2;
311
312  /*C->A*/(*d2Pij)[1][0] = -A*d2e3dl2;
313  /*C->C*/(*d2Pij)[1][1] = R*C/Y*d2e3dl2+T/Y*d2e2dl2;
314  /*C->G*/(*d2Pij)[1][2] = -G*d2e3dl2;
315  /*C->T*/(*d2Pij)[1][3] = R*T/Y*d2e3dl2-T/Y*d2e2dl2;
316
317  /*G->A*/(*d2Pij)[2][0] = Y*A/R*d2e3dl2-A/R*d2e1dl2;
318  /*G->C*/(*d2Pij)[2][1] = -C*d2e3dl2;
319  /*G->G*/(*d2Pij)[2][2] = Y*G/R*d2e3dl2+A/R*d2e1dl2;
320  /*G->T*/(*d2Pij)[2][3] = -T*d2e3dl2;
321
322  /*T->A*/(*d2Pij)[3][0] = -A*d2e3dl2;
323  /*T->C*/(*d2Pij)[3][1] = R*C/Y*d2e3dl2-C/Y*d2e2dl2;
324  /*T->G*/(*d2Pij)[3][2] = -G*d2e3dl2;
325  /*T->T*/(*d2Pij)[3][3] = R*T/Y*d2e3dl2+C/Y*d2e2dl2;
326
327}
328
329/*********************************************************/
330
331static int Matinv (double *x, int n, int m)
332{
333/* x[n*m]  ... m>=n
334*/
335   int i,j,k;
336   int *irow;
337   double ee, t,t1,xmax;
338   double det;
339   
340   ee = 1.0E-20;
341   det = 1.0;
342
343   irow = (int *)mCalloc(n,sizeof(int));
344
345
346   For (i,n) 
347     {
348       xmax = 0.;
349       for (j=i; j<n; j++)
350         if (xmax < fabs(x[j*m+i])) 
351           { 
352             xmax = fabs(x[j*m+i]); 
353             irow[i]=j; 
354           }
355
356      det *= xmax;
357      if (xmax < ee)   
358        {
359          printf("\nDet becomes zero at %3d!\t\n", i+1);
360          return(-1);
361        }
362      if (irow[i] != i) 
363        {
364          For (j,m) 
365            {
366              t = x[i*m+j];
367              x[i*m+j] = x[irow[i]*m+j];
368              x[irow[i]*m+j] = t;
369            }
370        }
371      t = 1./x[i*m+i];
372      For (j,n) 
373        {
374          if (j == i) continue;
375          t1 = t*x[j*m+i];
376          For(k,m)  x[j*m+k] -= t1*x[i*m+k];
377          x[j*m+i] = -t1;
378        }
379      For(j,m)   x[i*m+j] *= t;
380      x[i*m+i] = t;
381   }                            /* i  */
382   for (i=n-1; i>=0; i--) 
383     {
384       if (irow[i] == i) continue;
385       For(j,n) 
386         {
387           t = x[j*m+i];
388           x[j*m+i] = x[j*m + irow[i]];
389           x[j*m + irow[i]] = t;
390         }
391     }
392
393   free(irow);
394   return (0);
395}
396
397/********************************************************************/
398/* void PMat_Empirical(double l, model *mod, double ***Pij)         */
399/*                                                                  */
400/* Computes the substitution probability matrix                     */
401/* from the initial substitution rate matrix and frequency vector   */
402/* and one specific branch length                                   */
403/*                                                                  */
404/* input : l , branch length                                        */
405/* input : mod , choosen model parameters, mat_Q and pi             */
406/* ouput : Pij , substitution probability matrix                    */
407/*                                                                  */
408/* matrix P(l) is computed as follows :                             */
409/* P(l) = exp(Q*t) , where :                                        */
410/*                                                                  */
411/*   Q = substitution rate matrix = Vr*D*inverse(Vr) , where :      */
412/*                                                                  */
413/*     Vr = right eigenvector matrix for Q                          */
414/*     D  = diagonal matrix of eigenvalues for Q                    */
415/*                                                                  */
416/*   t = time interval = l / mr , where :                           */
417/*                                                                  */
418/*     mr = mean rate = branch length/time interval                 */
419/*        = sum(i)(pi[i]*p(i->j)) , where :                         */
420/*                                                                  */
421/*       pi = state frequency vector                                */
422/*       p(i->j) = subst. probability from i to a different state   */
423/*               = -Q[ii] , as sum(j)(Q[ij]) +Q[ii] =0              */
424/*                                                                  */
425/* the Taylor development of exp(Q*t) gives :                       */
426/* P(l) = Vr*exp(D*t)        *inverse(Vr)                           */
427/*      = Vr*pow(exp(D/mr),l)*inverse(Vr)                           */
428/*                                                                  */
429/* for performance we compute only once the following matrixes :    */
430/* Vr, inverse(Vr), exp(D/mr)                                       */
431/* thus each time we compute P(l) we only have to :                 */
432/* make 20 times the operation pow()                                */
433/* make 2 20x20 matrix multiplications , that is :                  */
434/*   16000 = 2x20x20x20 times the operation *                       */
435/*   16000 = 2x20x20x20 times the operation +                       */
436/*   which can be reduced to (the central matrix being diagonal) :  */
437/*   8400 = 20x20 + 20x20x20 times the operation *                  */
438/*   8000 = 20x20x20 times the operation +                          */
439/********************************************************************/
440
441void PMat_Empirical(double l, model *mod, double ***Pij)
442{
443  int n = mod->ns;
444  int i, j, k;
445  double *U,*V,*R;
446  double *expt  = (double*)calloc(n,sizeof(double));
447  double *uexpt = (double*)calloc(n*n,sizeof(double));
448
449  U = mod->mat_Vr;
450  V = mod->mat_Vi;
451  R = mod->vct_eDmr;
452
453  For (i,n) For (k,n)
454    (*Pij)[i][k] = .0;
455
456  /* compute pow(exp(D/mr),l) into mat_eDmrl */
457  For (k,n)  expt[k] = pow(R[k], l);
458 
459  /* multiply Vr*pow(exp(D/mr),l)*Vi into Pij */
460  For (i,n) For (k,n)
461    uexpt[i*n+k] = U[i*n+k] * expt[k];
462  For (i,n) 
463    {
464      For (j,n) 
465        {
466          For(k,n)
467            {
468              (*Pij)[i][j] += uexpt[i*n+k] * V[k*n+j];
469            }
470          if((*Pij)[i][j] < MDBL_MIN)
471            (*Pij)[i][j] = MDBL_MIN;
472        }
473    }
474
475  free(expt);
476  free(uexpt);
477}
478
479/*********************************************************/
480
481void PMat(double l, model *mod, double ***Pij)
482{
483    switch(mod->datatype)
484        {
485        case NT :
486            {
487                if(mod->whichmodel < 3)
488                    PMat_K80(l,mod->kappa,Pij);
489                else
490                    {
491                        if(mod->whichmodel < 7)
492                            PMat_TN93(l,mod,Pij);
493                        else
494                            {
495                                PMat_Empirical(l,mod,Pij);
496                            }
497                    }
498                break;
499            }
500        case AA :
501            PMat_Empirical(l,mod,Pij);
502        }
503}
504
505/*********************************************************/
506
507void dPMat(double l, double rr, model *mod, double ***dPij)
508{
509  if(mod->whichmodel < 3)
510    dPMat_K80(l,dPij,rr,mod->kappa);
511  else
512    dPMat_TN93(l,dPij,mod,rr);
513}
514
515/*********************************************************/
516
517void d2PMat(double l, double rr, model *mod, double ***d2Pij)
518{
519  if(mod->whichmodel < 3)
520    d2PMat_K80(l,d2Pij,rr,mod->kappa);
521  else
522    d2PMat_TN93(l,d2Pij,mod,rr);
523}
524
525/*********************************************************/
526
527int GetDaa (double *daa, double *pi, char *file_name)
528{
529/* Get the amino acid distance (or substitution rate) matrix
530   (grantham, dayhoff, jones, etc).
531*/
532   FILE * fdaa;
533   int i,j, naa;
534   double dmax,dmin;
535   double sum; 
536
537   naa = 20;
538   dmax = .0;
539   dmin = 1.E+40;
540
541   fdaa = (FILE *)Openfile(file_name,0);
542
543   for (i=0; i<naa; i++)  for (j=0; j<i; j++)  {
544     fscanf(fdaa, "%lf", &daa[i*naa+j]);
545      daa[j*naa+i]=daa[i*naa+j];
546      if (dmax<daa[i*naa+j]) dmax=daa[i*naa+j];
547      if (dmin>daa[i*naa+j]) dmin=daa[i*naa+j];
548   }
549
550   For(i,naa) {
551     if(fscanf(fdaa,"%lf",&pi[i])!=1) Exit("err aaRatefile");
552   }
553   sum = 0.0;
554   For(i, naa) sum += pi[i];
555   if (fabs(1-sum)>1e-4) {
556     printf("\nSum of freq. = %.6f != 1 in aaRateFile\n",sum); 
557     exit(-1);
558   }
559   
560   fclose (fdaa);
561
562   return (0);
563}
564
565/*********************************************************/
566
567int Init_Qmat_Dayhoff(double *daa, double *pi)
568{
569  /* Dayhoff's model data
570   * Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. (1978)
571   * "A model of evolutionary change in proteins."
572   * Dayhoff, M.O.(ed.) Atlas of Protein Sequence Structur., Vol5, Suppl3.
573   * National Biomedical Research Foundation, Washington DC, pp.345-352.
574   */
575
576  int i,j,naa;
577
578  naa = 20;
579
580
581  daa[ 1*20+ 0] =   27.00; daa[ 2*20+ 0] =   98.00; daa[ 2*20+ 1] =   32.00; daa[ 3*20+ 0] =  120.00;
582  daa[ 3*20+ 1] =    0.00; daa[ 3*20+ 2] =  905.00; daa[ 4*20+ 0] =   36.00; daa[ 4*20+ 1] =   23.00;
583  daa[ 4*20+ 2] =    0.00; daa[ 4*20+ 3] =    0.00; daa[ 5*20+ 0] =   89.00; daa[ 5*20+ 1] =  246.00;
584  daa[ 5*20+ 2] =  103.00; daa[ 5*20+ 3] =  134.00; daa[ 5*20+ 4] =    0.00; daa[ 6*20+ 0] =  198.00;
585  daa[ 6*20+ 1] =    1.00; daa[ 6*20+ 2] =  148.00; daa[ 6*20+ 3] = 1153.00; daa[ 6*20+ 4] =    0.00;
586  daa[ 6*20+ 5] =  716.00; daa[ 7*20+ 0] =  240.00; daa[ 7*20+ 1] =    9.00; daa[ 7*20+ 2] =  139.00;
587  daa[ 7*20+ 3] =  125.00; daa[ 7*20+ 4] =   11.00; daa[ 7*20+ 5] =   28.00; daa[ 7*20+ 6] =   81.00;
588  daa[ 8*20+ 0] =   23.00; daa[ 8*20+ 1] =  240.00; daa[ 8*20+ 2] =  535.00; daa[ 8*20+ 3] =   86.00;
589  daa[ 8*20+ 4] =   28.00; daa[ 8*20+ 5] =  606.00; daa[ 8*20+ 6] =   43.00; daa[ 8*20+ 7] =   10.00;
590  daa[ 9*20+ 0] =   65.00; daa[ 9*20+ 1] =   64.00; daa[ 9*20+ 2] =   77.00; daa[ 9*20+ 3] =   24.00;
591  daa[ 9*20+ 4] =   44.00; daa[ 9*20+ 5] =   18.00; daa[ 9*20+ 6] =   61.00; daa[ 9*20+ 7] =    0.00;
592  daa[ 9*20+ 8] =    7.00; daa[10*20+ 0] =   41.00; daa[10*20+ 1] =   15.00; daa[10*20+ 2] =   34.00;
593  daa[10*20+ 3] =    0.00; daa[10*20+ 4] =    0.00; daa[10*20+ 5] =   73.00; daa[10*20+ 6] =   11.00;
594  daa[10*20+ 7] =    7.00; daa[10*20+ 8] =   44.00; daa[10*20+ 9] =  257.00; daa[11*20+ 0] =   26.00;
595  daa[11*20+ 1] =  464.00; daa[11*20+ 2] =  318.00; daa[11*20+ 3] =   71.00; daa[11*20+ 4] =    0.00;
596  daa[11*20+ 5] =  153.00; daa[11*20+ 6] =   83.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   26.00;
597  daa[11*20+ 9] =   46.00; daa[11*20+10] =   18.00; daa[12*20+ 0] =   72.00; daa[12*20+ 1] =   90.00;
598  daa[12*20+ 2] =    1.00; daa[12*20+ 3] =    0.00; daa[12*20+ 4] =    0.00; daa[12*20+ 5] =  114.00;
599  daa[12*20+ 6] =   30.00; daa[12*20+ 7] =   17.00; daa[12*20+ 8] =    0.00; daa[12*20+ 9] =  336.00;
600  daa[12*20+10] =  527.00; daa[12*20+11] =  243.00; daa[13*20+ 0] =   18.00; daa[13*20+ 1] =   14.00;
601  daa[13*20+ 2] =   14.00; daa[13*20+ 3] =    0.00; daa[13*20+ 4] =    0.00; daa[13*20+ 5] =    0.00;
602  daa[13*20+ 6] =    0.00; daa[13*20+ 7] =   15.00; daa[13*20+ 8] =   48.00; daa[13*20+ 9] =  196.00;
603  daa[13*20+10] =  157.00; daa[13*20+11] =    0.00; daa[13*20+12] =   92.00; daa[14*20+ 0] =  250.00;
604  daa[14*20+ 1] =  103.00; daa[14*20+ 2] =   42.00; daa[14*20+ 3] =   13.00; daa[14*20+ 4] =   19.00;
605  daa[14*20+ 5] =  153.00; daa[14*20+ 6] =   51.00; daa[14*20+ 7] =   34.00; daa[14*20+ 8] =   94.00;
606  daa[14*20+ 9] =   12.00; daa[14*20+10] =   32.00; daa[14*20+11] =   33.00; daa[14*20+12] =   17.00;
607  daa[14*20+13] =   11.00; daa[15*20+ 0] =  409.00; daa[15*20+ 1] =  154.00; daa[15*20+ 2] =  495.00;
608  daa[15*20+ 3] =   95.00; daa[15*20+ 4] =  161.00; daa[15*20+ 5] =   56.00; daa[15*20+ 6] =   79.00;
609  daa[15*20+ 7] =  234.00; daa[15*20+ 8] =   35.00; daa[15*20+ 9] =   24.00; daa[15*20+10] =   17.00;
610  daa[15*20+11] =   96.00; daa[15*20+12] =   62.00; daa[15*20+13] =   46.00; daa[15*20+14] =  245.00;
611  daa[16*20+ 0] =  371.00; daa[16*20+ 1] =   26.00; daa[16*20+ 2] =  229.00; daa[16*20+ 3] =   66.00;
612  daa[16*20+ 4] =   16.00; daa[16*20+ 5] =   53.00; daa[16*20+ 6] =   34.00; daa[16*20+ 7] =   30.00;
613  daa[16*20+ 8] =   22.00; daa[16*20+ 9] =  192.00; daa[16*20+10] =   33.00; daa[16*20+11] =  136.00;
614  daa[16*20+12] =  104.00; daa[16*20+13] =   13.00; daa[16*20+14] =   78.00; daa[16*20+15] =  550.00;
615  daa[17*20+ 0] =    0.00; daa[17*20+ 1] =  201.00; daa[17*20+ 2] =   23.00; daa[17*20+ 3] =    0.00;
616  daa[17*20+ 4] =    0.00; daa[17*20+ 5] =    0.00; daa[17*20+ 6] =    0.00; daa[17*20+ 7] =    0.00;
617  daa[17*20+ 8] =   27.00; daa[17*20+ 9] =    0.00; daa[17*20+10] =   46.00; daa[17*20+11] =    0.00;
618  daa[17*20+12] =    0.00; daa[17*20+13] =   76.00; daa[17*20+14] =    0.00; daa[17*20+15] =   75.00;
619  daa[17*20+16] =    0.00; daa[18*20+ 0] =   24.00; daa[18*20+ 1] =    8.00; daa[18*20+ 2] =   95.00;
620  daa[18*20+ 3] =    0.00; daa[18*20+ 4] =   96.00; daa[18*20+ 5] =    0.00; daa[18*20+ 6] =   22.00;
621  daa[18*20+ 7] =    0.00; daa[18*20+ 8] =  127.00; daa[18*20+ 9] =   37.00; daa[18*20+10] =   28.00;
622  daa[18*20+11] =   13.00; daa[18*20+12] =    0.00; daa[18*20+13] =  698.00; daa[18*20+14] =    0.00;
623  daa[18*20+15] =   34.00; daa[18*20+16] =   42.00; daa[18*20+17] =   61.00; daa[19*20+ 0] =  208.00;
624  daa[19*20+ 1] =   24.00; daa[19*20+ 2] =   15.00; daa[19*20+ 3] =   18.00; daa[19*20+ 4] =   49.00;
625  daa[19*20+ 5] =   35.00; daa[19*20+ 6] =   37.00; daa[19*20+ 7] =   54.00; daa[19*20+ 8] =   44.00;
626  daa[19*20+ 9] =  889.00; daa[19*20+10] =  175.00; daa[19*20+11] =   10.00; daa[19*20+12] =  258.00;
627  daa[19*20+13] =   12.00; daa[19*20+14] =   48.00; daa[19*20+15] =   30.00; daa[19*20+16] =  157.00;
628  daa[19*20+17] =    0.00; daa[19*20+18] =   28.00;
629 
630  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
631 
632  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872;
633  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
634  pi[ 8] = 0.033618; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080482;
635  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
636  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718;
637 
638 
639   
640  return 1;
641}
642
643/*********************************************************/
644
645int Init_Qmat_DCMut(double *daa, double *pi)
646{
647  /*
648     DCMut : new implementation based on Dayhoff et al.'s raw data and amino acid mutabilities
649     C. Kosiol and N. Goldman
650     http://www.ebi.ac.uk/goldman-srv/dayhoff/
651  */
652
653  int i,j,naa;
654
655  naa = 20;
656
657  daa[ 1*20+ 0] =   26.78280; daa[ 2*20+ 0] =   98.44740; daa[ 2*20+ 1] =   32.70590; daa[ 3*20+ 0] =  119.98050; 
658  daa[ 3*20+ 1] =    0.00000; daa[ 3*20+ 2] =  893.15150; daa[ 4*20+ 0] =   36.00160; daa[ 4*20+ 1] =   23.23740; 
659  daa[ 4*20+ 2] =    0.00000; daa[ 4*20+ 3] =    0.00000; daa[ 5*20+ 0] =   88.77530; daa[ 5*20+ 1] =  243.99390; 
660  daa[ 5*20+ 2] =  102.85090; daa[ 5*20+ 3] =  134.85510; daa[ 5*20+ 4] =    0.00000; daa[ 6*20+ 0] =  196.11670; 
661  daa[ 6*20+ 1] =    0.00000; daa[ 6*20+ 2] =  149.34090; daa[ 6*20+ 3] = 1138.86590; daa[ 6*20+ 4] =    0.00000; 
662  daa[ 6*20+ 5] =  708.60220; daa[ 7*20+ 0] =  238.61110; daa[ 7*20+ 1] =    8.77910; daa[ 7*20+ 2] =  138.53520; 
663  daa[ 7*20+ 3] =  124.09810; daa[ 7*20+ 4] =   10.72780; daa[ 7*20+ 5] =   28.15810; daa[ 7*20+ 6] =   81.19070; 
664  daa[ 8*20+ 0] =   22.81160; daa[ 8*20+ 1] =  238.31480; daa[ 8*20+ 2] =  529.00240; daa[ 8*20+ 3] =   86.82410; 
665  daa[ 8*20+ 4] =   28.27290; daa[ 8*20+ 5] =  601.16130; daa[ 8*20+ 6] =   43.94690; daa[ 8*20+ 7] =   10.68020; 
666  daa[ 9*20+ 0] =   65.34160; daa[ 9*20+ 1] =   63.26290; daa[ 9*20+ 2] =   76.80240; daa[ 9*20+ 3] =   23.92480; 
667  daa[ 9*20+ 4] =   43.80740; daa[ 9*20+ 5] =   18.03930; daa[ 9*20+ 6] =   60.95260; daa[ 9*20+ 7] =    0.00000; 
668  daa[ 9*20+ 8] =    7.69810; daa[10*20+ 0] =   40.64310; daa[10*20+ 1] =   15.49240; daa[10*20+ 2] =   34.11130; 
669  daa[10*20+ 3] =    0.00000; daa[10*20+ 4] =    0.00000; daa[10*20+ 5] =   73.07720; daa[10*20+ 6] =   11.28800; 
670  daa[10*20+ 7] =    7.15140; daa[10*20+ 8] =   44.35040; daa[10*20+ 9] =  255.66850; daa[11*20+ 0] =   25.86350; 
671  daa[11*20+ 1] =  461.01240; daa[11*20+ 2] =  314.83710; daa[11*20+ 3] =   71.69130; daa[11*20+ 4] =    0.00000; 
672  daa[11*20+ 5] =  151.90780; daa[11*20+ 6] =   83.00780; daa[11*20+ 7] =   26.76830; daa[11*20+ 8] =   27.04750; 
673  daa[11*20+ 9] =   46.08570; daa[11*20+10] =   18.06290; daa[12*20+ 0] =   71.78400; daa[12*20+ 1] =   89.63210; 
674  daa[12*20+ 2] =    0.00000; daa[12*20+ 3] =    0.00000; daa[12*20+ 4] =    0.00000; daa[12*20+ 5] =  112.74990; 
675  daa[12*20+ 6] =   30.48030; daa[12*20+ 7] =   17.03720; daa[12*20+ 8] =    0.00000; daa[12*20+ 9] =  333.27320; 
676  daa[12*20+10] =  523.01150; daa[12*20+11] =  241.17390; daa[13*20+ 0] =   18.36410; daa[13*20+ 1] =   13.69060; 
677  daa[13*20+ 2] =   13.85030; daa[13*20+ 3] =    0.00000; daa[13*20+ 4] =    0.00000; daa[13*20+ 5] =    0.00000; 
678  daa[13*20+ 6] =    0.00000; daa[13*20+ 7] =   15.34780; daa[13*20+ 8] =   47.59270; daa[13*20+ 9] =  195.19510; 
679  daa[13*20+10] =  156.51600; daa[13*20+11] =    0.00000; daa[13*20+12] =   92.18600; daa[14*20+ 0] =  248.59200; 
680  daa[14*20+ 1] =  102.83130; daa[14*20+ 2] =   41.92440; daa[14*20+ 3] =   13.39400; daa[14*20+ 4] =   18.75500; 
681  daa[14*20+ 5] =  152.61880; daa[14*20+ 6] =   50.70030; daa[14*20+ 7] =   34.71530; daa[14*20+ 8] =   93.37090; 
682  daa[14*20+ 9] =   11.91520; daa[14*20+10] =   31.62580; daa[14*20+11] =   33.54190; daa[14*20+12] =   17.02050; 
683  daa[14*20+13] =   11.05060; daa[15*20+ 0] =  405.18700; daa[15*20+ 1] =  153.15900; daa[15*20+ 2] =  488.58920; 
684  daa[15*20+ 3] =   95.60970; daa[15*20+ 4] =  159.83560; daa[15*20+ 5] =   56.18280; daa[15*20+ 6] =   79.39990; 
685  daa[15*20+ 7] =  232.22430; daa[15*20+ 8] =   35.36430; daa[15*20+ 9] =   24.79550; daa[15*20+10] =   17.14320; 
686  daa[15*20+11] =   95.45570; daa[15*20+12] =   61.99510; daa[15*20+13] =   45.99010; daa[15*20+14] =  242.72020; 
687  daa[16*20+ 0] =  368.03650; daa[16*20+ 1] =   26.57450; daa[16*20+ 2] =  227.16970; daa[16*20+ 3] =   66.09300; 
688  daa[16*20+ 4] =   16.23660; daa[16*20+ 5] =   52.56510; daa[16*20+ 6] =   34.01560; daa[16*20+ 7] =   30.66620; 
689  daa[16*20+ 8] =   22.63330; daa[16*20+ 9] =  190.07390; daa[16*20+10] =   33.10900; daa[16*20+11] =  135.05990; 
690  daa[16*20+12] =  103.15340; daa[16*20+13] =   13.66550; daa[16*20+14] =   78.28570; daa[16*20+15] =  543.66740; 
691  daa[17*20+ 0] =    0.00000; daa[17*20+ 1] =  200.13750; daa[17*20+ 2] =   22.49680; daa[17*20+ 3] =    0.00000; 
692  daa[17*20+ 4] =    0.00000; daa[17*20+ 5] =    0.00000; daa[17*20+ 6] =    0.00000; daa[17*20+ 7] =    0.00000; 
693  daa[17*20+ 8] =   27.05640; daa[17*20+ 9] =    0.00000; daa[17*20+10] =   46.17760; daa[17*20+11] =    0.00000; 
694  daa[17*20+12] =    0.00000; daa[17*20+13] =   76.23540; daa[17*20+14] =    0.00000; daa[17*20+15] =   74.08190; 
695  daa[17*20+16] =    0.00000; daa[18*20+ 0] =   24.41390; daa[18*20+ 1] =    7.80120; daa[18*20+ 2] =   94.69400; 
696  daa[18*20+ 3] =    0.00000; daa[18*20+ 4] =   95.31640; daa[18*20+ 5] =    0.00000; daa[18*20+ 6] =   21.47170; 
697  daa[18*20+ 7] =    0.00000; daa[18*20+ 8] =  126.54000; daa[18*20+ 9] =   37.48340; daa[18*20+10] =   28.65720; 
698  daa[18*20+11] =   13.21420; daa[18*20+12] =    0.00000; daa[18*20+13] =  695.26290; daa[18*20+14] =    0.00000; 
699  daa[18*20+15] =   33.62890; daa[18*20+16] =   41.78390; daa[18*20+17] =   60.80700; daa[19*20+ 0] =  205.95640; 
700  daa[19*20+ 1] =   24.03680; daa[19*20+ 2] =   15.80670; daa[19*20+ 3] =   17.83160; daa[19*20+ 4] =   48.46780; 
701  daa[19*20+ 5] =   34.69830; daa[19*20+ 6] =   36.72500; daa[19*20+ 7] =   53.81650; daa[19*20+ 8] =   43.87150; 
702  daa[19*20+ 9] =  881.00380; daa[19*20+10] =  174.51560; daa[19*20+11] =   10.38500; daa[19*20+12] =  256.59550; 
703  daa[19*20+13] =   12.36060; daa[19*20+14] =   48.50260; daa[19*20+15] =   30.38360; daa[19*20+16] =  156.19970; 
704  daa[19*20+17] =    0.00000; daa[19*20+18] =   27.93790; 
705 
706
707
708  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
709
710  pi[ 0] = 0.087127; pi[ 1] = 0.040904; pi[ 2] = 0.040432; pi[ 3] = 0.046872; 
711  pi[ 4] = 0.033474; pi[ 5] = 0.038255; pi[ 6] = 0.049530; pi[ 7] = 0.088612;
712  pi[ 8] = 0.033619; pi[ 9] = 0.036886; pi[10] = 0.085357; pi[11] = 0.080481;
713  pi[12] = 0.014753; pi[13] = 0.039772; pi[14] = 0.050680; pi[15] = 0.069577;
714  pi[16] = 0.058542; pi[17] = 0.010494; pi[18] = 0.029916; pi[19] = 0.064718; 
715
716  return 1;
717
718
719}
720
721
722/*********************************************************/
723
724int Init_Qmat_JTT(double *daa, double *pi)
725{
726  int i,j,naa;
727
728  /* JTT's model data
729   * D.T.Jones, W.R.Taylor and J.M.Thornton
730   * "The rapid generation of mutation data matrices from protein sequences"
731   * CABIOS  vol.8 no.3 1992 pp275-282
732   */
733
734
735  naa = 20;
736
737  daa[ 1*20+ 0] =   58.00; daa[ 2*20+ 0] =   54.00; daa[ 2*20+ 1] =   45.00; daa[ 3*20+ 0] =   81.00;
738  daa[ 3*20+ 1] =   16.00; daa[ 3*20+ 2] =  528.00; daa[ 4*20+ 0] =   56.00; daa[ 4*20+ 1] =  113.00;
739  daa[ 4*20+ 2] =   34.00; daa[ 4*20+ 3] =   10.00; daa[ 5*20+ 0] =   57.00; daa[ 5*20+ 1] =  310.00;
740  daa[ 5*20+ 2] =   86.00; daa[ 5*20+ 3] =   49.00; daa[ 5*20+ 4] =    9.00; daa[ 6*20+ 0] =  105.00;
741  daa[ 6*20+ 1] =   29.00; daa[ 6*20+ 2] =   58.00; daa[ 6*20+ 3] =  767.00; daa[ 6*20+ 4] =    5.00;
742  daa[ 6*20+ 5] =  323.00; daa[ 7*20+ 0] =  179.00; daa[ 7*20+ 1] =  137.00; daa[ 7*20+ 2] =   81.00;
743  daa[ 7*20+ 3] =  130.00; daa[ 7*20+ 4] =   59.00; daa[ 7*20+ 5] =   26.00; daa[ 7*20+ 6] =  119.00;
744  daa[ 8*20+ 0] =   27.00; daa[ 8*20+ 1] =  328.00; daa[ 8*20+ 2] =  391.00; daa[ 8*20+ 3] =  112.00;
745  daa[ 8*20+ 4] =   69.00; daa[ 8*20+ 5] =  597.00; daa[ 8*20+ 6] =   26.00; daa[ 8*20+ 7] =   23.00;
746  daa[ 9*20+ 0] =   36.00; daa[ 9*20+ 1] =   22.00; daa[ 9*20+ 2] =   47.00; daa[ 9*20+ 3] =   11.00;
747  daa[ 9*20+ 4] =   17.00; daa[ 9*20+ 5] =    9.00; daa[ 9*20+ 6] =   12.00; daa[ 9*20+ 7] =    6.00;
748  daa[ 9*20+ 8] =   16.00; daa[10*20+ 0] =   30.00; daa[10*20+ 1] =   38.00; daa[10*20+ 2] =   12.00;
749  daa[10*20+ 3] =    7.00; daa[10*20+ 4] =   23.00; daa[10*20+ 5] =   72.00; daa[10*20+ 6] =    9.00;
750  daa[10*20+ 7] =    6.00; daa[10*20+ 8] =   56.00; daa[10*20+ 9] =  229.00; daa[11*20+ 0] =   35.00;
751  daa[11*20+ 1] =  646.00; daa[11*20+ 2] =  263.00; daa[11*20+ 3] =   26.00; daa[11*20+ 4] =    7.00;
752  daa[11*20+ 5] =  292.00; daa[11*20+ 6] =  181.00; daa[11*20+ 7] =   27.00; daa[11*20+ 8] =   45.00;
753  daa[11*20+ 9] =   21.00; daa[11*20+10] =   14.00; daa[12*20+ 0] =   54.00; daa[12*20+ 1] =   44.00;
754  daa[12*20+ 2] =   30.00; daa[12*20+ 3] =   15.00; daa[12*20+ 4] =   31.00; daa[12*20+ 5] =   43.00;
755  daa[12*20+ 6] =   18.00; daa[12*20+ 7] =   14.00; daa[12*20+ 8] =   33.00; daa[12*20+ 9] =  479.00;
756  daa[12*20+10] =  388.00; daa[12*20+11] =   65.00; daa[13*20+ 0] =   15.00; daa[13*20+ 1] =    5.00;
757  daa[13*20+ 2] =   10.00; daa[13*20+ 3] =    4.00; daa[13*20+ 4] =   78.00; daa[13*20+ 5] =    4.00;
758  daa[13*20+ 6] =    5.00; daa[13*20+ 7] =    5.00; daa[13*20+ 8] =   40.00; daa[13*20+ 9] =   89.00;
759  daa[13*20+10] =  248.00; daa[13*20+11] =    4.00; daa[13*20+12] =   43.00; daa[14*20+ 0] =  194.00;
760  daa[14*20+ 1] =   74.00; daa[14*20+ 2] =   15.00; daa[14*20+ 3] =   15.00; daa[14*20+ 4] =   14.00;
761  daa[14*20+ 5] =  164.00; daa[14*20+ 6] =   18.00; daa[14*20+ 7] =   24.00; daa[14*20+ 8] =  115.00;
762  daa[14*20+ 9] =   10.00; daa[14*20+10] =  102.00; daa[14*20+11] =   21.00; daa[14*20+12] =   16.00;
763  daa[14*20+13] =   17.00; daa[15*20+ 0] =  378.00; daa[15*20+ 1] =  101.00; daa[15*20+ 2] =  503.00;
764  daa[15*20+ 3] =   59.00; daa[15*20+ 4] =  223.00; daa[15*20+ 5] =   53.00; daa[15*20+ 6] =   30.00;
765  daa[15*20+ 7] =  201.00; daa[15*20+ 8] =   73.00; daa[15*20+ 9] =   40.00; daa[15*20+10] =   59.00;
766  daa[15*20+11] =   47.00; daa[15*20+12] =   29.00; daa[15*20+13] =   92.00; daa[15*20+14] =  285.00;
767  daa[16*20+ 0] =  475.00; daa[16*20+ 1] =   64.00; daa[16*20+ 2] =  232.00; daa[16*20+ 3] =   38.00;
768  daa[16*20+ 4] =   42.00; daa[16*20+ 5] =   51.00; daa[16*20+ 6] =   32.00; daa[16*20+ 7] =   33.00;
769  daa[16*20+ 8] =   46.00; daa[16*20+ 9] =  245.00; daa[16*20+10] =   25.00; daa[16*20+11] =  103.00;
770  daa[16*20+12] =  226.00; daa[16*20+13] =   12.00; daa[16*20+14] =  118.00; daa[16*20+15] =  477.00;
771  daa[17*20+ 0] =    9.00; daa[17*20+ 1] =  126.00; daa[17*20+ 2] =    8.00; daa[17*20+ 3] =    4.00;
772  daa[17*20+ 4] =  115.00; daa[17*20+ 5] =   18.00; daa[17*20+ 6] =   10.00; daa[17*20+ 7] =   55.00;
773  daa[17*20+ 8] =    8.00; daa[17*20+ 9] =    9.00; daa[17*20+10] =   52.00; daa[17*20+11] =   10.00;
774  daa[17*20+12] =   24.00; daa[17*20+13] =   53.00; daa[17*20+14] =    6.00; daa[17*20+15] =   35.00;
775  daa[17*20+16] =   12.00; daa[18*20+ 0] =   11.00; daa[18*20+ 1] =   20.00; daa[18*20+ 2] =   70.00;
776  daa[18*20+ 3] =   46.00; daa[18*20+ 4] =  209.00; daa[18*20+ 5] =   24.00; daa[18*20+ 6] =    7.00;
777  daa[18*20+ 7] =    8.00; daa[18*20+ 8] =  573.00; daa[18*20+ 9] =   32.00; daa[18*20+10] =   24.00;
778  daa[18*20+11] =    8.00; daa[18*20+12] =   18.00; daa[18*20+13] =  536.00; daa[18*20+14] =   10.00;
779  daa[18*20+15] =   63.00; daa[18*20+16] =   21.00; daa[18*20+17] =   71.00; daa[19*20+ 0] =  298.00;
780  daa[19*20+ 1] =   17.00; daa[19*20+ 2] =   16.00; daa[19*20+ 3] =   31.00; daa[19*20+ 4] =   62.00;
781  daa[19*20+ 5] =   20.00; daa[19*20+ 6] =   45.00; daa[19*20+ 7] =   47.00; daa[19*20+ 8] =   11.00;
782  daa[19*20+ 9] =  961.00; daa[19*20+10] =  180.00; daa[19*20+11] =   14.00; daa[19*20+12] =  323.00;
783  daa[19*20+13] =   62.00; daa[19*20+14] =   23.00; daa[19*20+15] =   38.00; daa[19*20+16] =  112.00;
784  daa[19*20+17] =   25.00; daa[19*20+18] =   16.00;
785 
786  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
787
788  pi[ 0] = 0.076748; pi[ 1] = 0.051691; pi[ 2] = 0.042645; pi[ 3] = 0.051544;
789  pi[ 4] = 0.019803; pi[ 5] = 0.040752; pi[ 6] = 0.061830; pi[ 7] = 0.073152;
790  pi[ 8] = 0.022944; pi[ 9] = 0.053761; pi[10] = 0.091904; pi[11] = 0.058676;
791  pi[12] = 0.023826; pi[13] = 0.040126; pi[14] = 0.050901; pi[15] = 0.068765;
792  pi[16] = 0.058565; pi[17] = 0.014261; pi[18] = 0.032102; pi[19] = 0.066005;
793 
794 return 1;
795}
796
797/*********************************************************/
798
799int Init_Qmat_MtREV(double *daa, double *pi)
800{
801  int i,j,naa;
802
803  naa = 20;
804
805
806  daa[ 1*20+ 0] =   23.18; daa[ 2*20+ 0] =   26.95; daa[ 2*20+ 1] =   13.24; daa[ 3*20+ 0] =   17.67;
807  daa[ 3*20+ 1] =    1.90; daa[ 3*20+ 2] =  794.38; daa[ 4*20+ 0] =   59.93; daa[ 4*20+ 1] =  103.33;
808  daa[ 4*20+ 2] =   58.94; daa[ 4*20+ 3] =    1.90; daa[ 5*20+ 0] =    1.90; daa[ 5*20+ 1] =  220.99;
809  daa[ 5*20+ 2] =  173.56; daa[ 5*20+ 3] =   55.28; daa[ 5*20+ 4] =   75.24; daa[ 6*20+ 0] =    9.77;
810  daa[ 6*20+ 1] =    1.90; daa[ 6*20+ 2] =   63.05; daa[ 6*20+ 3] =  583.55; daa[ 6*20+ 4] =    1.90;
811  daa[ 6*20+ 5] =  313.56; daa[ 7*20+ 0] =  120.71; daa[ 7*20+ 1] =   23.03; daa[ 7*20+ 2] =   53.30;
812  daa[ 7*20+ 3] =   56.77; daa[ 7*20+ 4] =   30.71; daa[ 7*20+ 5] =    6.75; daa[ 7*20+ 6] =   28.28;
813  daa[ 8*20+ 0] =   13.90; daa[ 8*20+ 1] =  165.23; daa[ 8*20+ 2] =  496.13; daa[ 8*20+ 3] =  113.99;
814  daa[ 8*20+ 4] =  141.49; daa[ 8*20+ 5] =  582.40; daa[ 8*20+ 6] =   49.12; daa[ 8*20+ 7] =    1.90;
815  daa[ 9*20+ 0] =   96.49; daa[ 9*20+ 1] =    1.90; daa[ 9*20+ 2] =   27.10; daa[ 9*20+ 3] =    4.34;
816  daa[ 9*20+ 4] =   62.73; daa[ 9*20+ 5] =    8.34; daa[ 9*20+ 6] =    3.31; daa[ 9*20+ 7] =    5.98;
817  daa[ 9*20+ 8] =   12.26; daa[10*20+ 0] =   25.46; daa[10*20+ 1] =   15.58; daa[10*20+ 2] =   15.16;
818  daa[10*20+ 3] =    1.90; daa[10*20+ 4] =   25.65; daa[10*20+ 5] =   39.70; daa[10*20+ 6] =    1.90;
819  daa[10*20+ 7] =    2.41; daa[10*20+ 8] =   11.49; daa[10*20+ 9] =  329.09; daa[11*20+ 0] =    8.36;
820  daa[11*20+ 1] =  141.40; daa[11*20+ 2] =  608.70; daa[11*20+ 3] =    2.31; daa[11*20+ 4] =    1.90;
821  daa[11*20+ 5] =  465.58; daa[11*20+ 6] =  313.86; daa[11*20+ 7] =   22.73; daa[11*20+ 8] =  127.67;
822  daa[11*20+ 9] =   19.57; daa[11*20+10] =   14.88; daa[12*20+ 0] =  141.88; daa[12*20+ 1] =    1.90;
823  daa[12*20+ 2] =   65.41; daa[12*20+ 3] =    1.90; daa[12*20+ 4] =    6.18; daa[12*20+ 5] =   47.37;
824  daa[12*20+ 6] =    1.90; daa[12*20+ 7] =    1.90; daa[12*20+ 8] =   11.97; daa[12*20+ 9] =  517.98;
825  daa[12*20+10] =  537.53; daa[12*20+11] =   91.37; daa[13*20+ 0] =    6.37; daa[13*20+ 1] =    4.69;
826  daa[13*20+ 2] =   15.20; daa[13*20+ 3] =    4.98; daa[13*20+ 4] =   70.80; daa[13*20+ 5] =   19.11;
827  daa[13*20+ 6] =    2.67; daa[13*20+ 7] =    1.90; daa[13*20+ 8] =   48.16; daa[13*20+ 9] =   84.67;
828  daa[13*20+10] =  216.06; daa[13*20+11] =    6.44; daa[13*20+12] =   90.82; daa[14*20+ 0] =   54.31;
829  daa[14*20+ 1] =   23.64; daa[14*20+ 2] =   73.31; daa[14*20+ 3] =   13.43; daa[14*20+ 4] =   31.26;
830  daa[14*20+ 5] =  137.29; daa[14*20+ 6] =   12.83; daa[14*20+ 7] =    1.90; daa[14*20+ 8] =   60.97;
831  daa[14*20+ 9] =   20.63; daa[14*20+10] =   40.10; daa[14*20+11] =   50.10; daa[14*20+12] =   18.84;
832  daa[14*20+13] =   17.31; daa[15*20+ 0] =  387.86; daa[15*20+ 1] =    6.04; daa[15*20+ 2] =  494.39;
833  daa[15*20+ 3] =   69.02; daa[15*20+ 4] =  277.05; daa[15*20+ 5] =   54.11; daa[15*20+ 6] =   54.71;
834  daa[15*20+ 7] =  125.93; daa[15*20+ 8] =   77.46; daa[15*20+ 9] =   47.70; daa[15*20+10] =   73.61;
835  daa[15*20+11] =  105.79; daa[15*20+12] =  111.16; daa[15*20+13] =   64.29; daa[15*20+14] =  169.90;
836  daa[16*20+ 0] =  480.72; daa[16*20+ 1] =    2.08; daa[16*20+ 2] =  238.46; daa[16*20+ 3] =   28.01;
837  daa[16*20+ 4] =  179.97; daa[16*20+ 5] =   94.93; daa[16*20+ 6] =   14.82; daa[16*20+ 7] =   11.17;
838  daa[16*20+ 8] =   44.78; daa[16*20+ 9] =  368.43; daa[16*20+10] =  126.40; daa[16*20+11] =  136.33;
839  daa[16*20+12] =  528.17; daa[16*20+13] =   33.85; daa[16*20+14] =  128.22; daa[16*20+15] =  597.21;
840  daa[17*20+ 0] =    1.90; daa[17*20+ 1] =   21.95; daa[17*20+ 2] =   10.68; daa[17*20+ 3] =   19.86;
841  daa[17*20+ 4] =   33.60; daa[17*20+ 5] =    1.90; daa[17*20+ 6] =    1.90; daa[17*20+ 7] =   10.92;
842  daa[17*20+ 8] =    7.08; daa[17*20+ 9] =    1.90; daa[17*20+10] =   32.44; daa[17*20+11] =   24.00;
843  daa[17*20+12] =   21.71; daa[17*20+13] =    7.84; daa[17*20+14] =    4.21; daa[17*20+15] =   38.58;
844  daa[17*20+16] =    9.99; daa[18*20+ 0] =    6.48; daa[18*20+ 1] =    1.90; daa[18*20+ 2] =  191.36;
845  daa[18*20+ 3] =   21.21; daa[18*20+ 4] =  254.77; daa[18*20+ 5] =   38.82; daa[18*20+ 6] =   13.12;
846  daa[18*20+ 7] =    3.21; daa[18*20+ 8] =  670.14; daa[18*20+ 9] =   25.01; daa[18*20+10] =   44.15;
847  daa[18*20+11] =   51.17; daa[18*20+12] =   39.96; daa[18*20+13] =  465.58; daa[18*20+14] =   16.21;
848  daa[18*20+15] =   64.92; daa[18*20+16] =   38.73; daa[18*20+17] =   26.25; daa[19*20+ 0] =  195.06;
849  daa[19*20+ 1] =    7.64; daa[19*20+ 2] =    1.90; daa[19*20+ 3] =    1.90; daa[19*20+ 4] =    1.90;
850  daa[19*20+ 5] =   19.00; daa[19*20+ 6] =   21.14; daa[19*20+ 7] =    2.53; daa[19*20+ 8] =    1.90;
851  daa[19*20+ 9] = 1222.94; daa[19*20+10] =   91.67; daa[19*20+11] =    1.90; daa[19*20+12] =  387.54;
852  daa[19*20+13] =    6.35; daa[19*20+14] =    8.23; daa[19*20+15] =    1.90; daa[19*20+16] =  204.54;
853  daa[19*20+17] =    5.37; daa[19*20+18] =    1.90;
854 
855  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
856
857  pi[ 0] = 0.072000; pi[ 1] = 0.019000; pi[ 2] = 0.039000; pi[ 3] = 0.019000;
858  pi[ 4] = 0.006000; pi[ 5] = 0.025000; pi[ 6] = 0.024000; pi[ 7] = 0.056000;
859  pi[ 8] = 0.028000; pi[ 9] = 0.088000; pi[10] = 0.169000; pi[11] = 0.023000;
860  pi[12] = 0.054000; pi[13] = 0.061000; pi[14] = 0.054000; pi[15] = 0.072000;
861  pi[16] = 0.086000; pi[17] = 0.029000; pi[18] = 0.033000; pi[19] = 0.043000;
862
863  return 1;
864}
865
866
867/*********************************************************/
868
869int Init_Qmat_WAG(double *daa, double *pi)
870{
871  int i,j,naa;
872
873  /* WAG's model data
874   * Simon Whelan and Nick Goldman
875   * 'A general empirical model of protein evolution derived from multiple
876   *  protein families using a maximum-likelihood approach'
877   * MBE (2001) 18:691-699
878   */
879
880  naa = 20;
881
882  daa[ 1*20+ 0] =  55.15710; daa[ 2*20+ 0] =  50.98480; daa[ 2*20+ 1] =  63.53460; 
883  daa[ 3*20+ 0] =  73.89980; daa[ 3*20+ 1] =  14.73040; daa[ 3*20+ 2] = 542.94200; 
884  daa[ 4*20+ 0] = 102.70400; daa[ 4*20+ 1] =  52.81910; daa[ 4*20+ 2] =  26.52560; 
885  daa[ 4*20+ 3] =   3.02949; daa[ 5*20+ 0] =  90.85980; daa[ 5*20+ 1] = 303.55000; 
886  daa[ 5*20+ 2] = 154.36400; daa[ 5*20+ 3] =  61.67830; daa[ 5*20+ 4] =   9.88179; 
887  daa[ 6*20+ 0] = 158.28500; daa[ 6*20+ 1] =  43.91570; daa[ 6*20+ 2] =  94.71980; 
888  daa[ 6*20+ 3] = 617.41600; daa[ 6*20+ 4] =   2.13520; daa[ 6*20+ 5] = 546.94700; 
889  daa[ 7*20+ 0] = 141.67200; daa[ 7*20+ 1] =  58.46650; daa[ 7*20+ 2] = 112.55600; 
890  daa[ 7*20+ 3] =  86.55840; daa[ 7*20+ 4] =  30.66740; daa[ 7*20+ 5] =  33.00520; 
891  daa[ 7*20+ 6] =  56.77170; daa[ 8*20+ 0] =  31.69540; daa[ 8*20+ 1] = 213.71500; 
892  daa[ 8*20+ 2] = 395.62900; daa[ 8*20+ 3] =  93.06760; daa[ 8*20+ 4] =  24.89720; 
893  daa[ 8*20+ 5] = 429.41100; daa[ 8*20+ 6] =  57.00250; daa[ 8*20+ 7] =  24.94100; 
894  daa[ 9*20+ 0] =  19.33350; daa[ 9*20+ 1] =  18.69790; daa[ 9*20+ 2] =  55.42360; 
895  daa[ 9*20+ 3] =   3.94370; daa[ 9*20+ 4] =  17.01350; daa[ 9*20+ 5] =  11.39170; 
896  daa[ 9*20+ 6] =  12.73950; daa[ 9*20+ 7] =   3.04501; daa[ 9*20+ 8] =  13.81900; 
897  daa[10*20+ 0] =  39.79150; daa[10*20+ 1] =  49.76710; daa[10*20+ 2] =  13.15280; 
898  daa[10*20+ 3] =   8.48047; daa[10*20+ 4] =  38.42870; daa[10*20+ 5] =  86.94890; 
899  daa[10*20+ 6] =  15.42630; daa[10*20+ 7] =   6.13037; daa[10*20+ 8] =  49.94620; 
900  daa[10*20+ 9] = 317.09700; daa[11*20+ 0] =  90.62650; daa[11*20+ 1] = 535.14200; 
901  daa[11*20+ 2] = 301.20100; daa[11*20+ 3] =  47.98550; daa[11*20+ 4] =   7.40339; 
902  daa[11*20+ 5] = 389.49000; daa[11*20+ 6] = 258.44300; daa[11*20+ 7] =  37.35580; 
903  daa[11*20+ 8] =  89.04320; daa[11*20+ 9] =  32.38320; daa[11*20+10] =  25.75550; 
904  daa[12*20+ 0] =  89.34960; daa[12*20+ 1] =  68.31620; daa[12*20+ 2] =  19.82210; 
905  daa[12*20+ 3] =  10.37540; daa[12*20+ 4] =  39.04820; daa[12*20+ 5] = 154.52600; 
906  daa[12*20+ 6] =  31.51240; daa[12*20+ 7] =  17.41000; daa[12*20+ 8] =  40.41410; 
907  daa[12*20+ 9] = 425.74600; daa[12*20+10] = 485.40200; daa[12*20+11] =  93.42760; 
908  daa[13*20+ 0] =  21.04940; daa[13*20+ 1] =  10.27110; daa[13*20+ 2] =   9.61621; 
909  daa[13*20+ 3] =   4.67304; daa[13*20+ 4] =  39.80200; daa[13*20+ 5] =   9.99208; 
910  daa[13*20+ 6] =   8.11339; daa[13*20+ 7] =   4.99310; daa[13*20+ 8] =  67.93710; 
911  daa[13*20+ 9] = 105.94700; daa[13*20+10] = 211.51700; daa[13*20+11] =   8.88360; 
912  daa[13*20+12] = 119.06300; daa[14*20+ 0] = 143.85500; daa[14*20+ 1] =  67.94890; 
913  daa[14*20+ 2] =  19.50810; daa[14*20+ 3] =  42.39840; daa[14*20+ 4] =  10.94040; 
914  daa[14*20+ 5] =  93.33720; daa[14*20+ 6] =  68.23550; daa[14*20+ 7] =  24.35700; 
915  daa[14*20+ 8] =  69.61980; daa[14*20+ 9] =   9.99288; daa[14*20+10] =  41.58440; 
916  daa[14*20+11] =  55.68960; daa[14*20+12] =  17.13290; daa[14*20+13] =  16.14440; 
917  daa[15*20+ 0] = 337.07900; daa[15*20+ 1] = 122.41900; daa[15*20+ 2] = 397.42300; 
918  daa[15*20+ 3] = 107.17600; daa[15*20+ 4] = 140.76600; daa[15*20+ 5] = 102.88700; 
919  daa[15*20+ 6] =  70.49390; daa[15*20+ 7] = 134.18200; daa[15*20+ 8] =  74.01690; 
920  daa[15*20+ 9] =  31.94400; daa[15*20+10] =  34.47390; daa[15*20+11] =  96.71300; 
921  daa[15*20+12] =  49.39050; daa[15*20+13] =  54.59310; daa[15*20+14] = 161.32800; 
922  daa[16*20+ 0] = 212.11100; daa[16*20+ 1] =  55.44130; daa[16*20+ 2] = 203.00600; 
923  daa[16*20+ 3] =  37.48660; daa[16*20+ 4] =  51.29840; daa[16*20+ 5] =  85.79280; 
924  daa[16*20+ 6] =  82.27650; daa[16*20+ 7] =  22.58330; daa[16*20+ 8] =  47.33070; 
925  daa[16*20+ 9] = 145.81600; daa[16*20+10] =  32.66220; daa[16*20+11] = 138.69800; 
926  daa[16*20+12] = 151.61200; daa[16*20+13] =  17.19030; daa[16*20+14] =  79.53840; 
927  daa[16*20+15] = 437.80200; daa[17*20+ 0] =  11.31330; daa[17*20+ 1] = 116.39200; 
928  daa[17*20+ 2] =   7.19167; daa[17*20+ 3] =  12.97670; daa[17*20+ 4] =  71.70700; 
929  daa[17*20+ 5] =  21.57370; daa[17*20+ 6] =  15.65570; daa[17*20+ 7] =  33.69830; 
930  daa[17*20+ 8] =  26.25690; daa[17*20+ 9] =  21.24830; daa[17*20+10] =  66.53090; 
931  daa[17*20+11] =  13.75050; daa[17*20+12] =  51.57060; daa[17*20+13] = 152.96400; 
932  daa[17*20+14] =  13.94050; daa[17*20+15] =  52.37420; daa[17*20+16] =  11.08640; 
933  daa[18*20+ 0] =  24.07350; daa[18*20+ 1] =  38.15330; daa[18*20+ 2] = 108.60000; 
934  daa[18*20+ 3] =  32.57110; daa[18*20+ 4] =  54.38330; daa[18*20+ 5] =  22.77100; 
935  daa[18*20+ 6] =  19.63030; daa[18*20+ 7] =  10.36040; daa[18*20+ 8] = 387.34400; 
936  daa[18*20+ 9] =  42.01700; daa[18*20+10] =  39.86180; daa[18*20+11] =  13.32640; 
937  daa[18*20+12] =  42.84370; daa[18*20+13] = 645.42800; daa[18*20+14] =  21.60460; 
938  daa[18*20+15] =  78.69930; daa[18*20+16] =  29.11480; daa[18*20+17] = 248.53900; 
939  daa[19*20+ 0] = 200.60100; daa[19*20+ 1] =  25.18490; daa[19*20+ 2] =  19.62460; 
940  daa[19*20+ 3] =  15.23350; daa[19*20+ 4] = 100.21400; daa[19*20+ 5] =  30.12810; 
941  daa[19*20+ 6] =  58.87310; daa[19*20+ 7] =  18.72470; daa[19*20+ 8] =  11.83580; 
942  daa[19*20+ 9] = 782.13000; daa[19*20+10] = 180.03400; daa[19*20+11] =  30.54340; 
943  daa[19*20+12] = 205.84500; daa[19*20+13] =  64.98920; daa[19*20+14] =  31.48870; 
944  daa[19*20+15] =  23.27390; daa[19*20+16] = 138.82300; daa[19*20+17] =  36.53690; 
945  daa[19*20+18] =  31.47300; 
946
947  for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
948
949  pi[0] = 0.0866279; pi[1] =  0.043972; pi[2] =  0.0390894; pi[3] =  0.0570451;
950  pi[4] =  0.0193078; pi[5] =  0.0367281; pi[6] =  0.0580589; pi[7] =  0.0832518;
951  pi[8] =  0.0244313; pi[9] =  0.048466; pi[10] =  0.086209; pi[11] = 0.0620286;
952  pi[12] = 0.0195027; pi[13] =  0.0384319; pi[14] =  0.0457631; pi[15] = 0.0695179;
953  pi[16] =  0.0610127; pi[17] =  0.0143859; pi[18] =  0.0352742; pi[19] =  0.0708956;
954
955  return 1;
956}
957
958/*********************************************************/
959
960int Init_Qmat_RtREV(double *daa, double *pi)
961{
962    /*
963       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
964       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
965    */
966   
967    /*
968    Dimmic M.W., J.S. Rest, D.P. Mindell, and D. Goldstein. 2002. RArtREV:
969    An amino acid substitution matrix for inference of retrovirus and
970    reverse transcriptase phylogeny. Journal of Molecular Evolution
971    55: 65-73.
972    */
973
974    int i,j,naa;
975    naa = 20;
976
977    daa[1*20+0]= 34;         daa[2*20+0]= 51;         daa[2*20+1]= 35;         daa[3*20+0]= 10;         
978    daa[3*20+1]= 30;         daa[3*20+2]= 384;        daa[4*20+0]= 439;        daa[4*20+1]= 92;         
979    daa[4*20+2]= 128;        daa[4*20+3]= 1;          daa[5*20+0]= 32;         daa[5*20+1]= 221;       
980    daa[5*20+2]= 236;        daa[5*20+3]= 78;         daa[5*20+4]= 70;         daa[6*20+0]= 81;         
981    daa[6*20+1]= 10;         daa[6*20+2]= 79;         daa[6*20+3]= 542;        daa[6*20+4]= 1;         
982    daa[6*20+5]= 372;        daa[7*20+0]= 135;        daa[7*20+1]= 41;         daa[7*20+2]= 94;         
983    daa[7*20+3]= 61;         daa[7*20+4]= 48;         daa[7*20+5]= 18;         daa[7*20+6]= 70;         
984    daa[8*20+0]= 30;         daa[8*20+1]= 90;         daa[8*20+2]= 320;        daa[8*20+3]= 91;         
985    daa[8*20+4]= 124;        daa[8*20+5]= 387;        daa[8*20+6]= 34;         daa[8*20+7]= 68;         
986    daa[9*20+0]= 1;          daa[9*20+1]= 24;         daa[9*20+2]= 35;         daa[9*20+3]= 1;         
987    daa[9*20+4]= 104;        daa[9*20+5]= 33;         daa[9*20+6]= 1;          daa[9*20+7]= 1;         
988    daa[9*20+8]= 34;         daa[10*20+0]= 45;        daa[10*20+1]= 18;        daa[10*20+2]= 15;       
989    daa[10*20+3]= 5;         daa[10*20+4]= 110;       daa[10*20+5]= 54;        daa[10*20+6]= 21;       
990    daa[10*20+7]= 3;         daa[10*20+8]= 51;        daa[10*20+9]= 385;       daa[11*20+0]= 38;       
991    daa[11*20+1]= 593;       daa[11*20+2]= 123;       daa[11*20+3]= 20;        daa[11*20+4]= 16;       
992    daa[11*20+5]= 309;       daa[11*20+6]= 141;       daa[11*20+7]= 30;        daa[11*20+8]= 76;       
993    daa[11*20+9]= 34;        daa[11*20+10]= 23;       daa[12*20+0]= 235;       daa[12*20+1]= 57;       
994    daa[12*20+2]= 1;         daa[12*20+3]= 1;         daa[12*20+4]= 156;       daa[12*20+5]= 158;       
995    daa[12*20+6]= 1;         daa[12*20+7]= 37;        daa[12*20+8]= 116;       daa[12*20+9]= 375;       
996    daa[12*20+10]= 581;      daa[12*20+11]= 134;      daa[13*20+0]= 1;         daa[13*20+1]= 7;         
997    daa[13*20+2]= 49;        daa[13*20+3]= 1;         daa[13*20+4]= 70;        daa[13*20+5]= 1;         
998    daa[13*20+6]= 1;         daa[13*20+7]= 7;         daa[13*20+8]= 141;       daa[13*20+9]= 64;       
999    daa[13*20+10]= 179;      daa[13*20+11]= 14;       daa[13*20+12]= 247;      daa[14*20+0]= 97;       
1000    daa[14*20+1]= 24;        daa[14*20+2]= 33;        daa[14*20+3]= 55;        daa[14*20+4]= 1;         
1001    daa[14*20+5]= 68;        daa[14*20+6]= 52;        daa[14*20+7]= 17;        daa[14*20+8]= 44;       
1002    daa[14*20+9]= 10;        daa[14*20+10]= 22;       daa[14*20+11]= 43;       daa[14*20+12]= 1;       
1003    daa[14*20+13]= 11;       daa[15*20+0]= 460;       daa[15*20+1]= 102;       daa[15*20+2]= 294;       
1004    daa[15*20+3]= 136;       daa[15*20+4]= 75;        daa[15*20+5]= 225;       daa[15*20+6]= 95;       
1005    daa[15*20+7]= 152;       daa[15*20+8]= 183;       daa[15*20+9]= 4;         daa[15*20+10]= 24;       
1006    daa[15*20+11]= 77;       daa[15*20+12]= 1;        daa[15*20+13]= 20;       daa[15*20+14]= 134;     
1007    daa[16*20+0]= 258;       daa[16*20+1]= 64;        daa[16*20+2]= 148;       daa[16*20+3]= 55;       
1008    daa[16*20+4]= 117;       daa[16*20+5]= 146;       daa[16*20+6]= 82;        daa[16*20+7]= 7;         
1009    daa[16*20+8]= 49;        daa[16*20+9]= 72;        daa[16*20+10]= 25;       daa[16*20+11]= 110;     
1010    daa[16*20+12]= 131;      daa[16*20+13]= 69;       daa[16*20+14]= 62;       daa[16*20+15]= 671;     
1011    daa[17*20+0]= 5;         daa[17*20+1]= 13;        daa[17*20+2]= 16;        daa[17*20+3]= 1;         
1012    daa[17*20+4]= 55;        daa[17*20+5]= 10;        daa[17*20+6]= 17;        daa[17*20+7]= 23;       
1013    daa[17*20+8]= 48;        daa[17*20+9]= 39;        daa[17*20+10]= 47;       daa[17*20+11]= 6;       
1014    daa[17*20+12]= 111;      daa[17*20+13]= 182;      daa[17*20+14]= 9;        daa[17*20+15]= 14;       
1015    daa[17*20+16]= 1;        daa[18*20+0]= 55;        daa[18*20+1]= 47;        daa[18*20+2]= 28;       
1016    daa[18*20+3]= 1;         daa[18*20+4]= 131;       daa[18*20+5]= 45;        daa[18*20+6]= 1;         
1017    daa[18*20+7]= 21;        daa[18*20+8]= 307;       daa[18*20+9]= 26;        daa[18*20+10]= 64;       
1018    daa[18*20+11]= 1;        daa[18*20+12]= 74;       daa[18*20+13]= 1017;     daa[18*20+14]= 14;       
1019    daa[18*20+15]= 31;       daa[18*20+16]= 34;       daa[18*20+17]= 176;      daa[19*20+0]= 197;       
1020    daa[19*20+1]= 29;        daa[19*20+2]= 21;        daa[19*20+3]= 6;         daa[19*20+4]= 295;       
1021    daa[19*20+5]= 36;        daa[19*20+6]= 35;        daa[19*20+7]= 3;         daa[19*20+8]= 1;         
1022    daa[19*20+9]= 1048;      daa[19*20+10]= 112;      daa[19*20+11]= 19;       daa[19*20+12]= 236;     
1023    daa[19*20+13]= 92;       daa[19*20+14]= 25;       daa[19*20+15]= 39;       daa[19*20+16]= 196;     
1024    daa[19*20+17]= 26;       daa[19*20+18]= 59;       
1025
1026    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1027
1028    pi[0]= 0.0646;           pi[1]= 0.0453;           pi[2]= 0.0376;           pi[3]= 0.0422;           
1029    pi[4]= 0.0114;           pi[5]= 0.0606;           pi[6]= 0.0607;           pi[7]= 0.0639;           
1030    pi[8]= 0.0273;           pi[9]= 0.0679;           pi[10]= 0.1018;          pi[11]= 0.0751;         
1031    pi[12]= 0.015;           pi[13]= 0.0287;          pi[14]= 0.0681;          pi[15]= 0.0488;         
1032    pi[16]= 0.0622;          pi[17]= 0.0251;          pi[18]= 0.0318;          pi[19]= 0.0619;         
1033    return 1;
1034}
1035
1036/*********************************************************/
1037
1038int Init_Qmat_CpREV(double *daa, double *pi)
1039{
1040    /*
1041       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
1042       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
1043    */
1044
1045    /*
1046      Adachi, J., P. Waddell, W. Martin, and M. Hasegawa. 2000. Plastid         
1047      genome phylogeny and a model of amino acid substitution for proteins   
1048      encoded by chloroplast DNA. Journal of Molecular Evolution             
1049      50:348-358.
1050    */
1051
1052    int i,j,naa;
1053    naa = 20;
1054
1055    daa[1*20+0]= 105;        daa[2*20+0]= 227;        daa[2*20+1]= 357;        daa[3*20+0]= 175;       
1056    daa[3*20+1]= 43;         daa[3*20+2]= 4435;       daa[4*20+0]= 669;        daa[4*20+1]= 823;       
1057    daa[4*20+2]= 538;        daa[4*20+3]= 10;         daa[5*20+0]= 157;        daa[5*20+1]= 1745;       
1058    daa[5*20+2]= 768;        daa[5*20+3]= 400;        daa[5*20+4]= 10;         daa[6*20+0]= 499;       
1059    daa[6*20+1]= 152;        daa[6*20+2]= 1055;       daa[6*20+3]= 3691;       daa[6*20+4]= 10;         
1060    daa[6*20+5]= 3122;       daa[7*20+0]= 665;        daa[7*20+1]= 243;        daa[7*20+2]= 653;       
1061    daa[7*20+3]= 431;        daa[7*20+4]= 303;        daa[7*20+5]= 133;        daa[7*20+6]= 379;       
1062    daa[8*20+0]= 66;         daa[8*20+1]= 715;        daa[8*20+2]= 1405;       daa[8*20+3]= 331;       
1063    daa[8*20+4]= 441;        daa[8*20+5]= 1269;       daa[8*20+6]= 162;        daa[8*20+7]= 19;         
1064    daa[9*20+0]= 145;        daa[9*20+1]= 136;        daa[9*20+2]= 168;        daa[9*20+3]= 10;         
1065    daa[9*20+4]= 280;        daa[9*20+5]= 92;         daa[9*20+6]= 148;        daa[9*20+7]= 40;         
1066    daa[9*20+8]= 29;         daa[10*20+0]= 197;       daa[10*20+1]= 203;       daa[10*20+2]= 113;       
1067    daa[10*20+3]= 10;        daa[10*20+4]= 396;       daa[10*20+5]= 286;       daa[10*20+6]= 82;       
1068    daa[10*20+7]= 20;        daa[10*20+8]= 66;        daa[10*20+9]= 1745;      daa[11*20+0]= 236;       
1069    daa[11*20+1]= 4482;      daa[11*20+2]= 2430;      daa[11*20+3]= 412;       daa[11*20+4]= 48;       
1070    daa[11*20+5]= 3313;      daa[11*20+6]= 2629;      daa[11*20+7]= 263;       daa[11*20+8]= 305;       
1071    daa[11*20+9]= 345;       daa[11*20+10]= 218;      daa[12*20+0]= 185;       daa[12*20+1]= 125;       
1072    daa[12*20+2]= 61;        daa[12*20+3]= 47;        daa[12*20+4]= 159;       daa[12*20+5]= 202;       
1073    daa[12*20+6]= 113;       daa[12*20+7]= 21;        daa[12*20+8]= 10;        daa[12*20+9]= 1772;     
1074    daa[12*20+10]= 1351;     daa[12*20+11]= 193;      daa[13*20+0]= 68;        daa[13*20+1]= 53;       
1075    daa[13*20+2]= 97;        daa[13*20+3]= 22;        daa[13*20+4]= 726;       daa[13*20+5]= 10;       
1076    daa[13*20+6]= 145;       daa[13*20+7]= 25;        daa[13*20+8]= 127;       daa[13*20+9]= 454;       
1077    daa[13*20+10]= 1268;     daa[13*20+11]= 72;       daa[13*20+12]= 327;      daa[14*20+0]= 490;       
1078    daa[14*20+1]= 87;        daa[14*20+2]= 173;       daa[14*20+3]= 170;       daa[14*20+4]= 285;       
1079    daa[14*20+5]= 323;       daa[14*20+6]= 185;       daa[14*20+7]= 28;        daa[14*20+8]= 152;       
1080    daa[14*20+9]= 117;       daa[14*20+10]= 219;      daa[14*20+11]= 302;      daa[14*20+12]= 100;     
1081    daa[14*20+13]= 43;       daa[15*20+0]= 2440;      daa[15*20+1]= 385;       daa[15*20+2]= 2085;     
1082    daa[15*20+3]= 590;       daa[15*20+4]= 2331;      daa[15*20+5]= 396;       daa[15*20+6]= 568;       
1083    daa[15*20+7]= 691;       daa[15*20+8]= 303;       daa[15*20+9]= 216;       daa[15*20+10]= 516;     
1084    daa[15*20+11]= 868;      daa[15*20+12]= 93;       daa[15*20+13]= 487;      daa[15*20+14]= 1202;     
1085    daa[16*20+0]= 1340;      daa[16*20+1]= 314;       daa[16*20+2]= 1393;      daa[16*20+3]= 266;       
1086    daa[16*20+4]= 576;       daa[16*20+5]= 241;       daa[16*20+6]= 369;       daa[16*20+7]= 92;       
1087    daa[16*20+8]= 32;        daa[16*20+9]= 1040;      daa[16*20+10]= 156;      daa[16*20+11]= 918;     
1088    daa[16*20+12]= 645;      daa[16*20+13]= 148;      daa[16*20+14]= 260;      daa[16*20+15]= 2151;     
1089    daa[17*20+0]= 14;        daa[17*20+1]= 230;       daa[17*20+2]= 40;        daa[17*20+3]= 18;       
1090    daa[17*20+4]= 435;       daa[17*20+5]= 53;        daa[17*20+6]= 63;        daa[17*20+7]= 82;       
1091    daa[17*20+8]= 69;        daa[17*20+9]= 42;        daa[17*20+10]= 159;      daa[17*20+11]= 10;       
1092    daa[17*20+12]= 86;       daa[17*20+13]= 468;      daa[17*20+14]= 49;       daa[17*20+15]= 73;       
1093    daa[17*20+16]= 29;       daa[18*20+0]= 56;        daa[18*20+1]= 323;       daa[18*20+2]= 754;       
1094    daa[18*20+3]= 281;       daa[18*20+4]= 1466;      daa[18*20+5]= 391;       daa[18*20+6]= 142;       
1095    daa[18*20+7]= 10;        daa[18*20+8]= 1971;      daa[18*20+9]= 89;        daa[18*20+10]= 189;     
1096    daa[18*20+11]= 247;      daa[18*20+12]= 215;      daa[18*20+13]= 2370;     daa[18*20+14]= 97;       
1097    daa[18*20+15]= 522;      daa[18*20+16]= 71;       daa[18*20+17]= 346;      daa[19*20+0]= 968;       
1098    daa[19*20+1]= 92;        daa[19*20+2]= 83;        daa[19*20+3]= 75;        daa[19*20+4]= 592;       
1099    daa[19*20+5]= 54;        daa[19*20+6]= 200;       daa[19*20+7]= 91;        daa[19*20+8]= 25;       
1100    daa[19*20+9]= 4797;      daa[19*20+10]= 865;      daa[19*20+11]= 249;      daa[19*20+12]= 475;     
1101    daa[19*20+13]= 317;      daa[19*20+14]= 122;      daa[19*20+15]= 167;      daa[19*20+16]= 760;     
1102    daa[19*20+17]= 10;       daa[19*20+18]= 119;     
1103
1104    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1105
1106    pi[0]= 0.076;            pi[1]= 0.062;            pi[2]= 0.041;            pi[3]= 0.037;           
1107    pi[4]= 0.009;            pi[5]= 0.038;            pi[6]= 0.049;            pi[7]= 0.084;           
1108    pi[8]= 0.025;            pi[9]= 0.081;            pi[10]= 0.101;           pi[11]= 0.05;           
1109    pi[12]= 0.022;           pi[13]= 0.051;           pi[14]= 0.043;           pi[15]= 0.062;           
1110    pi[16]= 0.054;           pi[17]= 0.018;           pi[18]= 0.031;           pi[19]= 0.066;           
1111    return 1;
1112}
1113
1114/*********************************************************/
1115
1116int Init_Qmat_VT(double *daa, double *pi)
1117{
1118    /*
1119       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
1120       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
1121    */
1122
1123    /*
1124      Muller, T., and M. Vingron. 2000. Modeling amino acid replacement.         
1125      Journal of Computational Biology 7:761-776.                             
1126    */
1127
1128    int i,j,naa;
1129    naa = 20;
1130
1131    daa[1*20+0]= 0.233108;   daa[2*20+0]= 0.199097;   daa[2*20+1]= 0.210797;   daa[3*20+0]= 0.265145;   
1132    daa[3*20+1]= 0.105191;   daa[3*20+2]= 0.883422;   daa[4*20+0]= 0.227333;   daa[4*20+1]= 0.031726;   
1133    daa[4*20+2]= 0.027495;   daa[4*20+3]= 0.010313;   daa[5*20+0]= 0.310084;   daa[5*20+1]= 0.493763;   
1134    daa[5*20+2]= 0.2757;     daa[5*20+3]= 0.205842;   daa[5*20+4]= 0.004315;   daa[6*20+0]= 0.567957;   
1135    daa[6*20+1]= 0.25524;    daa[6*20+2]= 0.270417;   daa[6*20+3]= 1.599461;   daa[6*20+4]= 0.005321;   
1136    daa[6*20+5]= 0.960976;   daa[7*20+0]= 0.876213;   daa[7*20+1]= 0.156945;   daa[7*20+2]= 0.362028;   
1137    daa[7*20+3]= 0.311718;   daa[7*20+4]= 0.050876;   daa[7*20+5]= 0.12866;    daa[7*20+6]= 0.250447;   
1138    daa[8*20+0]= 0.078692;   daa[8*20+1]= 0.213164;   daa[8*20+2]= 0.290006;   daa[8*20+3]= 0.134252;   
1139    daa[8*20+4]= 0.016695;   daa[8*20+5]= 0.315521;   daa[8*20+6]= 0.104458;   daa[8*20+7]= 0.058131;   
1140    daa[9*20+0]= 0.222972;   daa[9*20+1]= 0.08151;    daa[9*20+2]= 0.087225;   daa[9*20+3]= 0.01172;   
1141    daa[9*20+4]= 0.046398;   daa[9*20+5]= 0.054602;   daa[9*20+6]= 0.046589;   daa[9*20+7]= 0.051089;   
1142    daa[9*20+8]= 0.020039;   daa[10*20+0]= 0.42463;   daa[10*20+1]= 0.192364;  daa[10*20+2]= 0.069245; 
1143    daa[10*20+3]= 0.060863;  daa[10*20+4]= 0.091709;  daa[10*20+5]= 0.24353;   daa[10*20+6]= 0.151924; 
1144    daa[10*20+7]= 0.087056;  daa[10*20+8]= 0.103552;  daa[10*20+9]= 2.08989;   daa[11*20+0]= 0.393245; 
1145    daa[11*20+1]= 1.755838;  daa[11*20+2]= 0.50306;   daa[11*20+3]= 0.261101;  daa[11*20+4]= 0.004067; 
1146    daa[11*20+5]= 0.738208;  daa[11*20+6]= 0.88863;   daa[11*20+7]= 0.193243;  daa[11*20+8]= 0.153323; 
1147    daa[11*20+9]= 0.093181;  daa[11*20+10]= 0.201204; daa[12*20+0]= 0.21155;   daa[12*20+1]= 0.08793;   
1148    daa[12*20+2]= 0.05742;   daa[12*20+3]= 0.012182;  daa[12*20+4]= 0.02369;   daa[12*20+5]= 0.120801; 
1149    daa[12*20+6]= 0.058643;  daa[12*20+7]= 0.04656;   daa[12*20+8]= 0.021157;  daa[12*20+9]= 0.493845; 
1150    daa[12*20+10]= 1.105667; daa[12*20+11]= 0.096474; daa[13*20+0]= 0.116646;  daa[13*20+1]= 0.042569; 
1151    daa[13*20+2]= 0.039769;  daa[13*20+3]= 0.016577;  daa[13*20+4]= 0.051127;  daa[13*20+5]= 0.026235; 
1152    daa[13*20+6]= 0.028168;  daa[13*20+7]= 0.050143;  daa[13*20+8]= 0.079807;  daa[13*20+9]= 0.32102;   
1153    daa[13*20+10]= 0.946499; daa[13*20+11]= 0.038261; daa[13*20+12]= 0.173052; daa[14*20+0]= 0.399143; 
1154    daa[14*20+1]= 0.12848;   daa[14*20+2]= 0.083956;  daa[14*20+3]= 0.160063;  daa[14*20+4]= 0.011137; 
1155    daa[14*20+5]= 0.15657;   daa[14*20+6]= 0.205134;  daa[14*20+7]= 0.124492;  daa[14*20+8]= 0.078892; 
1156    daa[14*20+9]= 0.054797;  daa[14*20+10]= 0.169784; daa[14*20+11]= 0.212302; daa[14*20+12]= 0.010363; 
1157    daa[14*20+13]= 0.042564; daa[15*20+0]= 1.817198;  daa[15*20+1]= 0.292327;  daa[15*20+2]= 0.847049; 
1158    daa[15*20+3]= 0.461519;  daa[15*20+4]= 0.17527;   daa[15*20+5]= 0.358017;  daa[15*20+6]= 0.406035; 
1159    daa[15*20+7]= 0.612843;  daa[15*20+8]= 0.167406;  daa[15*20+9]= 0.081567;  daa[15*20+10]= 0.214977; 
1160    daa[15*20+11]= 0.400072; daa[15*20+12]= 0.090515; daa[15*20+13]= 0.138119; daa[15*20+14]= 0.430431; 
1161    daa[16*20+0]= 0.877877;  daa[16*20+1]= 0.204109;  daa[16*20+2]= 0.471268;  daa[16*20+3]= 0.178197; 
1162    daa[16*20+4]= 0.079511;  daa[16*20+5]= 0.248992;  daa[16*20+6]= 0.321028;  daa[16*20+7]= 0.136266; 
1163    daa[16*20+8]= 0.101117;  daa[16*20+9]= 0.376588;  daa[16*20+10]= 0.243227; daa[16*20+11]= 0.446646; 
1164    daa[16*20+12]= 0.184609; daa[16*20+13]= 0.08587;  daa[16*20+14]= 0.207143; daa[16*20+15]= 1.767766; 
1165    daa[17*20+0]= 0.030309;  daa[17*20+1]= 0.046417;  daa[17*20+2]= 0.010459;  daa[17*20+3]= 0.011393; 
1166    daa[17*20+4]= 0.007732;  daa[17*20+5]= 0.021248;  daa[17*20+6]= 0.018844;  daa[17*20+7]= 0.02399;   
1167    daa[17*20+8]= 0.020009;  daa[17*20+9]= 0.034954;  daa[17*20+10]= 0.083439; daa[17*20+11]= 0.023321; 
1168    daa[17*20+12]= 0.022019; daa[17*20+13]= 0.12805;  daa[17*20+14]= 0.014584; daa[17*20+15]= 0.035933; 
1169    daa[17*20+16]= 0.020437; daa[18*20+0]= 0.087061;  daa[18*20+1]= 0.09701;   daa[18*20+2]= 0.093268; 
1170    daa[18*20+3]= 0.051664;  daa[18*20+4]= 0.042823;  daa[18*20+5]= 0.062544;  daa[18*20+6]= 0.0552;   
1171    daa[18*20+7]= 0.037568;  daa[18*20+8]= 0.286027;  daa[18*20+9]= 0.086237;  daa[18*20+10]= 0.189842; 
1172    daa[18*20+11]= 0.068689; daa[18*20+12]= 0.073223; daa[18*20+13]= 0.898663; daa[18*20+14]= 0.032043; 
1173    daa[18*20+15]= 0.121979; daa[18*20+16]= 0.094617; daa[18*20+17]= 0.124746; daa[19*20+0]= 1.230985; 
1174    daa[19*20+1]= 0.113146;  daa[19*20+2]= 0.049824;  daa[19*20+3]= 0.048769;  daa[19*20+4]= 0.163831; 
1175    daa[19*20+5]= 0.112027;  daa[19*20+6]= 0.205868;  daa[19*20+7]= 0.082579;  daa[19*20+8]= 0.068575; 
1176    daa[19*20+9]= 3.65443;   daa[19*20+10]= 1.337571; daa[19*20+11]= 0.144587; daa[19*20+12]= 0.307309; 
1177    daa[19*20+13]= 0.247329; daa[19*20+14]= 0.129315; daa[19*20+15]= 0.1277;   daa[19*20+16]= 0.740372; 
1178    daa[19*20+17]= 0.022134; daa[19*20+18]= 0.125733; 
1179
1180    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1181
1182    pi[0]= 0.078837;         pi[1]= 0.051238;         pi[2]= 0.042313;         pi[3]= 0.053066;         
1183    pi[4]= 0.015175;         pi[5]= 0.036713;         pi[6]= 0.061924;         pi[7]= 0.070852;         
1184    pi[8]= 0.023082;         pi[9]= 0.062056;         pi[10]= 0.096371;        pi[11]= 0.057324;       
1185    pi[12]= 0.023771;        pi[13]= 0.043296;        pi[14]= 0.043911;        pi[15]= 0.063403;       
1186    pi[16]= 0.055897;        pi[17]= 0.013272;        pi[18]= 0.034399;        pi[19]= 0.073101;       
1187    return 1;
1188}
1189
1190/*********************************************************/
1191
1192int Init_Qmat_Blosum62 (double *daa, double *pi)
1193{
1194
1195    /*
1196       This model has been 'translated' from John Huelsenbeck and Fredrik Ronquist
1197       MrBayes program into PHYML format by Federico Abascal. Many thanks to them.
1198    */
1199
1200    /*
1201      Henikoff, S., and J. G. Henikoff. 1992. Amino acid substitution           
1202      matrices from protein blocks. Proc. Natl. Acad. Sci., U.S.A.           
1203      89:10915-10919.                                                         
1204    */
1205
1206    int i,j,naa;
1207    naa = 20;
1208
1209    daa[1*20+0]= 0.735790389698;  daa[2*20+0]= 0.485391055466;  daa[2*20+1]= 1.297446705134;  daa[3*20+0]= 0.543161820899; 
1210    daa[3*20+1]= 0.500964408555;  daa[3*20+2]= 3.180100048216;  daa[4*20+0]= 1.45999531047;   daa[4*20+1]= 0.227826574209; 
1211    daa[4*20+2]= 0.397358949897;  daa[4*20+3]= 0.240836614802;  daa[5*20+0]= 1.199705704602;  daa[5*20+1]= 3.020833610064; 
1212    daa[5*20+2]= 1.839216146992;  daa[5*20+3]= 1.190945703396;  daa[5*20+4]= 0.32980150463;   daa[6*20+0]= 1.1709490428;   
1213    daa[6*20+1]= 1.36057419042;   daa[6*20+2]= 1.24048850864;   daa[6*20+3]= 3.761625208368;  daa[6*20+4]= 0.140748891814; 
1214    daa[6*20+5]= 5.528919177928;  daa[7*20+0]= 1.95588357496;   daa[7*20+1]= 0.418763308518;  daa[7*20+2]= 1.355872344485; 
1215    daa[7*20+3]= 0.798473248968;  daa[7*20+4]= 0.418203192284;  daa[7*20+5]= 0.609846305383;  daa[7*20+6]= 0.423579992176; 
1216    daa[8*20+0]= 0.716241444998;  daa[8*20+1]= 1.456141166336;  daa[8*20+2]= 2.414501434208;  daa[8*20+3]= 0.778142664022; 
1217    daa[8*20+4]= 0.354058109831;  daa[8*20+5]= 2.43534113114;   daa[8*20+6]= 1.626891056982;  daa[8*20+7]= 0.539859124954; 
1218    daa[9*20+0]= 0.605899003687;  daa[9*20+1]= 0.232036445142;  daa[9*20+2]= 0.283017326278;  daa[9*20+3]= 0.418555732462; 
1219    daa[9*20+4]= 0.774894022794;  daa[9*20+5]= 0.236202451204;  daa[9*20+6]= 0.186848046932;  daa[9*20+7]= 0.189296292376; 
1220    daa[9*20+8]= 0.252718447885;  daa[10*20+0]= 0.800016530518; daa[10*20+1]= 0.622711669692; daa[10*20+2]= 0.211888159615; 
1221    daa[10*20+3]= 0.218131577594; daa[10*20+4]= 0.831842640142; daa[10*20+5]= 0.580737093181; daa[10*20+6]= 0.372625175087; 
1222    daa[10*20+7]= 0.217721159236; daa[10*20+8]= 0.348072209797; daa[10*20+9]= 3.890963773304; daa[11*20+0]= 1.295201266783; 
1223    daa[11*20+1]= 5.411115141489; daa[11*20+2]= 1.593137043457; daa[11*20+3]= 1.032447924952; daa[11*20+4]= 0.285078800906; 
1224    daa[11*20+5]= 3.945277674515; daa[11*20+6]= 2.802427151679; daa[11*20+7]= 0.752042440303; daa[11*20+8]= 1.022507035889; 
1225    daa[11*20+9]= 0.406193586642; daa[11*20+10]= 0.445570274261;daa[12*20+0]= 1.253758266664; daa[12*20+1]= 0.983692987457; 
1226    daa[12*20+2]= 0.648441278787; daa[12*20+3]= 0.222621897958; daa[12*20+4]= 0.76768882348;  daa[12*20+5]= 2.494896077113; 
1227    daa[12*20+6]= 0.55541539747;  daa[12*20+7]= 0.459436173579; daa[12*20+8]= 0.984311525359; daa[12*20+9]= 3.364797763104; 
1228    daa[12*20+10]= 6.030559379572;daa[12*20+11]= 1.073061184332;daa[13*20+0]= 0.492964679748; daa[13*20+1]= 0.371644693209; 
1229    daa[13*20+2]= 0.354861249223; daa[13*20+3]= 0.281730694207; daa[13*20+4]= 0.441337471187; daa[13*20+5]= 0.14435695975; 
1230    daa[13*20+6]= 0.291409084165; daa[13*20+7]= 0.368166464453; daa[13*20+8]= 0.714533703928; daa[13*20+9]= 1.517359325954; 
1231    daa[13*20+10]= 2.064839703237;daa[13*20+11]= 0.266924750511;daa[13*20+12]= 1.77385516883; daa[14*20+0]= 1.173275900924; 
1232    daa[14*20+1]= 0.448133661718; daa[14*20+2]= 0.494887043702; daa[14*20+3]= 0.730628272998; daa[14*20+4]= 0.356008498769; 
1233    daa[14*20+5]= 0.858570575674; daa[14*20+6]= 0.926563934846; daa[14*20+7]= 0.504086599527; daa[14*20+8]= 0.527007339151; 
1234    daa[14*20+9]= 0.388355409206; daa[14*20+10]= 0.374555687471;daa[14*20+11]= 1.047383450722;daa[14*20+12]= 0.454123625103;
1235    daa[14*20+13]= 0.233597909629;daa[15*20+0]= 4.325092687057; daa[15*20+1]= 1.12278310421;  daa[15*20+2]= 2.904101656456; 
1236    daa[15*20+3]= 1.582754142065; daa[15*20+4]= 1.197188415094; daa[15*20+5]= 1.934870924596; daa[15*20+6]= 1.769893238937; 
1237    daa[15*20+7]= 1.509326253224; daa[15*20+8]= 1.11702976291;  daa[15*20+9]= 0.35754441246;  daa[15*20+10]= 0.352969184527;
1238    daa[15*20+11]= 1.752165917819;daa[15*20+12]= 0.918723415746;daa[15*20+13]= 0.540027644824;daa[15*20+14]= 1.169129577716;
1239    daa[16*20+0]= 1.729178019485; daa[16*20+1]= 0.914665954563; daa[16*20+2]= 1.898173634533; daa[16*20+3]= 0.934187509431; 
1240    daa[16*20+4]= 1.119831358516; daa[16*20+5]= 1.277480294596; daa[16*20+6]= 1.071097236007; daa[16*20+7]= 0.641436011405; 
1241    daa[16*20+8]= 0.585407090225; daa[16*20+9]= 1.17909119726;  daa[16*20+10]= 0.915259857694;daa[16*20+11]= 1.303875200799;
1242    daa[16*20+12]= 1.488548053722;daa[16*20+13]= 0.488206118793;daa[16*20+14]= 1.005451683149;daa[16*20+15]= 5.15155629227; 
1243    daa[17*20+0]= 0.465839367725; daa[17*20+1]= 0.426382310122; daa[17*20+2]= 0.191482046247; daa[17*20+3]= 0.145345046279; 
1244    daa[17*20+4]= 0.527664418872; daa[17*20+5]= 0.758653808642; daa[17*20+6]= 0.407635648938; daa[17*20+7]= 0.508358924638; 
1245    daa[17*20+8]= 0.30124860078;  daa[17*20+9]= 0.34198578754;  daa[17*20+10]= 0.6914746346;  daa[17*20+11]= 0.332243040634;
1246    daa[17*20+12]= 0.888101098152;daa[17*20+13]= 2.074324893497;daa[17*20+14]= 0.252214830027;daa[17*20+15]= 0.387925622098;
1247    daa[17*20+16]= 0.513128126891;daa[18*20+0]= 0.718206697586; daa[18*20+1]= 0.720517441216; daa[18*20+2]= 0.538222519037; 
1248    daa[18*20+3]= 0.261422208965; daa[18*20+4]= 0.470237733696; daa[18*20+5]= 0.95898974285;  daa[18*20+6]= 0.596719300346; 
1249    daa[18*20+7]= 0.308055737035; daa[18*20+8]= 4.218953969389; daa[18*20+9]= 0.674617093228; daa[18*20+10]= 0.811245856323;
1250    daa[18*20+11]= 0.7179934869;  daa[18*20+12]= 0.951682162246;daa[18*20+13]= 6.747260430801;daa[18*20+14]= 0.369405319355;
1251    daa[18*20+15]= 0.796751520761;daa[18*20+16]= 0.801010243199;daa[18*20+17]= 4.054419006558;daa[19*20+0]= 2.187774522005; 
1252    daa[19*20+1]= 0.438388343772; daa[19*20+2]= 0.312858797993; daa[19*20+3]= 0.258129289418; daa[19*20+4]= 1.116352478606; 
1253    daa[19*20+5]= 0.530785790125; daa[19*20+6]= 0.524253846338; daa[19*20+7]= 0.25334079019;  daa[19*20+8]= 0.20155597175; 
1254    daa[19*20+9]= 8.311839405458; daa[19*20+10]= 2.231405688913;daa[19*20+11]= 0.498138475304;daa[19*20+12]= 2.575850755315;
1255    daa[19*20+13]= 0.838119610178;daa[19*20+14]= 0.496908410676;daa[19*20+15]= 0.561925457442;daa[19*20+16]= 2.253074051176;
1256    daa[19*20+17]= 0.266508731426;daa[19*20+18]= 1;             
1257
1258    for (i=0; i<naa; i++)  for (j=0; j<i; j++)  daa[j*naa+i] = daa[i*naa+j];
1259
1260    pi[0]= 0.074;                 pi[1]= 0.052;                 pi[2]= 0.045;                 pi[3]= 0.054;                 
1261    pi[4]= 0.025;                 pi[5]= 0.034;                 pi[6]= 0.054;                 pi[7]= 0.074;                 
1262    pi[8]= 0.026;                 pi[9]= 0.068;                 pi[10]= 0.099;                pi[11]= 0.058;               
1263    pi[12]= 0.025;                pi[13]= 0.047;                pi[14]= 0.039;                pi[15]= 0.057;               
1264    pi[16]= 0.051;                pi[17]= 0.013;                pi[18]= 0.032;                pi[19]= 0.073;               
1265
1266    return 1;
1267  }
1268
1269/*********************************************************/
1270
1271int Init_Qmat_MtMam(double *daa, double *pi)
1272{
1273    /*
1274       This model has been 'translated' from Ziheng Yang's PAML program
1275       into PHYML format by Federico Abascal. Many thanks to them.
1276    */
1277
1278    /*
1279      Cao, Y. et al. 1998 Conflict amongst individual mitochondrial
1280      proteins in resolving the phylogeny of eutherian orders. Journal
1281      of Molecular Evolution 15:1600-1611.
1282    */
1283
1284    int i,j,naa;
1285    naa = 20;
1286   
1287    daa[1*20+0]= 32;              daa[2*20+0]= 2;    daa[2*20+1]= 4;               daa[3*20+0]= 11;
1288    daa[3*20+1]= 0;               daa[3*20+2]= 864;  daa[4*20+0]= 0;               daa[4*20+1]= 186;
1289    daa[4*20+2]= 0;               daa[4*20+3]= 0;    daa[5*20+0]= 0;               daa[5*20+1]= 246;
1290    daa[5*20+2]= 8;               daa[5*20+3]= 49;   daa[5*20+4]= 0;               daa[6*20+0]= 0;
1291    daa[6*20+1]= 0;               daa[6*20+2]= 0;    daa[6*20+3]= 569;             daa[6*20+4]= 0;
1292    daa[6*20+5]= 274;             daa[7*20+0]= 78;   daa[7*20+1]= 18;              daa[7*20+2]= 47;
1293    daa[7*20+3]= 79;              daa[7*20+4]= 0;    daa[7*20+5]= 0;               daa[7*20+6]= 22;
1294    daa[8*20+0]= 8;               daa[8*20+1]= 232;  daa[8*20+2]= 458;             daa[8*20+3]= 11;
1295    daa[8*20+4]= 305;             daa[8*20+5]= 550;  daa[8*20+6]= 22;              daa[8*20+7]= 0;
1296    daa[9*20+0]= 75;              daa[9*20+1]= 0;    daa[9*20+2]= 19;              daa[9*20+3]= 0;
1297    daa[9*20+4]= 41;              daa[9*20+5]= 0;    daa[9*20+6]= 0;               daa[9*20+7]= 0;
1298    daa[9*20+8]= 0;               daa[10*20+0]= 21;  daa[10*20+1]= 6;              daa[10*20+2]= 0;
1299    daa[10*20+3]= 0;              daa[10*20+4]= 27;  daa[10*20+5]= 20;             daa[10*20+6]= 0;
1300    daa[10*20+7]= 0;              daa[10*20+8]= 26;  daa[10*20+9]= 232;            daa[11*20+0]= 0;
1301    daa[11*20+1]= 50;             daa[11*20+2]= 408; daa[11*20+3]= 0;              daa[11*20+4]= 0;
1302    daa[11*20+5]= 242;            daa[11*20+6]= 215; daa[11*20+7]= 0;              daa[11*20+8]= 0;
1303    daa[11*20+9]= 6;              daa[11*20+10]= 4;  daa[12*20+0]= 76;             daa[12*20+1]= 0;
1304    daa[12*20+2]= 21;             daa[12*20+3]= 0;   daa[12*20+4]= 0;              daa[12*20+5]= 22;
1305    daa[12*20+6]= 0;              daa[12*20+7]= 0;   daa[12*20+8]= 0;              daa[12*20+9]= 378;
1306    daa[12*20+10]= 609;           daa[12*20+11]= 59; daa[13*20+0]= 0;              daa[13*20+1]= 0;
1307    daa[13*20+2]= 6;              daa[13*20+3]= 5;   daa[13*20+4]= 7;              daa[13*20+5]= 0;
1308    daa[13*20+6]= 0;              daa[13*20+7]= 0;   daa[13*20+8]= 0;              daa[13*20+9]= 57;
1309    daa[13*20+10]= 246;           daa[13*20+11]= 0;  daa[13*20+12]= 11;            daa[14*20+0]= 53;
1310    daa[14*20+1]= 9;              daa[14*20+2]= 33;  daa[14*20+3]= 2;              daa[14*20+4]= 0;
1311    daa[14*20+5]= 51;             daa[14*20+6]= 0;   daa[14*20+7]= 0;              daa[14*20+8]= 53;
1312    daa[14*20+9]= 5;              daa[14*20+10]= 43; daa[14*20+11]= 18;            daa[14*20+12]= 0;
1313    daa[14*20+13]= 17;            daa[15*20+0]= 342; daa[15*20+1]= 3;              daa[15*20+2]= 446;
1314    daa[15*20+3]= 16;             daa[15*20+4]= 347; daa[15*20+5]= 30;             daa[15*20+6]= 21;
1315    daa[15*20+7]= 112;            daa[15*20+8]= 20;  daa[15*20+9]= 0;              daa[15*20+10]= 74;
1316    daa[15*20+11]= 65;            daa[15*20+12]= 47; daa[15*20+13]= 90;            daa[15*20+14]= 202;
1317    daa[16*20+0]= 681;            daa[16*20+1]= 0;   daa[16*20+2]= 110;            daa[16*20+3]= 0;
1318    daa[16*20+4]= 114;            daa[16*20+5]= 0;   daa[16*20+6]= 4;              daa[16*20+7]= 0;
1319    daa[16*20+8]= 1;              daa[16*20+9]= 360; daa[16*20+10]= 34;            daa[16*20+11]= 50;
1320    daa[16*20+12]= 691;           daa[16*20+13]= 8;  daa[16*20+14]= 78;            daa[16*20+15]= 614;
1321    daa[17*20+0]= 5;              daa[17*20+1]= 16;  daa[17*20+2]= 6;              daa[17*20+3]= 0;
1322    daa[17*20+4]= 65;             daa[17*20+5]= 0;   daa[17*20+6]= 0;              daa[17*20+7]= 0;
1323    daa[17*20+8]= 0;              daa[17*20+9]= 0;   daa[17*20+10]= 12;            daa[17*20+11]= 0;
1324    daa[17*20+12]= 13;            daa[17*20+13]= 0;  daa[17*20+14]= 7;             daa[17*20+15]= 17;
1325    daa[17*20+16]= 0;             daa[18*20+0]= 0;   daa[18*20+1]= 0;              daa[18*20+2]= 156;
1326    daa[18*20+3]= 0;              daa[18*20+4]= 530; daa[18*20+5]= 54;             daa[18*20+6]= 0;
1327    daa[18*20+7]= 1;              daa[18*20+8]= 1525;daa[18*20+9]= 16;             daa[18*20+10]= 25;
1328    daa[18*20+11]= 67;            daa[18*20+12]= 0;  daa[18*20+13]= 682;           daa[18*20+14]= 8;
1329    daa[18*20+15]= 107;           daa[18*20+16]= 0;  daa[18*20+17]= 14;            daa[19*20+0]= 398;
1330    daa[19*20+1]= 0;              daa[19*20+2]= 0;   daa[19*20+3]= 10;             daa[19*20+4]= 0;
1331    daa[19*20+5]= 33;             daa[19*20+6]= 20;  daa[19*20+7]= 5;              daa[19*20+8]= 0;
1332    daa[19*20+9]= 2220;           daa[19*20+10]= 100;daa[19*20+11]= 0;             daa[19*20+12]= 832;
1333    daa[19*20+13]= 6;             daa[19*20+14]= 0;  daa[19*20+15]= 0;             daa[19*20+16]= 237;
1334    daa[19*20+17]= 0;             daa[19*20+18]= 0;
1335   
1336    for (i=0; i<naa; i++) for (j=0; j<i; j++) daa[j*naa+i] = daa[i*naa+j];
1337   
1338    pi[0]= 0.0692;  pi[1]=  0.0184;  pi[2]= 0.04;    pi[3]= 0.0186;
1339    pi[4]= 0.0065;  pi[5]=  0.0238;  pi[6]= 0.0236;  pi[7]= 0.0557;
1340    pi[8]= 0.0277;  pi[9]=  0.0905;  pi[10]=0.1675;  pi[11]= 0.0221;
1341    pi[12]=0.0561;  pi[13]= 0.0611;  pi[14]=0.0536;  pi[15]= 0.0725;
1342    pi[16]=0.087;   pi[17]= 0.0293;  pi[18]=0.034;   pi[19]= 0.0428;
1343   
1344    return 1;
1345}
1346
1347                         
1348/*********************************************************/
1349
1350void Init_Model(allseq *data, model *mod)
1351{
1352  int i,j;
1353  int ns;
1354  double sum,aux;
1355  int result;
1356  double *dr, *di, *space;
1357
1358
1359  if(data) mod->seq_len = data->init_len;
1360
1361   
1362  For(i,mod->n_catg) 
1363    {
1364      mod->rr[i]      = 1.0;
1365      mod->r_proba[i] = 1.0;
1366    }
1367
1368 
1369  if(!mod->invar) For(i,data->crunch_len) data->invar[i] = 0;
1370
1371  ns = mod->ns;
1372  mod->datatype = (mod->whichmodel>=10)?((mod->whichmodel>20)?(0):(1)):(0);
1373
1374
1375  dr      = (double *)mCalloc(  ns,sizeof(double));
1376  di      = (double *)mCalloc(  ns,sizeof(double));
1377  space   = (double *)mCalloc(2*ns,sizeof(double));
1378
1379 
1380
1381  For(i,ns)
1382      mod->pi[i] = data->b_frq[i];
1383
1384  if(mod->datatype == NT) /* Nucleotides */
1385    { 
1386      if(mod->whichmodel < 40)
1387        {
1388          /* init for nucleotides */
1389          mod->lambda = 1.;
1390         
1391          if((mod->whichmodel==1) || (mod->whichmodel==2))
1392            {
1393              mod->pi[0] = mod->pi[1] = mod->pi[2] = mod->pi[3] = .25;
1394            }
1395         
1396          if((mod->whichmodel==1) || (mod->whichmodel==3) || (mod->whichmodel==7) || (mod->whichmodel==8))
1397            {
1398              mod->kappa = 1.;
1399            }
1400         
1401          if(mod->whichmodel == 5)
1402            {
1403              aux = ((mod->pi[0]+mod->pi[2])-(mod->pi[1]+mod->pi[3]))/(2.*mod->kappa);
1404              mod->lambda = ((mod->pi[1]+mod->pi[3]) + aux)/((mod->pi[0]+mod->pi[2]) - aux); 
1405            }
1406         
1407          if(mod->whichmodel == 7)
1408            {             
1409              mod->custom_mod_string[0] = '0';
1410              mod->custom_mod_string[1] = '1';
1411              mod->custom_mod_string[2] = '2';
1412              mod->custom_mod_string[3] = '3';
1413              mod->custom_mod_string[4] = '4';
1414              mod->custom_mod_string[5] = '5';
1415              Translate_Custom_Mod_String(mod);
1416              Update_Qmat_GTR(mod);
1417            }
1418
1419          if(mod->whichmodel == 8)
1420              {
1421                 
1422              if(mod->user_b_freq[0] > -1.)
1423                  {
1424                      For(i,4) 
1425                          {
1426                              mod->pi[i] = mod->user_b_freq[i];
1427                          }
1428                  }
1429              Update_Qmat_GTR(mod);
1430            }
1431         
1432        }
1433      else
1434        {
1435          /* init for codon model */
1436          For(i,64) mod->pi[i] = 1./61;
1437        }
1438    }
1439  else 
1440    { /* init for amino-acids */
1441      /* see comments of PMat_Empirical for details */
1442      /* read pi and Q from file */
1443     
1444     
1445      switch(mod->whichmodel)
1446        {
1447        case 11 : 
1448          {
1449            Init_Qmat_Dayhoff(mod->mat_Q,mod->pi);
1450            break;
1451          }
1452        case 12 : 
1453          {
1454            Init_Qmat_JTT(mod->mat_Q,mod->pi);
1455            break;
1456          }
1457        case 13 : 
1458          {
1459            Init_Qmat_MtREV(mod->mat_Q,mod->pi);
1460            break;
1461          }
1462        case 14 : 
1463          {
1464            Init_Qmat_WAG(mod->mat_Q,mod->pi);
1465            break;
1466          }
1467        case 15 : 
1468          {
1469            Init_Qmat_DCMut(mod->mat_Q,mod->pi);
1470            break;
1471          }
1472        case 16 : 
1473          {
1474            Init_Qmat_RtREV(mod->mat_Q,mod->pi);
1475            break;
1476          }
1477        case 17 : 
1478          {
1479            Init_Qmat_CpREV(mod->mat_Q,mod->pi);
1480            break;
1481          }
1482        case 18 : 
1483          {
1484            Init_Qmat_VT(mod->mat_Q,mod->pi);
1485            break;
1486          }
1487        case 19 : 
1488          {
1489            Init_Qmat_Blosum62(mod->mat_Q,mod->pi);
1490            break;
1491          }
1492        case 20 : 
1493          {
1494            Init_Qmat_MtMam(mod->mat_Q,mod->pi);
1495            break;
1496          }
1497        default : break;
1498        }
1499
1500
1501      /* multiply the nth col of Q by the nth term of pi/100 just as in PAML */
1502      For(i,ns) For(j,ns) mod->mat_Q[i*ns+j] *= mod->pi[j] / 100.0;
1503
1504     
1505      /* compute diagonal terms of Q and mean rate mr = l/t */
1506      mod->mr = .0;
1507      For (i,ns)
1508        {
1509          sum=.0;
1510          For(j, ns) sum += mod->mat_Q[i*ns+j];
1511          mod->mat_Q[i*ns+i] = -sum;
1512          mod->mr += mod->pi[i] * sum;
1513        }
1514
1515      /* scale instantaneous rate matrix so that mu=1 */
1516      For (i,ns*ns) mod->mat_Q[i] /= mod->mr;
1517
1518
1519      /* compute eigenvectors/values */
1520      result = 0;
1521      if (result==eigen(1, mod->mat_Q,ns,dr,di,mod->mat_Vr, mod->mat_Vi, space))
1522        {
1523          /* compute inverse(Vr) into Vi */
1524          For (i,ns*ns) mod->mat_Vi[i] = mod->mat_Vr[i];
1525          Matinv(mod->mat_Vi, ns, ns);
1526         
1527          /* compute the diagonal terms of exp(D) */
1528          For (i,ns)  mod->vct_eDmr[i] = exp(dr[i]);
1529        }
1530      else
1531        {
1532          if (result==-1)
1533            printf("\n. Eigenvalues/vectors computation does not converge : computation cancelled");
1534          else if (result==1)
1535            printf("\n. Complex eigenvalues/vectors : computation cancelled");
1536        }
1537    }
1538
1539  mod->alpha_old  = mod->alpha;
1540  mod->kappa_old  = mod->kappa;
1541  mod->lambda_old = mod->lambda;
1542  mod->pinvar_old = mod->pinvar;
1543
1544
1545  free(dr);free(di);free(space);
1546}
1547
1548/*********************************************************/
1549
1550void Update_Qmat_GTR(model *mod)
1551{
1552  int result;
1553  int i,j,ns;
1554  double *di,*space,*mat_buff;
1555
1556  ns = 4;
1557 
1558  di    = (double *)mCalloc(ns,sizeof(double));
1559  space = (double *)mCalloc(2*ns,sizeof(double));
1560  mat_buff = (double *)mCalloc(ns*ns,sizeof(double));
1561
1562
1563  mod->mat_Q[0*4+1] = *(mod->rr_param[0])*mod->pi[1];
1564  mod->mat_Q[0*4+2] = *(mod->rr_param[1])*mod->pi[2];
1565  mod->mat_Q[0*4+3] = *(mod->rr_param[2])*mod->pi[3];
1566 
1567  mod->mat_Q[1*4+0] = *(mod->rr_param[0])*mod->pi[0];
1568  mod->mat_Q[1*4+2] = *(mod->rr_param[3])*mod->pi[2];
1569  mod->mat_Q[1*4+3] = *(mod->rr_param[4])*mod->pi[3];
1570
1571  mod->mat_Q[2*4+0] = *(mod->rr_param[1])*mod->pi[0];
1572  mod->mat_Q[2*4+1] = *(mod->rr_param[3])*mod->pi[1];
1573  mod->mat_Q[2*4+3] =              1.0*mod->pi[3];
1574
1575  mod->mat_Q[3*4+0] = *(mod->rr_param[2])*mod->pi[0];
1576  mod->mat_Q[3*4+1] = *(mod->rr_param[4])*mod->pi[1];
1577  mod->mat_Q[3*4+2] =              1.0*mod->pi[2];
1578
1579  mod->mat_Q[0*4+0] = -(*(mod->rr_param[0])*mod->pi[1]+*(mod->rr_param[1])*mod->pi[2]+*(mod->rr_param[2])*mod->pi[3]);
1580  mod->mat_Q[1*4+1] = -(*(mod->rr_param[0])*mod->pi[0]+*(mod->rr_param[3])*mod->pi[2]+*(mod->rr_param[4])*mod->pi[3]);
1581  mod->mat_Q[2*4+2] = -(*(mod->rr_param[1])*mod->pi[0]+*(mod->rr_param[3])*mod->pi[1]+1.0*mod->pi[3]);
1582  mod->mat_Q[3*4+3] = -(*(mod->rr_param[2])*mod->pi[0]+*(mod->rr_param[4])*mod->pi[1]+1.0*mod->pi[2]);
1583
1584
1585
1586  /* compute diagonal terms of Q and mean rate mr = l/t */
1587  mod->mr = .0;
1588  For (i,ns) mod->mr += mod->pi[i] * (-mod->mat_Q[i*ns+i]);
1589  For(i,ns*ns) mod->mat_Q[i] /= mod->mr;
1590/*   printf("mr=%f\n",mod->mr); */
1591 
1592  /* compute eigenvectors/values */
1593  result = 0;
1594  For(i,ns) For(j,ns) mat_buff[i*ns+j] = mod->mat_Q[i*ns+j];
1595 
1596  if (result==eigen(1,mat_buff,ns,mod->vct_ev,di,mod->mat_Vr, mod->mat_Vi, space))
1597    {
1598    /* compute inverse(Vr) into Vi */
1599      For (i,ns*ns) mod->mat_Vi[i] = mod->mat_Vr[i];
1600      Matinv(mod->mat_Vi, ns, ns);
1601     
1602      /* compute the diagonal terms of exp(D) */
1603      For (i,ns)  mod->vct_eDmr[i] = exp(mod->vct_ev[i]);
1604    }
1605  else
1606    {
1607      if (result==-1)
1608        printf("eigenvalues/vectors computation does not converge : computation cancelled");
1609      else if (result==1)
1610        printf("complex eigenvalues/vectors : computation cancelled");
1611    }
1612
1613 
1614  Free(di);
1615  Free(space);
1616  Free(mat_buff);
1617}
1618
1619/*********************************************************/
1620
1621void Translate_Custom_Mod_String(model *mod)
1622{
1623  int identity;
1624  int i,j;
1625  int *n_diff_param;
1626  int *mod_code;
1627
1628  mod_code = (int *)mCalloc(6,sizeof(int));
1629
1630  n_diff_param = &(mod->n_diff_rr_param);
1631
1632  *n_diff_param=0;
1633
1634  For(i,6)
1635    {
1636
1637      identity = -1;
1638
1639      For(j,i)
1640        {
1641          if((mod->custom_mod_string[i] == 
1642              mod->custom_mod_string[j]))
1643            {
1644              identity = j;
1645              break;
1646            }
1647        }
1648
1649      if(identity == -1)
1650        {
1651          mod->rr_param_num[*n_diff_param][0] = i;
1652          mod->n_rr_param_per_cat[*n_diff_param] = 1;
1653          mod_code[i] = *n_diff_param;
1654          (*n_diff_param)++;
1655        }
1656      else
1657        {
1658          mod_code[i] = mod_code[identity];
1659          mod->rr_param_num[mod_code[i]]
1660            [mod->n_rr_param_per_cat[mod_code[i]]] = i;
1661          mod->n_rr_param_per_cat[mod_code[i]] += 1;
1662        }
1663    }
1664
1665  Free(mod_code);
1666
1667}
1668
1669/*********************************************************/
1670
1671void Set_Model_Parameters(arbre *tree)
1672{
1673  double sum;
1674  int i;
1675
1676
1677  DiscreteGamma (tree->mod->r_proba, tree->mod->rr, tree->mod->alpha,
1678                 tree->mod->alpha,tree->mod->n_catg,0);
1679   
1680
1681  if((tree->mod->whichmodel < 10) && (tree->mod->s_opt->opt_bfreq))
1682    {
1683      sum = .0;
1684      For(i,4) sum += tree->mod->pi[i];
1685      For(i,4) 
1686        {
1687          tree->mod->pi[i] /= sum;
1688/*        printf("pi[%d]->%f\n",i+1,tree->mod->pi[i]); */
1689        }
1690    }
1691
1692
1693  if(tree->mod->whichmodel >= 7)
1694    {
1695        For(i,6) 
1696            {
1697                tree->mod->rr_param_values[i] /= tree->mod->rr_param_values[5];
1698            }
1699    }
1700
1701  if(tree->mod->update_eigen) 
1702    {
1703      if(tree->mod->datatype == NT)
1704        {
1705          if(tree->mod->whichmodel < 20) Update_Qmat_GTR(tree->mod);
1706        }
1707    }
1708
1709}
1710
1711/*********************************************************/
Note: See TracBrowser for help on using the repository browser.