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

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