← Back to branch summary

~ubuntu-branches/ubuntu/utopic/paml/utopic

« back to all changes in this revision

Viewing changes to Technical/Pt/testPMat.c

  • Committer: Bazaar Package Importer
  • Author(s): Pjotr Prins
  • Date: 2010-09-11 23:01:37 UTC
  • Revision ID: james.westby@ubuntu.com-20100911230137-jjf5d0blx5p0m9ba
Tags: upstream-4.4c
Import upstream version 4.4c

Show diffs side-by-side

added added

removed removed

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