source: branches/ali/GDE/PHYML20130708/phyml/src/models.c

Last change on this file was 10307, checked in by aboeckma, 11 years ago

added most recent version of phyml
File size: 34.8 KB
Line 
1/*
2
3PHYML :  a program that  computes maximum likelihood  phyLOGenies from
4DNA or AA homoLOGous sequences
5
6Copyright (C) Stephane Guindon. Oct 2003 onward
7
8All parts of  the source except where indicated  are distributed under
9the GNU public licence.  See http://www.opensource.org for details.
10
11*/
12
13#include "models.h"
14
15
16//////////////////////////////////////////////////////////////
17//////////////////////////////////////////////////////////////
18
19/* Handle any number of states (>1) */
20void PMat_JC69(phydbl l, int pos, phydbl *Pij, t_mod *mod)
21{
22  int ns;
23  int i,j;
24
25  ns = mod->ns;
26
27
28  For(i,ns) Pij[pos+ ns*i+i] = 1. - ((ns - 1.)/ns)*(1. - EXP(-ns*l/(ns - 1.)));
29  For(i,ns-1) 
30    for(j=i+1;j<ns;j++) 
31      {
32        Pij[pos+ ns*i+j] = (1./ns)*(1. - EXP(-ns*l/(ns - 1.)));
33        if(Pij[pos+ns*i+j] < SMALL_PIJ) Pij[pos+ns*i+j] = SMALL_PIJ;
34        Pij[pos+ ns*j+i] = Pij[pos+ ns*i+j];
35      }
36}
37
38
39//////////////////////////////////////////////////////////////
40//////////////////////////////////////////////////////////////
41
42void PMat_K80(phydbl l, phydbl kappa, int pos, phydbl *Pij)
43{
44  phydbl Ts,Tv,e1,e2,aux;
45  int i,j;
46  /*0 => A*/
47  /*1 => C*/
48  /*2 => G*/
49  /*3 => T*/
50
51  /* Ts -> transition*/
52  /* Tv -> transversion*/
53
54
55  aux = -2*l/(kappa+2);
56  e1 = (phydbl)EXP(aux *2);
57 
58  e2 = (phydbl)EXP(aux *(kappa+1));
59  Tv = .25*(1-e1);
60  Ts = .25*(1+e1-2*e2);
61
62  Pij[pos+ 4*0+0] = Pij[pos+ 4*1+1] = 
63  Pij[pos+ 4*2+2] = Pij[pos+ 4*3+3] = 1.-Ts-2.*Tv;
64
65  Pij[pos+ 4*0+1] = Pij[pos+ 4*1+0] = Tv;
66  Pij[pos+ 4*0+2] = Pij[pos+ 4*2+0] = Ts;
67  Pij[pos+ 4*0+3] = Pij[pos+ 4*3+0] = Tv;
68
69  Pij[pos+ 4*1+2] = Pij[pos+ 4*2+1] = Tv;
70  Pij[pos+ 4*1+3] = Pij[pos+ 4*3+1] = Ts;
71
72  Pij[pos+ 4*2+3] = Pij[pos+ 4*3+2] = Tv;
73
74  For(i,4) For(j,4)
75    if(Pij[pos + 4*i+j] < SMALL_PIJ) Pij[pos + 4*i+j] = SMALL_PIJ;
76
77}
78
79//////////////////////////////////////////////////////////////
80//////////////////////////////////////////////////////////////
81
82
83
84void PMat_TN93(phydbl l, t_mod *mod, int pos, phydbl *Pij)
85{
86  int i,j;
87  phydbl e1,e2,e3;
88  phydbl a1t,a2t,bt;
89  phydbl A,C,G,T,R,Y;
90  phydbl kappa1,kappa2;
91
92  A = mod->e_frq->pi->v[0]; C = mod->e_frq->pi->v[1]; G = mod->e_frq->pi->v[2]; T = mod->e_frq->pi->v[3];
93  R = A+G;  Y = T+C;
94
95  if(mod->kappa->v < .0) mod->kappa->v = 1.0e-5;
96
97  if((mod->whichmodel != F84) && (mod->whichmodel != TN93)) mod->lambda->v = 1.; 
98  else if(mod->whichmodel == F84)
99    {
100      mod->lambda->v = Get_Lambda_F84(mod->e_frq->pi->v,&mod->kappa->v);
101    }
102
103  kappa2 = mod->kappa->v*2./(1.+mod->lambda->v);
104  kappa1 = kappa2 * mod->lambda->v;
105
106 
107  bt = l/(2.*(A*G*kappa1+C*T*kappa2+R*Y));
108
109  a1t = kappa1;
110  a2t = kappa2;
111  a1t*=bt; a2t*=bt;
112
113  e1 = (phydbl)EXP(-a1t*R-bt*Y);
114  e2 = (phydbl)EXP(-a2t*Y-bt*R);
115  e3 = (phydbl)EXP(-bt);
116
117
118  /*A->A*/Pij[pos + 4*0+0] = A+Y*A/R*e3+G/R*e1; 
119  /*A->C*/Pij[pos + 4*0+1] = C*(1-e3);
120  /*A->G*/Pij[pos + 4*0+2] = G+Y*G/R*e3-G/R*e1;
121  /*A->T*/Pij[pos + 4*0+3] = T*(1-e3);
122
123  /*C->A*/Pij[pos + 4*1+0] = A*(1-e3);
124  /*C->C*/Pij[pos + 4*1+1] = C+R*C/Y*e3+T/Y*e2;
125  /*C->G*/Pij[pos + 4*1+2] = G*(1-e3);
126  /*C->T*/Pij[pos + 4*1+3] = T+R*T/Y*e3-T/Y*e2;
127
128  /*G->A*/Pij[pos + 4*2+0] = A+Y*A/R*e3-A/R*e1;
129  /*G->C*/Pij[pos + 4*2+1] = C*(1-e3);
130  /*G->G*/Pij[pos + 4*2+2] = G+Y*G/R*e3+A/R*e1;
131  /*G->T*/Pij[pos + 4*2+3] = T*(1-e3);
132
133  /*T->A*/Pij[pos + 4*3+0] = A*(1-e3);
134  /*T->C*/Pij[pos + 4*3+1] = C+R*C/Y*e3-C/Y*e2;
135  /*T->G*/Pij[pos + 4*3+2] = G*(1-e3);
136  /*T->T*/Pij[pos + 4*3+3] = T+R*T/Y*e3+C/Y*e2;
137 
138  For(i,4) For(j,4)
139    if(Pij[pos + 4*i+j] < SMALL_PIJ) Pij[pos + 4*i+j] = SMALL_PIJ;
140
141/*   /\*A->A*\/(*Pij)[0][0] = A+Y*A/R*e3+G/R*e1;  */
142/*   /\*A->C*\/(*Pij)[0][1] = C*(1-e3); */
143/*   /\*A->G*\/(*Pij)[0][2] = G+Y*G/R*e3-G/R*e1; */
144/*   /\*A->T*\/(*Pij)[0][3] = T*(1-e3); */
145
146/*   /\*C->A*\/(*Pij)[1][0] = A*(1-e3); */
147/*   /\*C->C*\/(*Pij)[1][1] = C+R*C/Y*e3+T/Y*e2; */
148/*   /\*C->G*\/(*Pij)[1][2] = G*(1-e3); */
149/*   /\*C->T*\/(*Pij)[1][3] = T+R*T/Y*e3-T/Y*e2; */
150
151/*   /\*G->A*\/(*Pij)[2][0] = A+Y*A/R*e3-A/R*e1; */
152/*   /\*G->C*\/(*Pij)[2][1] = C*(1-e3); */
153/*   /\*G->G*\/(*Pij)[2][2] = G+Y*G/R*e3+A/R*e1; */
154/*   /\*G->T*\/(*Pij)[2][3] = T*(1-e3); */
155
156/*   /\*T->A*\/(*Pij)[3][0] = A*(1-e3); */
157/*   /\*T->C*\/(*Pij)[3][1] = C+R*C/Y*e3-C/Y*e2; */
158/*   /\*T->G*\/(*Pij)[3][2] = G*(1-e3); */
159/*   /\*T->T*\/(*Pij)[3][3] = T+R*T/Y*e3+C/Y*e2; */
160 
161/*   For(i,4) For(j,4) */
162/*     if((*Pij)[i][j] < SMALL) (*Pij)[i][j] = SMALL; */
163
164}
165
166//////////////////////////////////////////////////////////////
167//////////////////////////////////////////////////////////////
168
169
170phydbl Get_Lambda_F84(phydbl *pi, phydbl *kappa)
171{
172  phydbl A,C,G,T,R,Y,lambda;
173  int kappa_has_changed;
174
175  A = pi[0]; C = pi[1]; G = pi[2]; T = pi[3];
176  R = A+G;  Y = T+C;
177
178  if(*kappa < .0) *kappa = 1.0e-5;
179
180  kappa_has_changed = NO;
181
182  do
183    {
184      lambda = (Y+(R-Y)/(2.*(*kappa)))/(R-(R-Y)/(2.*(*kappa)));
185     
186      if(lambda < .0)
187        {
188          *kappa += *kappa/10.;
189          kappa_has_changed = YES;
190        }
191    }while(lambda < .0);
192
193  if(kappa_has_changed)
194    {
195      PhyML_Printf("\n. WARNING: This transition/transversion ratio\n");
196      PhyML_Printf("  is impossible with these base frequencies!\n");
197      PhyML_Printf("  The ratio is now set to %.3f\n",*kappa);
198    }
199
200  return lambda;
201}
202
203
204//////////////////////////////////////////////////////////////
205//////////////////////////////////////////////////////////////
206
207
208
209/********************************************************************/
210/* void PMat_Empirical(phydbl l, t_mod *mod, phydbl ***Pij)         */
211/*                                                                  */
212/* Computes the substitution probability matrix                     */
213/* from the initial substitution rate matrix and frequency vector   */
214/* and one specific branch length                                   */
215/*                                                                  */
216/* input : l , branch length                                        */
217/* input : mod , choosen model parameters, qmat and pi             */
218/* ouput : Pij , substitution probability matrix                    */
219/*                                                                  */
220/* matrix P(l) is computed as follows :                             */
221/* P(l) = EXP(Q*t) , where :                                        */
222/*                                                                  */
223/*   Q = substitution rate matrix = Vr*D*inverse(Vr) , where :      */
224/*                                                                  */
225/*     Vr = right eigenvector matrix for Q                          */
226/*     D  = diagonal matrix of eigenvalues for Q                    */
227/*                                                                  */
228/*   t = time interval = l / mr , where :                           */
229/*                                                                  */
230/*     mr = mean rate = branch length/time interval                 */
231/*        = sum(i)(pi[i]*p(i->j)) , where :                         */
232/*                                                                  */
233/*       pi = state frequency vector                                */
234/*       p(i->j) = subst. probability from i to a different state   */
235/*               = -Q[ii] , as sum(j)(Q[ij]) +Q[ii] =0              */
236/*                                                                  */
237/* the Taylor development of EXP(Q*t) gives :                       */
238/* P(l) = Vr*EXP(D*t)        *inverse(Vr)                           */
239/*      = Vr*POW(EXP(D/mr),l)*inverse(Vr)                           */
240/*                                                                  */
241/* for performance we compute only once the following matrixes :    */
242/* Vr, inverse(Vr), EXP(D/mr)                                       */
243/* thus each time we compute P(l) we only have to :                 */
244/* make 20 times the operation POW()                                */
245/* make 2 20x20 matrix multiplications , that is :                  */
246/*   16000 = 2x20x20x20 times the operation *                       */
247/*   16000 = 2x20x20x20 times the operation +                       */
248/*   which can be reduced to (the central matrix being diagonal) :  */
249/*   8400 = 20x20 + 20x20x20 times the operation *                  */
250/*   8000 = 20x20x20 times the operation +                          */
251/********************************************************************/
252
253void PMat_Empirical(phydbl l, t_mod *mod, int pos, phydbl *Pij)
254{
255  int n = mod->ns;
256  int i, j, k;
257  phydbl *U,*V,*R;
258  phydbl *expt; 
259  phydbl *uexpt;
260
261  expt  = mod->eigen->e_val_im;
262  uexpt = mod->eigen->r_e_vect_im;
263  U     = mod->eigen->r_e_vect;
264  V     = mod->eigen->l_e_vect;
265  R     = mod->eigen->e_val; /* exponential of the eigen value matrix */
266
267  For(i,n) For(k,n) Pij[pos+mod->ns*i+k] = .0;
268
269  /* compute POW(EXP(D/mr),l) into mat_eDmrl */
270  For(k,n) expt[k] = (phydbl)POW(R[k],l);
271 
272  /* multiply Vr*POW(EXP(D/mr),l)*Vi into Pij */
273  For (i,n) For (k,n) uexpt[i*n+k] = U[i*n+k] * expt[k];
274
275  For (i,n) 
276    {
277      For (j,n) 
278        {
279          For(k,n)
280            {
281              Pij[pos+mod->ns*i+j] += (uexpt[i*n+k] * V[k*n+j]);
282            }
283/*        if(Pij[pos+mod->ns*i+j] < SMALL) Pij[pos+mod->ns*i+j] = SMALL; */
284          if(Pij[pos+mod->ns*i+j] < SMALL_PIJ) Pij[pos+mod->ns*i+j] = SMALL_PIJ;
285        }
286
287#ifndef PHYML
288      phydbl sum;
289      sum = .0;
290      For (j,n) sum += Pij[pos+mod->ns*i+j];
291      if((sum > 1.+.0001) || (sum < 1.-.0001))
292        {
293          PhyML_Printf("\n");
294          PhyML_Printf("\n. Q\n");
295          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",mod->eigen->q[i*n+j]); PhyML_Printf("\n"); }
296          PhyML_Printf("\n. U\n");
297          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",U[i*n+j]); PhyML_Printf("\n"); }
298          PhyML_Printf("\n");
299          PhyML_Printf("\n. V\n");
300          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",V[i*n+j]); PhyML_Printf("\n"); }
301          PhyML_Printf("\n");
302          PhyML_Printf("\n. Eigen\n");
303          For(i,n)  PhyML_Printf("%E ",expt[i]);
304          PhyML_Printf("\n");
305          PhyML_Printf("\n. Pij\n");
306          For(i,n) { For (j,n) PhyML_Printf("%f ",Pij[pos+mod->ns*i+j]); PhyML_Printf("\n"); }
307          PhyML_Printf("\n. sum = %f",sum);
308          if(mod->m4mod)
309            {
310              int i;
311              PhyML_Printf("\n. mod->m4mod->alpha = %f",mod->m4mod->alpha);
312              PhyML_Printf("\n. mod->m4mod->delta = %f",mod->m4mod->delta);
313              For(i,mod->m4mod->n_h)
314                {
315                  PhyML_Printf("\n. mod->m4mod->multipl[%d] = %f",i,mod->m4mod->multipl[i]);
316                }
317            }
318          PhyML_Printf("\n. l=%f",l);
319          PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
320          Warn_And_Exit("");
321        }
322#endif
323    }
324}
325
326//////////////////////////////////////////////////////////////
327//////////////////////////////////////////////////////////////
328
329
330void PMat_Gamma(phydbl l, t_mod *mod, int pos, phydbl *Pij)
331{
332  int n;
333  int i, j, k;
334  phydbl *U,*V,*R;
335  phydbl *expt; 
336  phydbl *uexpt;
337  phydbl shape;
338
339 
340  n     = mod->ns;
341  expt  = mod->eigen->e_val_im;
342  uexpt = mod->eigen->r_e_vect_im;
343  U     = mod->eigen->r_e_vect;
344  V     = mod->eigen->l_e_vect;
345  R     = mod->eigen->e_val; /* exponential of the eigen value matrix */
346 
347  if(mod->ras->n_catg == 1) shape = 1.E+4;
348  else                 shape = mod->ras->alpha->v;
349
350
351  For(i,n) For(k,n) Pij[pos+mod->ns*i+k] = .0;
352 
353  if(shape < 1.E-10) 
354    {
355      PhyML_Printf("\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
356      Warn_And_Exit("");
357    }
358
359  /* Formula 13.42, page 220 of Felsenstein's book ``Inferring Phylogenies'' */ 
360  For(k,n) expt[k] = POW(shape/(shape-LOG(R[k])*l),shape);
361
362  /* multiply Vr*expt*Vi into Pij */
363  For(i,n) For(k,n) uexpt[i*n+k] = U[i*n+k] * expt[k];
364
365  For (i,n) 
366    {
367      For (j,n) 
368        {
369          For(k,n)
370            {
371              Pij[pos+mod->ns*i+j] += (uexpt[i*n+k] * V[k*n+j]);
372            }
373          if(Pij[pos+mod->ns*i+j] < SMALL_PIJ) Pij[pos+mod->ns*i+j] = SMALL_PIJ;
374        }
375
376#ifdef DEBUG
377      phydbl sum;
378      sum = .0;
379      For (j,n) sum += Pij[pos+mod->ns*i+j];
380      if((sum > 1.+.0001) || (sum < 1.-.0001))
381        {
382          PhyML_Printf("\n");
383          PhyML_Printf("\n. Q\n");
384          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",mod->eigen->q[i*n+j]); PhyML_Printf("\n"); }
385          PhyML_Printf("\n. U\n");
386          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",U[i*n+j]); PhyML_Printf("\n"); }
387          PhyML_Printf("\n");
388          PhyML_Printf("\n. V\n");
389          For(i,n) { For(j,n) PhyML_Printf("%7.3f ",V[i*n+j]); PhyML_Printf("\n"); }
390          PhyML_Printf("\n");
391          PhyML_Printf("\n. Eigen\n");
392          For(i,n)  PhyML_Printf("%E ",expt[i]);
393          PhyML_Printf("\n");
394          PhyML_Printf("\n. Pij\n");
395          For(i,n) { For (j,n) PhyML_Printf("%f ",Pij[pos+mod->ns*i+j]); PhyML_Printf("\n"); }
396          PhyML_Printf("\n. sum = %f",sum);
397          if(mod->m4mod)
398            {
399              int i;
400              PhyML_Printf("\n. mod->m4mod->ras->alpha = %f",mod->m4mod->alpha);
401              PhyML_Printf("\n. mod->m4mod->delta = %f",mod->m4mod->delta);
402              For(i,mod->m4mod->n_h)
403                {
404                  PhyML_Printf("\n. mod->m4mod->multipl[%d] = %f",i,mod->m4mod->multipl[i]);
405                }
406            }
407          PhyML_Printf("\n. l=%f",l);
408          PhyML_Printf("\n. Err in file %s at line %d\n\n",__FILE__,__LINE__);
409          Warn_And_Exit("");
410        }
411#endif
412    }
413}
414
415//////////////////////////////////////////////////////////////
416//////////////////////////////////////////////////////////////
417
418
419void PMat_Zero_Br_Len(t_mod *mod, int pos, phydbl *Pij)
420{
421  int n = mod->ns;
422  int i, j;
423
424  For (i,n) For (j,n) Pij[pos+mod->ns*i+j] = .0;
425  For (i,n) Pij[pos+mod->ns*i+i] = 1.0;
426
427}
428
429//////////////////////////////////////////////////////////////
430//////////////////////////////////////////////////////////////
431
432
433void PMat(phydbl l, t_mod *mod, int pos, phydbl *Pij)
434{
435  /* Warning: l is never the og of branch length here */
436  if(l < 0.0)
437    {
438      PMat_Zero_Br_Len(mod,pos,Pij);
439    }
440  else
441    {
442      switch(mod->io->datatype)
443        {
444        case NT :
445          {
446            if(mod->use_m4mod)
447              {
448                PMat_Empirical(l,mod,pos,Pij);
449              }
450            else
451              {
452                if((mod->whichmodel == JC69) || 
453                   (mod->whichmodel == K80)) 
454                  {
455/*                  PMat_JC69(l,pos,Pij,mod); */
456                    PMat_K80(l,mod->kappa->v,pos,Pij);
457                  }
458                else
459                  {
460                    if(
461                       (mod->whichmodel == F81)   ||
462                       (mod->whichmodel == HKY85) ||
463                       (mod->whichmodel == F84)   ||
464                       (mod->whichmodel == TN93))
465                      {
466                        PMat_TN93(l,mod,pos,Pij);
467                      }
468                    else
469                      {
470                        PMat_Empirical(l,mod,pos,Pij);
471                      }
472                  }
473                break;
474              }
475          case AA : 
476            {
477              PMat_Empirical(l,mod,pos,Pij);
478              break;
479            }
480          default:
481            {
482              PMat_JC69(l,pos,Pij,mod);
483              break;
484/*            PhyML_Printf("\n. Not implemented yet.\n"); */
485/*            PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__); */
486/*            Warn_And_Exit(""); */
487/*            break; */
488            }
489          }
490        }
491    }
492}
493
494//////////////////////////////////////////////////////////////
495//////////////////////////////////////////////////////////////
496
497
498int GetDaa (phydbl *daa, phydbl *pi, char *file_name)
499{
500/* Get the amino acid distance (or substitution rate) matrix
501   (grantham, dayhoff, jones, etc).
502*/
503   FILE * fdaa;
504   int i,j, naa;
505   phydbl dmax,dmin;
506   phydbl sum; 
507   double val;
508
509   naa = 20;
510   dmax = .0;
511   dmin = 1.E+40;
512
513   fdaa = (FILE *)Openfile(file_name,0);
514
515   for(i=0; i<naa; i++) 
516     for(j=0; j<i; j++) 
517       {
518/*       if(fscanf(fdaa, "%lf", &daa[i*naa+j])) Exit("\n. err aaRatefile"); */
519         if(fscanf(fdaa, "%lf", &val)) Exit("\n. err aaRatefile");
520         daa[i*naa+j] = (phydbl)val;
521         daa[j*naa+i]=daa[i*naa+j];
522         if (dmax<daa[i*naa+j]) dmax=daa[i*naa+j];
523         if (dmin>daa[i*naa+j]) dmin=daa[i*naa+j];
524       }
525   
526   For(i,naa) 
527     {
528/*        if(fscanf(fdaa,"%lf",&pi[i])!=1) Exit("\n. err aaRatefile"); */
529       if(fscanf(fdaa,"%lf",&val)!=1) Exit("\n. err aaRatefile");
530       pi[i] = (phydbl)val;
531     }
532   sum = 0.0;
533   For(i, naa) sum += pi[i];
534   if (FABS(1-sum)>1e-4) {
535     PhyML_Printf("\nSum of freq. = %.6f != 1 in aaRateFile\n",sum); 
536     exit(-1);
537   }
538   
539   fclose (fdaa);
540
541   return (0);
542}
543
544//////////////////////////////////////////////////////////////
545//////////////////////////////////////////////////////////////
546
547
548void Update_Qmat_Generic(phydbl *rr, phydbl *pi, int ns, phydbl *qmat)
549{
550  int i,j;
551  phydbl sum,mr;
552
553  For(i,ns*ns) qmat[i] = .0;
554 
555  if(rr[(int)(ns*(ns-1)/2)-1] < 0.00001) 
556    {
557      PhyML_Printf("\n== rr[%d]=%f",(int)(ns*(ns-1)/2)-1,rr[(int)(ns*(ns-1)/2)-1]);
558      PhyML_Printf("\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
559      Warn_And_Exit("");
560    }
561
562  /* PhyML_Printf("\n"); */
563  /* For(i,(int)(ns*(ns-1)/2)) */
564  /*   { */
565  /*     PhyML_Printf("\n> rr %d = %f",i,rr[i]); */
566  /*   } */
567
568  For(i,(int)(ns*(ns-1)/2)) 
569    { 
570      rr[i] /= rr[(int)(ns*(ns-1)/2)-1];
571    }
572
573  /* Fill the non-diagonal parts */
574  For(i,ns)
575    {
576      for(j=i+1;j<ns;j++)
577        {
578          qmat[i*ns+j] = rr[MIN(i,j) * ns + MAX(i,j) -
579                            (MIN(i,j)+1+(int)POW(MIN(i,j)+1,2))/2];
580          qmat[j*ns+i] = qmat[i*ns+j];
581        }
582    }
583
584
585  /* Multiply by pi */
586  For(i,ns)
587    {
588      For(j,ns)
589        {
590          qmat[i*ns+j] *= pi[j];
591        }
592    }
593
594  /* Compute diagonal elements */
595  mr = .0;
596  For(i,ns)
597    {
598      sum = .0; 
599      For(j,ns) {sum += qmat[i*ns+j];}
600      qmat[i*ns+i] = -sum;
601      mr += sum * pi[i];
602    }
603 
604  /* For(i,ns) For(j,ns) qmat[i*ns+j] /= mr; */
605}
606
607//////////////////////////////////////////////////////////////
608//////////////////////////////////////////////////////////////
609
610
611void Update_Qmat_GTR(phydbl *rr, phydbl *rr_val, int *rr_num, phydbl *pi, phydbl *qmat)
612{
613  int i;
614  phydbl mr;
615 
616  For(i,6) rr[i] = rr_val[rr_num[i]];
617  For(i,6) 
618    if(rr[i] < 0.0) 
619      {
620        PhyML_Printf("\n. rr%d: %f",i,rr[i]);
621        PhyML_Printf("\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
622        Exit("");
623      }
624
625
626  For(i,6) rr[i] /= rr[5];
627
628  qmat[0*4+1] = (rr[0]*pi[1]);
629  qmat[0*4+2] = (rr[1]*pi[2]);
630  qmat[0*4+3] = (rr[2]*pi[3]);
631 
632  qmat[1*4+0] = (rr[0]*pi[0]);
633  qmat[1*4+2] = (rr[3]*pi[2]);
634  qmat[1*4+3] = (rr[4]*pi[3]);
635
636  qmat[2*4+0] = (rr[1]*pi[0]);
637  qmat[2*4+1] = (rr[3]*pi[1]);
638  qmat[2*4+3] = (rr[5]*pi[3]);
639
640  qmat[3*4+0] = (rr[2]*pi[0]);
641  qmat[3*4+1] = (rr[4]*pi[1]);
642  qmat[3*4+2] = (rr[5]*pi[2]);
643
644  qmat[0*4+0] = -(rr[0]*pi[1]+rr[1]*pi[2]+rr[2]*pi[3]);
645  qmat[1*4+1] = -(rr[0]*pi[0]+rr[3]*pi[2]+rr[4]*pi[3]);
646  qmat[2*4+2] = -(rr[1]*pi[0]+rr[3]*pi[1]+rr[5]*pi[3]);
647  qmat[3*4+3] = -(rr[2]*pi[0]+rr[4]*pi[1]+rr[5]*pi[2]);
648
649  /* compute diagonal terms of Q and mean rate mr = l/t */
650  mr = .0;
651  For (i,4) mr += pi[i] * (-qmat[i*4+i]);
652  For(i,16) qmat[i] /= mr;
653
654  /* int j; */
655  /* printf("\n"); */
656  /* printf("\n. rr -- "); */
657  /* For(i,5) printf(" %15f ",rr[i]); */
658  /* printf("\n"); */
659  /* For(i,4) */
660  /*   { */
661  /*     printf("\n. [%15f] \t ",pi[i]); */
662  /*     For(j,4) */
663  /*    { */
664  /*      printf(" %15f ",qmat[i*4+j]); */
665  /*    } */
666  /*   } */
667}
668
669//////////////////////////////////////////////////////////////
670//////////////////////////////////////////////////////////////
671
672
673void Update_Qmat_HKY(phydbl kappa, phydbl *pi, phydbl *qmat)
674{
675  int i;
676  phydbl mr;
677 
678  /* A -> C */ qmat[0*4+1] = (phydbl)(pi[1]);
679  /* A -> G */ qmat[0*4+2] = (phydbl)(kappa*pi[2]);
680  /* A -> T */ qmat[0*4+3] = (phydbl)(pi[3]);
681 
682  /* C -> A */ qmat[1*4+0] = (phydbl)(pi[0]);
683  /* C -> G */ qmat[1*4+2] = (phydbl)(pi[2]);
684  /* C -> T */ qmat[1*4+3] = (phydbl)(kappa*pi[3]);
685
686  /* G -> A */ qmat[2*4+0] = (phydbl)(kappa*pi[0]);
687  /* G -> C */ qmat[2*4+1] = (phydbl)(pi[1]);
688  /* G -> T */ qmat[2*4+3] = (phydbl)(pi[3]);
689
690  /* T -> A */ qmat[3*4+0] = (phydbl)(pi[0]);
691  /* T -> C */ qmat[3*4+1] = (phydbl)(kappa*pi[1]);
692  /* T -> G */ qmat[3*4+2] = (phydbl)(pi[2]);
693
694  qmat[0*4+0] = (phydbl)(-(qmat[0*4+1]+qmat[0*4+2]+qmat[0*4+3]));
695  qmat[1*4+1] = (phydbl)(-(qmat[1*4+0]+qmat[1*4+2]+qmat[1*4+3]));
696  qmat[2*4+2] = (phydbl)(-(qmat[2*4+0]+qmat[2*4+1]+qmat[2*4+3]));
697  qmat[3*4+3] = (phydbl)(-(qmat[3*4+0]+qmat[3*4+1]+qmat[3*4+2]));
698
699  /* compute diagonal terms of Q and mean rate mr = l/t */
700  mr = .0;
701  For (i,4) mr += pi[i] * (-qmat[i*4+i]);
702  For(i,16) qmat[i] /= mr;
703}
704
705//////////////////////////////////////////////////////////////
706//////////////////////////////////////////////////////////////
707
708
709void Translate_Custom_Mod_String(t_mod *mod)
710{
711  int i,j;
712
713  For(i,6) mod->r_mat->n_rr_per_cat->v[i] = 0;
714
715  mod->r_mat->n_diff_rr = 0;
716
717  For(i,6)
718    {
719      For(j,i)
720        {
721          if((mod->custom_mod_string->s[i] == mod->custom_mod_string->s[j]))
722            {
723              break;
724            }
725        }
726
727      if(i == j)
728        {
729          mod->r_mat->rr_num->v[i] = mod->r_mat->n_diff_rr;
730          mod->r_mat->n_diff_rr++;
731        }
732      else
733        {
734          mod->r_mat->rr_num->v[i] = mod->r_mat->rr_num->v[j];
735        }
736
737      mod->r_mat->n_rr_per_cat->v[mod->r_mat->rr_num->v[j]]++;
738    }
739
740/*   PhyML_Printf("\n"); */
741/*   For(i,6) PhyML_Printf("%d ",mod->rr_param_num[i]); */
742/*   For(i,mod->n_diff_rr_param) PhyML_Printf("\n. Class %d size %d",i+1,mod->r_mat->n_rr_param_per_cat[i]); */
743}
744
745//////////////////////////////////////////////////////////////
746//////////////////////////////////////////////////////////////
747
748// Update rate across sites distribution settings.
749
750void Update_RAS(t_mod *mod)
751{
752  phydbl sum;
753  int i;
754
755  if(mod->ras->free_mixt_rates == NO) DiscreteGamma(mod->ras->gamma_r_proba->v, 
756                                                    mod->ras->gamma_rr->v, 
757                                                    mod->ras->alpha->v, 
758                                                    mod->ras->alpha->v, 
759                                                    mod->ras->n_catg, 
760                                                    mod->ras->gamma_median);
761  else
762    {
763
764#if (!defined PHYML)
765
766      Qksort(mod->ras->gamma_r_proba_unscaled->v,NULL,0,mod->ras->n_catg-1); // Unscaled class frequencies sorted in increasing order
767
768      // Update class frequencies
769      For(i,mod->ras->n_catg)
770        {
771          if(!i)
772            mod->ras->gamma_r_proba->v[i] =
773              mod->ras->gamma_r_proba_unscaled->v[i] / (mod->ras->gamma_r_proba_unscaled->v[mod->ras->n_catg-1]);
774          else
775            mod->ras->gamma_r_proba->v[i] =
776              (mod->ras->gamma_r_proba_unscaled->v[i] - mod->ras->gamma_r_proba_unscaled->v[i-1]) /
777              (mod->ras->gamma_r_proba_unscaled->v[mod->ras->n_catg-1]) ;
778        }
779
780#else
781
782      sum = 0.0;
783      For(i,mod->ras->n_catg) sum += mod->ras->gamma_r_proba_unscaled->v[i];
784      For(i,mod->ras->n_catg) mod->ras->gamma_r_proba->v[i] = mod->ras->gamma_r_proba_unscaled->v[i] / sum;
785
786
787#endif
788
789      do
790        {
791          sum = .0;
792          For(i,mod->ras->n_catg)
793            {
794              if(mod->ras->gamma_r_proba->v[i] < 0.01) mod->ras->gamma_r_proba->v[i]=0.01;
795              if(mod->ras->gamma_r_proba->v[i] > 0.99) mod->ras->gamma_r_proba->v[i]=0.99;
796              sum += mod->ras->gamma_r_proba->v[i];
797            }
798          For(i,mod->ras->n_catg) mod->ras->gamma_r_proba->v[i]/=sum;
799        }
800      while((sum > 1.01) || (sum < 0.99));
801
802
803      // Update class rates
804      sum = .0;
805      For(i,mod->ras->n_catg) sum += mod->ras->gamma_r_proba->v[i] * mod->ras->gamma_rr_unscaled->v[i];
806
807      if(mod->ras->normalise_rr == YES)
808        For(i,mod->ras->n_catg) 
809          mod->ras->gamma_rr->v[i] = mod->ras->gamma_rr_unscaled->v[i]/sum;
810      else
811        For(i,mod->ras->n_catg) 
812          mod->ras->gamma_rr->v[i] = mod->ras->gamma_rr_unscaled->v[i] * mod->ras->free_rate_mr->v;
813
814      /* printf("\n"); */
815      /* For(i,mod->ras->n_catg) */
816      /*   printf("\nx %12f %12f xx %12f %12f", */
817      /*          mod->ras->gamma_r_proba->v[i], */
818      /*          mod->ras->gamma_rr->v[i], */
819      /*          mod->ras->gamma_r_proba_unscaled->v[i], */
820      /*          mod->ras->gamma_rr_unscaled->v[i]); */
821    } 
822
823}
824
825//////////////////////////////////////////////////////////////
826//////////////////////////////////////////////////////////////
827
828void Update_Efrq(t_mod *mod)
829{
830  phydbl sum;
831  int i;
832
833  if((mod->io->datatype == NT) && (mod->s_opt->opt_state_freq))
834    {
835      sum = .0;
836      For(i,mod->ns) sum += FABS(mod->e_frq->pi_unscaled->v[i]);
837      For(i,mod->ns) mod->e_frq->pi->v[i] = FABS(mod->e_frq->pi_unscaled->v[i])/sum;
838     
839      do
840        {
841          sum = .0;
842          For(i,mod->ns)
843            {
844              if(mod->e_frq->pi->v[i] < 0.01) mod->e_frq->pi->v[i]=0.01;
845              if(mod->e_frq->pi->v[i] > 0.99) mod->e_frq->pi->v[i]=0.99;
846              sum += mod->e_frq->pi->v[i];
847            }
848          For(i,mod->ns) mod->e_frq->pi->v[i]/=sum;
849        }
850      while((sum > 1.01) || (sum < 0.99));
851    }
852 
853 
854
855}
856
857//////////////////////////////////////////////////////////////
858//////////////////////////////////////////////////////////////
859
860void Set_Model_Parameters(t_mod *mod)
861{
862  Update_RAS(mod);
863  Update_Efrq(mod);
864  Update_Eigen(mod);
865}
866
867//////////////////////////////////////////////////////////////
868//////////////////////////////////////////////////////////////
869
870void Update_Eigen(t_mod *mod)
871{
872  int result, n_iter;
873  phydbl scalar;
874  int i;
875 
876
877  if(mod->is_mixt_mod) 
878    {
879      MIXT_Update_Eigen(mod);
880      return;
881    }
882
883  if(mod->update_eigen) 
884    {
885      if(mod->use_m4mod == NO)
886        {
887          if(mod->io->datatype == NT)
888            {
889              if(mod->whichmodel == GTR)
890                Update_Qmat_GTR(mod->r_mat->rr->v, mod->r_mat->rr_val->v, mod->r_mat->rr_num->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
891              else if(mod->whichmodel == CUSTOM) 
892                Update_Qmat_GTR(mod->r_mat->rr->v, mod->r_mat->rr_val->v, mod->r_mat->rr_num->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
893              else if(mod->whichmodel == HKY85) 
894                Update_Qmat_HKY(mod->kappa->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
895              else /* Any other nucleotide-based t_mod */
896                Update_Qmat_HKY(mod->kappa->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
897            }
898        }
899      else 
900        {
901          M4_Update_Qmat(mod->m4mod,mod);
902        }
903
904      scalar   = 1.0;
905      n_iter   = 0;
906      result   = 0;
907           
908      For(i,mod->ns*mod->ns) mod->r_mat->qmat_buff->v[i] = mod->r_mat->qmat->v[i];
909
910      /* compute eigenvectors/values */
911      /*       if(!EigenRealGeneral(mod->eigen->size,mod->r_mat->qmat,mod->eigen->e_val, */
912      /*                          mod->eigen->e_val_im, mod->eigen->r_e_vect, */
913      /*                          mod->eigen->space_int,mod->eigen->space)) */
914     
915      if(!Eigen(1,mod->r_mat->qmat_buff->v,mod->eigen->size,mod->eigen->e_val,
916                mod->eigen->e_val_im,mod->eigen->r_e_vect,
917                mod->eigen->r_e_vect_im,mod->eigen->space))
918        {
919          /* compute inverse(Vr) into Vi */
920          For (i,mod->ns*mod->ns) mod->eigen->l_e_vect[i] = mod->eigen->r_e_vect[i];
921          while(!Matinv(mod->eigen->l_e_vect, mod->eigen->size, mod->eigen->size,YES))
922            {
923              PhyML_Printf("\n. Trying Q<-Q*scalar and then Root<-Root/scalar to fix this...\n");
924              scalar += scalar / 3.;
925              For(i,mod->eigen->size*mod->eigen->size) mod->r_mat->qmat_buff->v[i]  = mod->r_mat->qmat->v[i];
926              For(i,mod->eigen->size*mod->eigen->size) mod->r_mat->qmat_buff->v[i] *= scalar;
927              result = Eigen(1,mod->r_mat->qmat_buff->v,mod->eigen->size,mod->eigen->e_val,
928                             mod->eigen->e_val_im,mod->eigen->r_e_vect,
929                             mod->eigen->r_e_vect_im,mod->eigen->space);
930              if (result == -1)
931                Exit("\n. Eigenvalues/vectors computation did not converge: computation cancelled\n");
932              else if (result == 1)
933                Exit("\n. Complex eigenvalues/vectors: computation cancelled\n");
934             
935              For (i,mod->eigen->size*mod->eigen->size) mod->eigen->l_e_vect[i] = mod->eigen->r_e_vect[i];
936              n_iter++;
937              if(n_iter > 100) Exit("\n. Cannot work out eigen vectors\n");
938            };
939          For(i,mod->eigen->size) mod->eigen->e_val[i] /= scalar;
940
941          /* compute the diagonal terms of EXP(D) */
942          For(i,mod->ns) mod->eigen->e_val[i] = (phydbl)EXP(mod->eigen->e_val[i]);
943
944
945          /* int j; */
946          /* double *U,*V,*R; */
947          /* double *expt; */
948          /* double *uexpt; */
949          /* int n; */
950
951          /* expt  = mod->eigen->e_val_im; */
952          /* uexpt = mod->eigen->r_e_vect_im; */
953          /* U     = mod->eigen->r_e_vect; */
954          /* V     = mod->eigen->l_e_vect; */
955          /* R     = mod->eigen->e_val; /\* exponential of the eigen value matrix *\/ */
956          /* n     = mod->ns; */
957
958          /* PhyML_Printf("\n"); */
959          /* PhyML_Printf("\n. Q\n"); */
960          /* For(i,n) { For(j,n) PhyML_Printf("%7.3f ",mod->eigen->q[i*n+j]); PhyML_Printf("\n"); } */
961          /* PhyML_Printf("\n. U\n"); */
962          /* For(i,n) { For(j,n) PhyML_Printf("%7.3f ",U[i*n+j]); PhyML_Printf("\n"); } */
963          /* PhyML_Printf("\n"); */
964          /* PhyML_Printf("\n. V\n"); */
965          /* For(i,n) { For(j,n) PhyML_Printf("%7.3f ",V[i*n+j]); PhyML_Printf("\n"); } */
966          /* PhyML_Printf("\n"); */
967          /* PhyML_Printf("\n. Eigen\n"); */
968          /* For(i,n)  PhyML_Printf("%E ",mod->eigen->e_val[i]); */
969          /* PhyML_Printf("\n"); */
970         
971/*        Exit("\n"); */
972
973        }
974      else
975        {
976          PhyML_Printf("\n. Eigenvalues/vectors computation does not converge : computation cancelled");
977          Warn_And_Exit("\n");
978        }
979    }
980
981}
982
983//////////////////////////////////////////////////////////////
984//////////////////////////////////////////////////////////////
985
986
987void Switch_From_M4mod_To_Mod(t_mod *mod)
988{
989  int i;
990
991  mod->use_m4mod = 0;
992  mod->ns = mod->m4mod->n_o;
993  For(i,mod->ns) mod->e_frq->pi->v[i] = mod->m4mod->o_fq[i];
994  mod->eigen->size = mod->ns;
995  Switch_Eigen(YES,mod);
996}
997
998//////////////////////////////////////////////////////////////
999//////////////////////////////////////////////////////////////
1000
1001
1002void Switch_From_Mod_To_M4mod(t_mod *mod)
1003{
1004  int i;
1005  mod->use_m4mod = 1;
1006  mod->ns = mod->m4mod->n_o * mod->m4mod->n_h;
1007  For(i,mod->ns) mod->e_frq->pi->v[i] = mod->m4mod->o_fq[i%mod->m4mod->n_o] * mod->m4mod->h_fq[i/mod->m4mod->n_o];
1008  mod->eigen->size = mod->ns;
1009  Switch_Eigen(YES,mod);
1010}
1011
1012//////////////////////////////////////////////////////////////
1013//////////////////////////////////////////////////////////////
1014
1015
1016phydbl General_Dist(phydbl *F, t_mod *mod, eigen *eigen_struct)
1017{
1018  phydbl *pi,*mod_pi;
1019  int i,j,k;
1020  phydbl dist;
1021  phydbl sum;
1022  phydbl sum_ev;
1023  phydbl *F_phydbl;
1024
1025
1026  /* TO DO : call eigen decomposition function for all nt models */
1027
1028  F_phydbl = (phydbl *)mCalloc(eigen_struct->size*eigen_struct->size,sizeof(phydbl));
1029  pi       = (phydbl *)mCalloc(eigen_struct->size,sizeof(phydbl));
1030  mod_pi   = (phydbl *)mCalloc(eigen_struct->size,sizeof(phydbl));
1031
1032  For(i,mod->ns) mod_pi[i] = mod->e_frq->pi->v[i];
1033
1034  sum = .0;
1035  For(i,eigen_struct->size) 
1036    {
1037      For(j,eigen_struct->size)
1038        {
1039          pi[i] += (F[eigen_struct->size*i+j] + F[eigen_struct->size*j+i])/2.;
1040          sum += F[eigen_struct->size*i+j];
1041        }
1042    }
1043 
1044  Make_Symmetric(&F,eigen_struct->size);
1045  Divide_Mat_By_Vect(&F,mod->e_frq->pi->v,eigen_struct->size);
1046
1047  /* Eigen decomposition of pi^{-1} x F */
1048  For(i,eigen_struct->size) For(j,eigen_struct->size) F_phydbl[eigen_struct->size*i+j] = F[eigen_struct->size*i+j];
1049
1050  if(Eigen(1,F_phydbl,mod->eigen->size,mod->eigen->e_val,
1051           mod->eigen->e_val_im,mod->eigen->r_e_vect,
1052           mod->eigen->r_e_vect_im,mod->eigen->space))
1053    {
1054      For(i,mod->ns) mod->e_frq->pi->v[i] = mod_pi[i];
1055      Update_Qmat_GTR(mod->r_mat->rr->v, mod->r_mat->rr_val->v, mod->r_mat->rr_num->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
1056      Free(pi); 
1057      Free(mod_pi); 
1058      return -1.;
1059    }
1060
1061  /* Get the left eigen vector of pi^{-1} x F */
1062  For(i,eigen_struct->size*eigen_struct->size) eigen_struct->l_e_vect[i] = eigen_struct->r_e_vect[i];
1063  if(!Matinv(eigen_struct->l_e_vect,eigen_struct->size,eigen_struct->size,YES)<0) 
1064    {
1065      For(i,mod->ns) mod->e_frq->pi->v[i] = mod_pi[i];
1066      Update_Qmat_GTR(mod->r_mat->rr->v, mod->r_mat->rr_val->v, mod->r_mat->rr_num->v, mod->e_frq->pi->v, mod->r_mat->qmat->v);
1067      Free(pi); 
1068      Free(mod_pi); 
1069      return -1.;
1070    }
1071
1072  /* LOG of eigen values */
1073  For(i,eigen_struct->size) 
1074    {
1075/*       if(eigen_struct->e_val[i] < 0.0) eigen_struct->e_val[i] = 0.0001; */
1076      eigen_struct->e_val[i] = (phydbl)LOG(eigen_struct->e_val[i]);
1077     }
1078 
1079  /* Matrix multiplications LOG(pi^{-1} x F) */
1080  For(i,eigen_struct->size) For(j,eigen_struct->size)
1081    eigen_struct->r_e_vect[eigen_struct->size*i+j] = 
1082    eigen_struct->r_e_vect[eigen_struct->size*i+j] * 
1083    eigen_struct->e_val[j];
1084  For(i,eigen_struct->size) For(j,eigen_struct->size) F[eigen_struct->size*i+j] = .0;
1085  For(i,eigen_struct->size) For(j,eigen_struct->size) For(k,eigen_struct->size)
1086    F[eigen_struct->size*i+j] += eigen_struct->r_e_vect[eigen_struct->size*i+k] * eigen_struct->l_e_vect[eigen_struct->size*k+j];
1087
1088
1089  /* Trace */
1090  dist = .0;
1091  For(i,eigen_struct->size) dist+=F[eigen_struct->size*i+i];
1092
1093  sum_ev = .0;
1094  For(i,mod->ns) sum_ev += mod->eigen->e_val[i];
1095
1096/*   dist /= sum_ev; */
1097  dist /= -4.;
1098
1099 
1100/*   For(i,mod->ns) mod->e_frq->pi->v[i] = mod_pi[i]; */
1101/*   Update_Qmat_GTR(mod); */
1102  Free(pi); 
1103  Free(mod_pi); 
1104  Free(F_phydbl);
1105
1106  if(isnan(dist)) return -1.;
1107  return dist;
1108}
1109
1110//////////////////////////////////////////////////////////////
1111//////////////////////////////////////////////////////////////
1112
1113phydbl GTR_Dist(phydbl *F, phydbl alpha, eigen *eigen_struct)
1114{
1115  phydbl *pi;
1116  int i,j,k;
1117  phydbl dist;
1118  phydbl sum;
1119  phydbl *F_phydbl;
1120
1121  pi       = (phydbl *)mCalloc(eigen_struct->size,sizeof(phydbl));
1122  F_phydbl = (phydbl *)mCalloc(eigen_struct->size*eigen_struct->size,sizeof(phydbl));
1123
1124/*   /\* Waddell and Steel's example *\/ */
1125/*   F[4*0+0] = 1415./4898.; F[4*0+1] = 8./4898.;    F[4*0+2] = 55./4898.;  F[4*0+3] = 2./4898.; */
1126/*   F[4*1+0] = 4./4898.;    F[4*1+1] = 1371./4898.; F[4*1+2] = 1./4898.;   F[4*1+3] = 144./4898.; */
1127/*   F[4*2+0] = 73./4898.;   F[4*2+1] = 0./4898.;    F[4*2+2] = 578./4898.; F[4*2+3] = 0./4898.; */
1128/*   F[4*3+0] = 3./4898.;    F[4*3+1] = 117./4898.;  F[4*3+2] = 1./4898.;   F[4*3+3] = 1126./4898.; */
1129
1130
1131  For(i,eigen_struct->size) 
1132    {
1133      For(j,eigen_struct->size)
1134        {
1135          pi[i] += (F[eigen_struct->size*i+j] + F[eigen_struct->size*j+i])/2.;
1136          sum += F[eigen_struct->size*i+j];
1137        }
1138    }
1139
1140/* /\*   Jukes and Cantor correction *\/ */
1141/*   sum = .0; */
1142/*   For(i,eigen_struct->size) sum += F[eigen_struct->size*i+i]; */
1143/*   sum = 1.-sum; */
1144/*   For(i,eigen_struct->size*eigen_struct->size) F[i] = sum/12.; */
1145/*   For(i,eigen_struct->size) F[eigen_struct->size*i+i] = (1.-sum)/4.; */
1146/*   For(i,eigen_struct->size) pi[i] = 1./(phydbl)eigen_struct->size; */
1147
1148
1149  Make_Symmetric(&F,eigen_struct->size);
1150  Divide_Mat_By_Vect(&F,pi,eigen_struct->size);
1151
1152
1153  /* Eigen decomposition of pi^{-1} x F */
1154  For(i,eigen_struct->size) For(j,eigen_struct->size) F_phydbl[eigen_struct->size*i+j] = F[eigen_struct->size*i+j];
1155  if(Eigen(1,F_phydbl,eigen_struct->size,eigen_struct->e_val,
1156           eigen_struct->e_val_im,eigen_struct->r_e_vect,
1157           eigen_struct->r_e_vect_im,eigen_struct->space))
1158    {
1159      Free(pi); 
1160      return -1.;
1161    }
1162
1163  /* Get the left eigen vector of pi^{-1} x F */
1164  For(i,eigen_struct->size*eigen_struct->size) eigen_struct->l_e_vect[i] = eigen_struct->r_e_vect[i];
1165  if(!Matinv(eigen_struct->l_e_vect,eigen_struct->size,eigen_struct->size,YES)<0) {Free(pi); return -1.;}
1166
1167  /* Equation (3) + inverse of the moment generating function for the gamma distribution (see Waddell & Steel, 1997) */
1168  For(i,eigen_struct->size) 
1169    {
1170      if(eigen_struct->e_val[i] < 0.0) 
1171        {
1172          eigen_struct->e_val[i] = 0.0001;
1173        }
1174      if(alpha < .0)
1175        eigen_struct->e_val[i] = (phydbl)LOG(eigen_struct->e_val[i]);
1176      else
1177        eigen_struct->e_val[i] = alpha * (1. - (phydbl)POW(eigen_struct->e_val[i],-1./alpha));
1178     }
1179 
1180  /* Matrix multiplications pi x LOG(pi^{-1} x F) */
1181  For(i,eigen_struct->size) For(j,eigen_struct->size)
1182    eigen_struct->r_e_vect[eigen_struct->size*i+j] = 
1183    eigen_struct->r_e_vect[eigen_struct->size*i+j] * eigen_struct->e_val[j];
1184  For(i,eigen_struct->size) For(j,eigen_struct->size) F[eigen_struct->size*i+j] = .0;
1185  For(i,eigen_struct->size) For(j,eigen_struct->size) For(k,eigen_struct->size)
1186    F[eigen_struct->size*i+j] += eigen_struct->r_e_vect[eigen_struct->size*i+k] * eigen_struct->l_e_vect[eigen_struct->size*k+j];
1187  For(i,eigen_struct->size) For(j,eigen_struct->size) F[eigen_struct->size*i+j] *= pi[i];
1188
1189  /* Trace */
1190  dist = .0;
1191  For(i,eigen_struct->size) dist-=F[eigen_struct->size*i+i];
1192
1193/*   PhyML_Printf("\nDIST = %f\n",dist); Exit("\n"); */
1194 
1195  Free(pi);
1196  Free(F_phydbl);
1197
1198  if(isnan(dist)) return -1.;
1199  return dist;
1200}
1201
1202
Note: See TracBrowser for help on using the repository browser.