ARB
trnsprob.cxx
Go to the documentation of this file.
1 // =============================================================
2 /* */
3 // File : trnsprob.c
4 // Purpose :
5 /* */
6 // Institute of Microbiology (Technical University Munich)
7 // http://www.arb-home.de/
8 /* */
9 // =============================================================
10 
11 #include "trnsprob.h"
12 #include "defs.h"
13 #include "mem.h"
14 
15 #include <arb_simple_assert.h>
16 #include <stdio.h>
17 
18 #define ih_assert(bed) arb_assert(bed)
19 #define EPS 0.00001
20 
21 //============================================================================
22 
23 static void identity(double **i) {
24  i[T][T]=1; i[T][C]=0; i[T][A]=0; i[T][G]=0;
25  i[C][T]=0; i[C][C]=1; i[C][A]=0; i[C][G]=0;
26  i[A][T]=0; i[A][C]=0; i[A][A]=1; i[A][G]=0;
27  i[G][T]=0; i[G][C]=0; i[G][A]=0; i[G][G]=1;
28 }
29 
30 //............................................................................
31 
32 static void copy(double **i,double **j) {
33  i[T][T]=j[T][T]; i[T][C]=j[T][C]; i[T][A]=j[T][A]; i[T][G]=j[T][G];
34  i[C][T]=j[C][T]; i[C][C]=j[C][C]; i[C][A]=j[C][A]; i[C][G]=j[C][G];
35  i[A][T]=j[A][T]; i[A][C]=j[A][C]; i[A][A]=j[A][A]; i[A][G]=j[A][G];
36  i[G][T]=j[G][T]; i[G][C]=j[G][C]; i[G][A]=j[G][A]; i[G][G]=j[G][G];
37 }
38 
39 //............................................................................
40 
41 static void ipol(double **i,double **j,double **k,double f) { double g; g=1.0-f;
42  i[T][T]=g*j[T][T]+f*k[T][T]; i[T][C]=g*j[T][C]+f*k[T][C]; i[T][A]=g*j[T][A]+f*k[T][A]; i[T][G]=g*j[T][G]+f*k[T][G];
43  i[C][T]=g*j[C][T]+f*k[C][T]; i[C][C]=g*j[C][C]+f*k[C][C]; i[C][A]=g*j[C][A]+f*k[C][A]; i[C][G]=g*j[C][G]+f*k[C][G];
44  i[A][T]=g*j[A][T]+f*k[A][T]; i[A][C]=g*j[A][C]+f*k[A][C]; i[A][A]=g*j[A][A]+f*k[A][A]; i[A][G]=g*j[A][G]+f*k[A][G];
45  i[G][T]=g*j[G][T]+f*k[G][T]; i[G][C]=g*j[G][C]+f*k[G][C]; i[G][A]=g*j[G][A]+f*k[G][A]; i[G][G]=g*j[G][G]+f*k[G][G];
46 }
47 
48 //............................................................................
49 
50 static void addmul(double **i,double **j,double f) {
51  i[T][T]+=j[T][T]*f; i[T][C]+=j[T][C]*f; i[T][A]+=j[T][A]*f; i[T][G]+=j[T][G]*f;
52  i[C][T]+=j[C][T]*f; i[C][C]+=j[C][C]*f; i[C][A]+=j[C][A]*f; i[C][G]+=j[C][G]*f;
53  i[A][T]+=j[A][T]*f; i[A][C]+=j[A][C]*f; i[A][A]+=j[A][A]*f; i[A][G]+=j[A][G]*f;
54  i[G][T]+=j[G][T]*f; i[G][C]+=j[G][C]*f; i[G][A]+=j[G][A]*f; i[G][G]+=j[G][G]*f;
55 }
56 
57 //............................................................................
58 
59 static void dot(double **i,double **j,double **k) {
60  i[T][T]=j[T][T]*k[T][T]+j[T][C]*k[C][T]+j[T][A]*k[A][T]+j[T][G]*k[G][T];
61  i[T][C]=j[T][T]*k[T][C]+j[T][C]*k[C][C]+j[T][A]*k[A][C]+j[T][G]*k[G][C];
62  i[T][A]=j[T][T]*k[T][A]+j[T][C]*k[C][A]+j[T][A]*k[A][A]+j[T][G]*k[G][A];
63  i[T][G]=j[T][T]*k[T][G]+j[T][C]*k[C][G]+j[T][A]*k[A][G]+j[T][G]*k[G][G];
64  i[C][T]=j[C][T]*k[T][T]+j[C][C]*k[C][T]+j[C][A]*k[A][T]+j[C][G]*k[G][T];
65  i[C][C]=j[C][T]*k[T][C]+j[C][C]*k[C][C]+j[C][A]*k[A][C]+j[C][G]*k[G][C];
66  i[C][A]=j[C][T]*k[T][A]+j[C][C]*k[C][A]+j[C][A]*k[A][A]+j[C][G]*k[G][A];
67  i[C][G]=j[C][T]*k[T][G]+j[C][C]*k[C][G]+j[C][A]*k[A][G]+j[C][G]*k[G][G];
68  i[A][T]=j[A][T]*k[T][T]+j[A][C]*k[C][T]+j[A][A]*k[A][T]+j[A][G]*k[G][T];
69  i[A][C]=j[A][T]*k[T][C]+j[A][C]*k[C][C]+j[A][A]*k[A][C]+j[A][G]*k[G][C];
70  i[A][A]=j[A][T]*k[T][A]+j[A][C]*k[C][A]+j[A][A]*k[A][A]+j[A][G]*k[G][A];
71  i[A][G]=j[A][T]*k[T][G]+j[A][C]*k[C][G]+j[A][A]*k[A][G]+j[A][G]*k[G][G];
72  i[G][T]=j[G][T]*k[T][T]+j[G][C]*k[C][T]+j[G][A]*k[A][T]+j[G][G]*k[G][T];
73  i[G][C]=j[G][T]*k[T][C]+j[G][C]*k[C][C]+j[G][A]*k[A][C]+j[G][G]*k[G][C];
74  i[G][A]=j[G][T]*k[T][A]+j[G][C]*k[C][A]+j[G][A]*k[A][A]+j[G][G]*k[G][A];
75  i[G][G]=j[G][T]*k[T][G]+j[G][C]*k[C][G]+j[G][A]*k[A][G]+j[G][G]*k[G][G];
76 }
77 
78 //============================================================================
79 
80 char encodeBase(char c) {
81 
82  switch(c) {
83  case 'U': return T;
84  case 'T': return T;
85  case 'C': return C;
86  case 'A': return A;
87  case 'G': return G;
88  default : return N;
89  }
90 
91 }
92 
93 //............................................................................
94 
95 char decodeBase(char c) {
96 
97  switch(c) {
98  case T: return 'T';
99  case C: return 'C';
100  case A: return 'A';
101  case G: return 'G';
102  case N: return '?';
103  default : return '-';
104  }
105 
106 }
107 
108 //============================================================================
109 
110 void normalizeBaseFreqs(
111  double *F,double fT,double fC,double fA,double fG
112  ) {
113  double s;
114 
115  s=fT+fC+fA+fG;
116 
117  if(s<REAL_MIN) {fT=1.; fC=1.; fA=1.; fG=1.; s=4.;}
118 
119  fT/=s; fC/=s; fA/=s; fG/=s;
120 
121  if(fT<ATOLPARAM) fT=ATOLPARAM;
122  if(fC<ATOLPARAM) fC=ATOLPARAM;
123  if(fA<ATOLPARAM) fA=ATOLPARAM;
124  if(fG<ATOLPARAM) fG=ATOLPARAM;
125 
126  F[T]=fT; F[C]=fC; F[A]=fA; F[G]=fG;
127 
128 }
129 
130 //............................................................................
131 
132 void normalizeRateParams(
133  double *X,double a,double b,double c,double d,double e,double f
134  ) { // TC TA TG CA CG AG
135  double s;
136 
137  s=a+b+c+d+e+f;
138 
139  if(s<REAL_MIN) {a=1.; b=1.; c=1.; d=1.; e=1.; s=6.;}
140 
141  a/=s; b/=s; c/=s; d/=s; e/=s;
142 
143  if(a<ATOLPARAM) a=ATOLPARAM;
144  if(b<ATOLPARAM) b=ATOLPARAM;
145  if(c<ATOLPARAM) c=ATOLPARAM;
146  if(d<ATOLPARAM) d=ATOLPARAM;
147  if(e<ATOLPARAM) e=ATOLPARAM;
148 
149  X[0]=a; X[1]=b; X[2]=c; X[3]=d; X[4]=e;
150 
151 }
152 
153 //============================================================================
154 
155 void initTrnsprob(double ****PP) {
156  double ***P; short j;
157 
158  P=(double***)newVector(128,sizeof(double **));
159  for(j=0;j<128;j++) P[j]=(double **)newMatrix(N,N,sizeof(double));
160 
161  *PP=P;
162 
163 }
164 
165 //............................................................................
166 
167 void uninitTrnsprob(double ****PP) {
168  double ***P; short j;
169 
170  P=*PP;
171 
172  for(j=0;j<128;j++) freeMatrix(&P[j]);
173  freeBlock(PP);
174 
175 }
176 
177 //============================================================================
178 
179 void getTrnsprob(double **P,double ***PP,double d) {
180  long I,J;
181  double *Y[N],YT[N],YC[N],YA[N],YG[N];
182  double *Z[N],ZT[N],ZC[N],ZA[N],ZG[N];
183  double *X[N],XT[N],XC[N],XA[N],XG[N];
184 
185  d/=EPS; I=(long)d;
186 
187  if( I < 0 ) {copy(P,PP[0]); return;}
188 
189  X[T]=XT; X[C]=XC; X[A]=XA; X[G]=XG;
190 
191  if( I < 32 ) {ipol(X,PP[I],PP[I==31?33:I+1],d-I); copy(P,X); return;}
192 
193  if( I < 1024 ) {
194  J=I%32; ipol(X,PP[J],PP[J==31?33:J+1],d-I);
195  dot(P,PP[ I/32 + 32 ],X);
196  return;
197  }
198 
199  Y[T]=YT; Y[C]=YC; Y[A]=YA; Y[G]=YG;
200 
201  if( I < 32768 ) {
202  J=I%32; ipol(X,PP[J],PP[J==31?33:J+1],d-I);
203  J=I%1024; dot(Y,PP[ J/32 + 32 ],X);
204  dot(P,PP[ I/1024 + 64 ],Y);
205  return;
206  }
207 
208  Z[T]=ZT; Z[C]=ZC; Z[A]=ZA; Z[G]=ZG;
209 
210  if( I < 1048576 ) {
211  J=I%32; ipol(X,PP[J],PP[J==31?33:J+1],d-I);
212  J=I%1024; dot(Y,PP[ J/32 + 32 ],X);
213  J=I%32768; dot(Z,PP[ J/1024 + 64 ],Y);
214  dot(P,PP[ I/32768 + 96 ],Z);
215  return;
216  }
217 
218  dot(Y,PP[ 31 + 32 ],PP[31]);
219  dot(Z,PP[ 31 + 64 ],Y);
220  dot(P,PP[ 31 + 96 ],Z);
221 
222 }
223 
224 //============================================================================
225 
226 void updateTrnsprob(double ***PP,double *F,double *X,short M) {
227  double s,a,b,c,d,e,f;
228  double *Q[N],QT[N],QC[N],QA[N],QG[N];
229  double *R[N],RT[N],RC[N],RA[N],RG[N];
230  double *S[N],ST[N],SC[N],SA[N],SG[N];
231 
232  Q[T]=QT; Q[C]=QC; Q[A]=QA; Q[G]=QG;
233  R[T]=RT; R[C]=RC; R[A]=RA; R[G]=RG;
234  S[T]=ST; S[C]=SC; S[A]=SA; S[G]=SG;
235 
236  switch(M) {
237  case REV: a = X[0]; b=X[1]; c=X[2]; d=X[3]; e=X[4]; f=1.0-(a+b+c+d+e); break;
238  case TN93: a = X[0]; b=c=d=e=X[1]; f=1.0-(a+b+c+d+e); break;
239  case HKY85: a = f=X[0]; b=c=d=e=0.25*(1.0-(a+f)); break;
240  default : ih_assert(0); break;
241  }
242 
243  s=0.5/(a*F[T]*F[C]+b*F[T]*F[A]+c*F[T]*F[G]+d*F[C]*F[A]+e*F[C]*F[G]+f*F[A]*F[G]);
244 
245  Q[T][C]=a*F[C]*s; Q[T][A]=b*F[A]*s; Q[T][G]=c*F[G]*s;
246  Q[C][T]=a*F[T]*s; Q[C][A]=d*F[A]*s; Q[C][G]=e*F[G]*s;
247  Q[A][T]=b*F[T]*s; Q[A][C]=d*F[C]*s; Q[A][G]=f*F[G]*s;
248  Q[G][T]=c*F[T]*s; Q[G][C]=e*F[C]*s; Q[G][A]=f*F[A]*s;
249 
250  Q[T][T] = - ( Q[T][C] + Q[T][A] + Q[T][G] );
251  Q[C][C] = - ( Q[C][T] + Q[C][A] + Q[C][G] );
252  Q[A][A] = - ( Q[A][T] + Q[A][C] + Q[A][G] );
253  Q[G][G] = - ( Q[G][T] + Q[G][C] + Q[G][A] );
254 
255  identity(PP[1]);
256  addmul(PP[1],Q,EPS);
257  dot(R,Q,Q);
258  addmul(PP[1],R,EPS*EPS/2.0);
259  dot(S,R,Q);
260  addmul(PP[1],S,EPS*EPS*EPS/6.0);
261  dot(S,R,R);
262  addmul(PP[1],S,EPS*EPS*EPS*EPS/24.0);
263 
264  identity(PP[ 0]);
265  identity(PP[32]);
266  identity(PP[64]);
267  identity(PP[96]);
268 
269  dot(PP[ 2 ],PP[ 1 ],PP[ 1 ]);
270  dot(PP[ 4 ],PP[ 2 ],PP[ 2 ]);
271  dot(PP[ 8 ],PP[ 4 ],PP[ 4 ]);
272  dot(PP[ 16 ],PP[ 8 ],PP[ 8 ]);
273  dot(PP[ 1 + 32 ],PP[ 16 ],PP[ 16 ]);
274  dot(PP[ 2 + 32 ],PP[ 1 + 32 ],PP[ 1 + 32 ]);
275  dot(PP[ 4 + 32 ],PP[ 2 + 32 ],PP[ 2 + 32 ]);
276  dot(PP[ 8 + 32 ],PP[ 4 + 32 ],PP[ 4 + 32 ]);
277  dot(PP[ 16 + 32 ],PP[ 8 + 32 ],PP[ 8 + 32 ]);
278  dot(PP[ 1 + 64 ],PP[ 16 + 32 ],PP[ 16 + 32 ]);
279  dot(PP[ 2 + 64 ],PP[ 1 + 64 ],PP[ 1 + 64 ]);
280  dot(PP[ 4 + 64 ],PP[ 2 + 64 ],PP[ 2 + 64 ]);
281  dot(PP[ 8 + 64 ],PP[ 4 + 64 ],PP[ 4 + 64 ]);
282  dot(PP[ 16 + 64 ],PP[ 8 + 64 ],PP[ 8 + 64 ]);
283  dot(PP[ 1 + 96 ],PP[ 16 + 64 ],PP[ 16 + 64 ]);
284  dot(PP[ 2 + 96 ],PP[ 1 + 96 ],PP[ 1 + 96 ]);
285  dot(PP[ 4 + 96 ],PP[ 2 + 96 ],PP[ 2 + 96 ]);
286  dot(PP[ 8 + 96 ],PP[ 4 + 96 ],PP[ 4 + 96 ]);
287  dot(PP[ 16 + 96 ],PP[ 8 + 96 ],PP[ 8 + 96 ]);
288 
289  dot(PP[ 3 ],PP[ 2 ],PP[ 1 ]);
290  dot(PP[ 6 ],PP[ 3 ],PP[ 3 ]);
291  dot(PP[ 12 ],PP[ 6 ],PP[ 6 ]);
292  dot(PP[ 24 ],PP[ 12 ],PP[ 12 ]);
293  dot(PP[ 3 + 32 ],PP[ 2 + 32 ],PP[ 1 + 32 ]);
294  dot(PP[ 6 + 32 ],PP[ 3 + 32 ],PP[ 3 + 32 ]);
295  dot(PP[ 12 + 32 ],PP[ 6 + 32 ],PP[ 6 + 32 ]);
296  dot(PP[ 24 + 32 ],PP[ 12 + 32 ],PP[ 12 + 32 ]);
297  dot(PP[ 3 + 64 ],PP[ 2 + 64 ],PP[ 1 + 64 ]);
298  dot(PP[ 6 + 64 ],PP[ 3 + 64 ],PP[ 3 + 64 ]);
299  dot(PP[ 12 + 64 ],PP[ 6 + 64 ],PP[ 6 + 64 ]);
300  dot(PP[ 24 + 64 ],PP[ 12 + 64 ],PP[ 12 + 64 ]);
301  dot(PP[ 3 + 96 ],PP[ 2 + 96 ],PP[ 1 + 96 ]);
302  dot(PP[ 6 + 96 ],PP[ 3 + 96 ],PP[ 3 + 96 ]);
303  dot(PP[ 12 + 96 ],PP[ 6 + 96 ],PP[ 6 + 96 ]);
304  dot(PP[ 24 + 96 ],PP[ 12 + 96 ],PP[ 12 + 96 ]);
305 
306  dot(PP[ 5 ],PP[ 3 ],PP[ 2 ]);
307  dot(PP[ 10 ],PP[ 5 ],PP[ 5 ]);
308  dot(PP[ 20 ],PP[ 10 ],PP[ 10 ]);
309  dot(PP[ 5 + 32 ],PP[ 3 + 32 ],PP[ 2 + 32 ]);
310  dot(PP[ 10 + 32 ],PP[ 5 + 32 ],PP[ 5 + 32 ]);
311  dot(PP[ 20 + 32 ],PP[ 10 + 32 ],PP[ 10 + 32 ]);
312  dot(PP[ 5 + 64 ],PP[ 3 + 64 ],PP[ 2 + 64 ]);
313  dot(PP[ 10 + 64 ],PP[ 5 + 64 ],PP[ 5 + 64 ]);
314  dot(PP[ 20 + 64 ],PP[ 10 + 64 ],PP[ 10 + 64 ]);
315  dot(PP[ 5 + 96 ],PP[ 3 + 96 ],PP[ 2 + 96 ]);
316  dot(PP[ 10 + 96 ],PP[ 5 + 96 ],PP[ 5 + 96 ]);
317  dot(PP[ 20 + 96 ],PP[ 10 + 96 ],PP[ 10 + 96 ]);
318 
319  dot(PP[ 7 ],PP[ 4 ],PP[ 3 ]);
320  dot(PP[ 14 ],PP[ 7 ],PP[ 7 ]);
321  dot(PP[ 28 ],PP[ 14 ],PP[ 14 ]);
322  dot(PP[ 7 + 32 ],PP[ 4 + 32 ],PP[ 3 + 32 ]);
323  dot(PP[ 14 + 32 ],PP[ 7 + 32 ],PP[ 7 + 32 ]);
324  dot(PP[ 28 + 32 ],PP[ 14 + 32 ],PP[ 14 + 32 ]);
325  dot(PP[ 7 + 64 ],PP[ 4 + 64 ],PP[ 3 + 64 ]);
326  dot(PP[ 14 + 64 ],PP[ 7 + 64 ],PP[ 7 + 64 ]);
327  dot(PP[ 28 + 64 ],PP[ 14 + 64 ],PP[ 14 + 64 ]);
328  dot(PP[ 7 + 96 ],PP[ 4 + 96 ],PP[ 3 + 96 ]);
329  dot(PP[ 14 + 96 ],PP[ 7 + 96 ],PP[ 7 + 96 ]);
330  dot(PP[ 28 + 96 ],PP[ 14 + 96 ],PP[ 14 + 96 ]);
331 
332  dot(PP[ 9 ],PP[ 5 ],PP[ 4 ]);
333  dot(PP[ 18 ],PP[ 9 ],PP[ 9 ]);
334  dot(PP[ 9 + 32 ],PP[ 5 + 32 ],PP[ 4 + 32 ]);
335  dot(PP[ 18 + 32 ],PP[ 9 + 32 ],PP[ 9 + 32 ]);
336  dot(PP[ 9 + 64 ],PP[ 5 + 64 ],PP[ 4 + 64 ]);
337  dot(PP[ 18 + 64 ],PP[ 9 + 64 ],PP[ 9 + 64 ]);
338  dot(PP[ 9 + 96 ],PP[ 5 + 96 ],PP[ 4 + 96 ]);
339  dot(PP[ 18 + 96 ],PP[ 9 + 96 ],PP[ 9 + 96 ]);
340 
341  dot(PP[ 11 ],PP[ 6 ],PP[ 5 ]);
342  dot(PP[ 22 ],PP[ 11 ],PP[ 11 ]);
343  dot(PP[ 11 + 32 ],PP[ 6 + 32 ],PP[ 5 + 32 ]);
344  dot(PP[ 22 + 32 ],PP[ 11 + 32 ],PP[ 11 + 32 ]);
345  dot(PP[ 11 + 64 ],PP[ 6 + 64 ],PP[ 5 + 64 ]);
346  dot(PP[ 22 + 64 ],PP[ 11 + 64 ],PP[ 11 + 64 ]);
347  dot(PP[ 11 + 96 ],PP[ 6 + 96 ],PP[ 5 + 96 ]);
348  dot(PP[ 22 + 96 ],PP[ 11 + 96 ],PP[ 11 + 96 ]);
349 
350  dot(PP[ 13 ],PP[ 7 ],PP[ 6 ]);
351  dot(PP[ 26 ],PP[ 13 ],PP[ 13 ]);
352  dot(PP[ 13 + 32 ],PP[ 7 + 32 ],PP[ 6 + 32 ]);
353  dot(PP[ 26 + 32 ],PP[ 13 + 32 ],PP[ 13 + 32 ]);
354  dot(PP[ 13 + 64 ],PP[ 7 + 64 ],PP[ 6 + 64 ]);
355  dot(PP[ 26 + 64 ],PP[ 13 + 64 ],PP[ 13 + 64 ]);
356  dot(PP[ 13 + 96 ],PP[ 7 + 96 ],PP[ 6 + 96 ]);
357  dot(PP[ 26 + 96 ],PP[ 13 + 96 ],PP[ 13 + 96 ]);
358 
359  dot(PP[ 15 ],PP[ 8 ],PP[ 7 ]);
360  dot(PP[ 30 ],PP[ 15 ],PP[ 15 ]);
361  dot(PP[ 15 + 32 ],PP[ 8 + 32 ],PP[ 7 + 32 ]);
362  dot(PP[ 30 + 32 ],PP[ 15 + 32 ],PP[ 15 + 32 ]);
363  dot(PP[ 15 + 64 ],PP[ 8 + 64 ],PP[ 7 + 64 ]);
364  dot(PP[ 30 + 64 ],PP[ 15 + 64 ],PP[ 15 + 64 ]);
365  dot(PP[ 15 + 96 ],PP[ 8 + 96 ],PP[ 7 + 96 ]);
366  dot(PP[ 30 + 96 ],PP[ 15 + 96 ],PP[ 15 + 96 ]);
367 
368  dot(PP[ 17 ],PP[ 9 ],PP[ 8 ]);
369  dot(PP[ 17 + 32 ],PP[ 9 + 32 ],PP[ 8 + 32 ]);
370  dot(PP[ 17 + 64 ],PP[ 9 + 64 ],PP[ 8 + 64 ]);
371  dot(PP[ 17 + 96 ],PP[ 9 + 96 ],PP[ 8 + 96 ]);
372 
373  dot(PP[ 19 ],PP[ 10 ],PP[ 9 ]);
374  dot(PP[ 19 + 32 ],PP[ 10 + 32 ],PP[ 9 + 32 ]);
375  dot(PP[ 19 + 64 ],PP[ 10 + 64 ],PP[ 9 + 64 ]);
376  dot(PP[ 19 + 96 ],PP[ 10 + 96 ],PP[ 9 + 96 ]);
377 
378  dot(PP[ 21 ],PP[ 11 ],PP[ 10 ]);
379  dot(PP[ 21 + 32 ],PP[ 11 + 32 ],PP[ 10 + 32 ]);
380  dot(PP[ 21 + 64 ],PP[ 11 + 64 ],PP[ 10 + 64 ]);
381  dot(PP[ 21 + 96 ],PP[ 11 + 96 ],PP[ 10 + 96 ]);
382 
383  dot(PP[ 23 ],PP[ 12 ],PP[ 11 ]);
384  dot(PP[ 23 + 32 ],PP[ 12 + 32 ],PP[ 11 + 32 ]);
385  dot(PP[ 23 + 64 ],PP[ 12 + 64 ],PP[ 11 + 64 ]);
386  dot(PP[ 23 + 96 ],PP[ 12 + 96 ],PP[ 11 + 96 ]);
387 
388  dot(PP[ 25 ],PP[ 13 ],PP[ 12 ]);
389  dot(PP[ 25 + 32 ],PP[ 13 + 32 ],PP[ 12 + 32 ]);
390  dot(PP[ 25 + 64 ],PP[ 13 + 64 ],PP[ 12 + 64 ]);
391  dot(PP[ 25 + 96 ],PP[ 13 + 96 ],PP[ 12 + 96 ]);
392 
393  dot(PP[ 27 ],PP[ 14 ],PP[ 13 ]);
394  dot(PP[ 27 + 32 ],PP[ 14 + 32 ],PP[ 13 + 32 ]);
395  dot(PP[ 27 + 64 ],PP[ 14 + 64 ],PP[ 13 + 64 ]);
396  dot(PP[ 27 + 96 ],PP[ 14 + 96 ],PP[ 13 + 96 ]);
397 
398  dot(PP[ 29 ],PP[ 15 ],PP[ 14 ]);
399  dot(PP[ 29 + 32 ],PP[ 15 + 32 ],PP[ 14 + 32 ]);
400  dot(PP[ 29 + 64 ],PP[ 15 + 64 ],PP[ 14 + 64 ]);
401  dot(PP[ 29 + 96 ],PP[ 15 + 96 ],PP[ 14 + 96 ]);
402 
403  dot(PP[ 31 ],PP[ 16 ],PP[ 15 ]);
404  dot(PP[ 31 + 32 ],PP[ 16 + 32 ],PP[ 15 + 32 ]);
405  dot(PP[ 31 + 64 ],PP[ 16 + 64 ],PP[ 15 + 64 ]);
406  dot(PP[ 31 + 96 ],PP[ 16 + 96 ],PP[ 15 + 96 ]);
407 
408 }
409 
initTrnsprob
void initTrnsprob(double ****PP)
Definition: trnsprob.cxx:155
I
static int I
Definition: align.cxx:489
encodeBase
char encodeBase(char c)
Definition: trnsprob.cxx:80
defs.h
REV
Definition: trnsprob.h:21
uninitTrnsprob
void uninitTrnsprob(double ****PP)
Definition: trnsprob.cxx:167
dot
static void dot(double **i, double **j, double **k)
Definition: trnsprob.cxx:59
normalizeBaseFreqs
void normalizeBaseFreqs(double *F, double fT, double fC, double fA, double fG)
Definition: trnsprob.cxx:110
A
Definition: trnsprob.h:20
long
long
Definition: AW_awar.cxx:152
ih_assert
#define ih_assert(bed)
Definition: trnsprob.cxx:18
N
Definition: trnsprob.h:20
EPS
#define EPS
Definition: trnsprob.cxx:19
updateTrnsprob
void updateTrnsprob(double ***PP, double *F, double *X, short M)
Definition: trnsprob.cxx:226
decodeBase
char decodeBase(char c)
Definition: trnsprob.cxx:95
Z
static Island * Z
Definition: align.cxx:71
REAL_MIN
#define REAL_MIN
Definition: defs.h:26
identity
static void identity(double **i)
Definition: trnsprob.cxx:23
newVector
#define newVector(n, s)
Definition: mem.h:23
HKY85
Definition: trnsprob.h:21
M
static void * M
Definition: mem.cxx:18
ipol
static void ipol(double **i, double **j, double **k, double f)
Definition: trnsprob.cxx:41
J
static int J
Definition: align.cxx:489
copy
static void copy(double **i, double **j)
Definition: trnsprob.cxx:32
C
Definition: trnsprob.h:20
freeBlock
#define freeBlock(v)
Definition: mem.h:25
addmul
static void addmul(double **i, double **j, double f)
Definition: trnsprob.cxx:50
ATOLPARAM
#define ATOLPARAM
Definition: trnsprob.h:17
newMatrix
void ** newMatrix(size_t nrow, size_t ncol, size_t s)
Definition: mem.cxx:87
getTrnsprob
void getTrnsprob(double **P, double ***PP, double d)
Definition: trnsprob.cxx:179
arb_simple_assert.h
T
Definition: trnsprob.h:20
normalizeRateParams
void normalizeRateParams(double *X, double a, double b, double c, double d, double e, double f)
Definition: trnsprob.cxx:132
G
Definition: trnsprob.h:20
TN93
Definition: trnsprob.h:21
P
static Params P
Definition: arb_probe.cxx:81
freeMatrix
#define freeMatrix(m)
Definition: mem.h:40
s
GB_write_int const char s
Definition: AW_awar.cxx:154
Generated by   doxygen 1.8.8