Technical/Pt/testPMat.c

   1 /* testPMat.c\r
   2 \r
   3    November 2003, Z Yang\r
   4 \r
   5    This small program tests two methods or routines for calculating the \r
   6    transition probability matrix for a time reversible Markov rate matrix.  \r
   7    The first is the routine matexp(), which uses repeated squaring and works \r
   8    on a general matrix (not necesarily time reversible).  For large distances, \r
   9    more squaring should be used, but I have not tested this extensively.  You\r
  10    can change TimeSquare to see its effect.  \r
  11    The second method, PMatQRev, uses a routine for eigen values and eigen \r
  12    vectors for a symmetrical matrix.  It involves some copying and moving \r
  13    around when some frequencies are 0.  \r
  14    For each algorithm, this program generates random frequencies and random rate \r
  15    matrix and then calls an algorithm to calculate P(t) = exp(Q*t).\r
  16 \r
  17      cl -O2 -Ot -W4 testPMat.c tools.c\r
  18      cc -O3 testPMat.c tools.c -lm\r
  19      testPMat\r
  20 */\r
  21 \r
  22 #include "paml.h"\r
  23 \r
  24 int main(void)\r
  25 {\r
  26    int n=400, noisy=0, i,j;\r
  27    int nr=10, ir, TimeSquare=10, algorithm; /* TimeSquare should be larger for large t */\r
  28    double t=5, *Q, *pi, *space, s;\r
  29    char timestr[96], *AlgStr[2]={"repeated squaring", "eigensolution"};\r
  30    \r
  31    if((Q=(double*)malloc(n*n*5*sizeof(double))) ==NULL) error2("oom");\r
  32    pi=Q+n*n; space=pi+n;\r
  33 \r
  34    for(algorithm=0; algorithm<2; algorithm++) {\r
  35       starttime();\r
  36       SetSeed(1234567);\r
  37       for (i=0; i<n; i++)  pi[i]=rndu();\r
  38       s=sum(pi,n);\r
  39       for (i=0; i<n; i++)  pi[i]/=s;\r
  40 \r
  41       for(ir=0; ir<nr; ir++) {\r
  42          printf("Replicate %d/%d ", ir+1,nr);\r
  43 \r
  44            for (i=0; i<n; i++)  \r
  45             for (j=0,Q[i*n+i]=0; j<i; j++)\r
  46                Q[i*n+j]=Q[j*n+i] = square(rndu());\r
  47          for (i=0; i<n; i++)\r
  48             for (j=0; j<n; j++)\r
  49                Q[i*n+j] *= pi[j];\r
  50          for(i=0,s=0; i<n; i++)  {  /* rescaling Q so that average rate is 1 */\r
  51             Q[i*n+i]=0; Q[i*n+i]=-sum(Q+i*n, n); \r
  52             s-=pi[i]*Q[i*n+i];\r
  53          }\r
  54 \r
  55          if(noisy) {\r
  56             matout(stdout, pi, 1, n);\r
  57             matout(stdout, Q, n, n);\r
  58          }\r
  59 \r
  60          if(algorithm==0) \r
  61             matexp(Q, 1, n, TimeSquare, space);\r
  62          else \r
  63             PMatQRev(Q, pi, 1, n, space);\r
  64          printf("%s, time: %s\n", AlgStr[algorithm], printtime(timestr));\r
  65          if(noisy) \r
  66             matout(stdout, Q, n, n);\r
  67       }\r
  68    }\r
  69    return (0);\r
  70 }\r