Context Navigation

source: branches/port5/GDE/RAxML/models.c

Visit:

Last change on this file was 5262, checked in by westram, 17 years ago
update to RAxML 7.0.3
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 86.8 KB

Line
1	/* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2	* Copyright August 2006 by Alexandros Stamatakis
3	*
4	* Partially derived from
5	* fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6	*
7	* and
8	*
9	* Programs of the PHYLIP package by Joe Felsenstein.
10	*
11	* This program is free software; you may redistribute it and/or modify its
12	* under the terms of the GNU General Public License as published by the Free
13	* Software Foundation; either version 2 of the License, or (at your option)
14	* any later version.
15	*
16	* This program is distributed in the hope that it will be useful, but
17	* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18	* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19	* for more details.
20	*
21	*
22	* For any other enquiries send an Email to Alexandros Stamatakis
23	* Alexandros.Stamatakis@epfl.ch
24	*
25	* When publishing work that is based on the results from RAxML-VI-HPC please cite:
26	*
27	* Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands
28	* of taxa and mixed models".
29	* Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30	*/
31
32	#ifndef WIN32
33	#include <sys/times.h>
34	#include <sys/types.h>
35	#include <sys/time.h>
36	#include <unistd.h>
37	#endif
38
39	#include <math.h>
40	#include <time.h>
41	#include <stdlib.h>
42	#include <stdio.h>
43	#include <ctype.h>
44	#include <string.h>
45
46
47	#include "axml.h"
48
49	extern int optimizeRatesInvocations;
50	extern int optimizeRateCategoryInvocations;
51	extern int optimizeAlphaInvocations;
52	extern int optimizeTTRatioInvocations;
53	extern int optimizeInvarInvocations;
54
55	extern int NumberOfThreads;
56
57
58
59	void makeVal(char code, double *val)
60	{
61	double i_value;
62	int j;
63
64	assert(code >= 0 && code <= 22);
65
66	if(code < 22)
67	i_value = 0.0;
68	else
69	i_value = 1.0;
70
71	for(j = 0; j < 20; j++)
72	val[j] = i_value;
73
74	if(code < 20)
75	val[code] = 1.0;
76	else
77	{
78	if(code == 20)
79	val[2] = val[3] = 1.0;
80	if(code == 21)
81	val[5] = val[6] = 1.0;
82	}
83	}
84
85
86	void baseFrequenciesGTR(rawdata rdta, cruncheddata cdta, tree tr, analdef adef)
87	{
88	double sum, suma, sumc, sumg, sumt, wj, fa, fc, fg, ft, freqa, freqc, freqg, freqt, pfreqs[20], sumf[20], val[20], temp[20], acc;
89	int i, j, k, l, code;
90	char *yptr;
91	int model;
92	int lower;
93	int upper;
94
95	if(!adef->useMultipleModel)
96	{
97	assert(tr->partitionData[0].lower == 0);
98	assert(tr->partitionData[0].upper == tr->cdta->endsite);
99	}
100
101
102	for(model = 0; model < tr->NumberOfModels; model++)
103	{
104	/*lower = tr->modelIndices[model][0];
105	upper = tr->modelIndices[model][1];*/
106	lower = tr->partitionData[model].lower;
107	upper = tr->partitionData[model].upper;
108
109	switch(/tr->dataType[model]/tr->partitionData[model].dataType)
110	{
111	case AA_DATA:
112	for(l = 0; l < 20; l++)
113	pfreqs[l] = 0.05;
114
115	for (k = 1; k <= 8; k++)
116	{
117	for(l = 0; l < 20; l++)
118	sumf[l] = 0.0;
119
120	for (i = 0; i < rdta->numsp; i++)
121	{
122	yptr = &(rdta->y0[i * tr->originalCrunchedLength]);
123	for(j = lower; j < upper; j++)
124	{
125	makeVal(yptr[j], val);
126	for(l = 0; l < 20; l++)
127	{
128	if(val[l] != 0.0)
129	temp[l] = pfreqs[l] * val[l];
130	else
131	temp[l] = 0.0;
132	}
133
134	acc = 0;
135	for(l = 0; l < 20; l++)
136	{
137	if(temp[l] != 0.0)
138	acc += temp[l];
139	}
140
141	wj = ((double)cdta->aliaswgt[j]) / acc;
142
143	for(l = 0; l < 20; l++)
144	{
145	if(temp[l] != 0.0)
146	{
147	sumf[l] += wj * temp[l];
148	}
149	}
150	}
151	}
152
153	acc = 0;
154	for(l = 0; l < 20; l++)
155	{
156	if(sumf[l] != 0.0)
157	acc += sumf[l];
158	}
159
160	for(l = 0; l < 20; l++)
161	pfreqs[l] = sumf[l] / acc;
162	}
163
164	{
165	int countZeros = 0;
166
167	for(l = 0; l < 20; l++)
168	{
169	if(pfreqs[l] == 0.0)
170	countZeros++;
171	}
172
173	if(countZeros > 0)
174	{
175	double correction;
176	correction = (FREQ_MIN * (double)countZeros)/ ((double)(20 - countZeros));
177
178	for(l = 0; l < 20; l++)
179	{
180	if(pfreqs[l] == 0.0)
181	pfreqs[l] = FREQ_MIN;
182	else
183	pfreqs[l] -= correction;
184
185	tr->frequencies_AA[model * 20 + l] = pfreqs[l];
186	}
187	}
188	else
189	{
190	for(l = 0; l < 20; l++)
191	tr->frequencies_AA[model * 20 + l] = pfreqs[l];
192	}
193
194	}
195	break;
196	case DNA_DATA:
197	freqa = 0.25;
198	freqc = 0.25;
199	freqg = 0.25;
200	freqt = 0.25;
201	for (k = 1; k <= 8; k++)
202	{
203	suma = 0.0;
204	sumc = 0.0;
205	sumg = 0.0;
206	sumt = 0.0;
207	for (i = 0; i < rdta->numsp; i++)
208	{
209	yptr = &(rdta->y0[i * tr->originalCrunchedLength + lower]);
210	for (j = lower; j < upper; j++)
211	{
212	code = *yptr++;
213	if(code > 0)
214	{
215	fa = freqa * ( code & 1);
216	fc = freqc * ((code >> 1) & 1);
217	fg = freqg * ((code >> 2) & 1);
218	ft = freqt * ((code >> 3) & 1);
219	wj = cdta->aliaswgt[j] / (fa + fc + fg + ft);
220	suma += wj * fa;
221	sumc += wj * fc;
222	sumg += wj * fg;
223	sumt += wj * ft;
224	}
225	}
226	}
227	sum = suma + sumc + sumg + sumt;
228	freqa = suma / sum;
229	freqc = sumc / sum;
230	freqg = sumg / sum;
231	freqt = sumt / sum;
232	}
233
234	if(freqa == 0.0 \|\| freqc == 0.0 \|\| freqg == 0.0 \|\| freqt == 0.0)
235	{
236	printf("ERROR: one of the DNA base freqs is equal to zero\n");
237	exit(-1);
238	}
239
240	tr->frequencies_DNA[model * 4] = freqa;
241	tr->frequencies_DNA[model * 4 + 1] = freqc;
242	tr->frequencies_DNA[model * 4 + 2] = freqg;
243	tr->frequencies_DNA[model * 4 + 3] = freqt;
244	break;
245	default:
246	assert(0);
247	}
248	}
249
250	return;
251	}
252
253
254
255	static void initProtMat(int model, double f[20], int proteinMatrix, double *ext_initialRates,
256	boolean userProteinModel, double *externalAAMatrix)
257	{
258	double q[20][20];
259	double daa[400], max, temp;
260	int i, j, r;
261	double initialRates = &(ext_initialRates[model 190]);
262	double scaler;
263
264	if(userProteinModel)
265	{
266	assert(externalAAMatrix);
267	memcpy(daa, externalAAMatrix, 400 * sizeof(double));
268	memcpy(f, &(externalAAMatrix[400]), 20 * sizeof(double));
269	}
270	else
271	{
272	switch(proteinMatrix)
273	{
274	case DAYHOFF:
275	{
276	daa[ 120+ 0] = 27.00; daa[ 220+ 0] = 98.00; daa[ 220+ 1] = 32.00; daa[ 320+ 0] = 120.00;
277	daa[ 320+ 1] = 0.00; daa[ 320+ 2] = 905.00; daa[ 420+ 0] = 36.00; daa[ 420+ 1] = 23.00;
278	daa[ 420+ 2] = 0.00; daa[ 420+ 3] = 0.00; daa[ 520+ 0] = 89.00; daa[ 520+ 1] = 246.00;
279	daa[ 520+ 2] = 103.00; daa[ 520+ 3] = 134.00; daa[ 520+ 4] = 0.00; daa[ 620+ 0] = 198.00;
280	daa[ 620+ 1] = 1.00; daa[ 620+ 2] = 148.00; daa[ 620+ 3] = 1153.00; daa[ 620+ 4] = 0.00;
281	daa[ 620+ 5] = 716.00; daa[ 720+ 0] = 240.00; daa[ 720+ 1] = 9.00; daa[ 720+ 2] = 139.00;
282	daa[ 720+ 3] = 125.00; daa[ 720+ 4] = 11.00; daa[ 720+ 5] = 28.00; daa[ 720+ 6] = 81.00;
283	daa[ 820+ 0] = 23.00; daa[ 820+ 1] = 240.00; daa[ 820+ 2] = 535.00; daa[ 820+ 3] = 86.00;
284	daa[ 820+ 4] = 28.00; daa[ 820+ 5] = 606.00; daa[ 820+ 6] = 43.00; daa[ 820+ 7] = 10.00;
285	daa[ 920+ 0] = 65.00; daa[ 920+ 1] = 64.00; daa[ 920+ 2] = 77.00; daa[ 920+ 3] = 24.00;
286	daa[ 920+ 4] = 44.00; daa[ 920+ 5] = 18.00; daa[ 920+ 6] = 61.00; daa[ 920+ 7] = 0.00;
287	daa[ 920+ 8] = 7.00; daa[1020+ 0] = 41.00; daa[1020+ 1] = 15.00; daa[1020+ 2] = 34.00;
288	daa[1020+ 3] = 0.00; daa[1020+ 4] = 0.00; daa[1020+ 5] = 73.00; daa[1020+ 6] = 11.00;
289	daa[1020+ 7] = 7.00; daa[1020+ 8] = 44.00; daa[1020+ 9] = 257.00; daa[1120+ 0] = 26.00;
290	daa[1120+ 1] = 464.00; daa[1120+ 2] = 318.00; daa[1120+ 3] = 71.00; daa[1120+ 4] = 0.00;
291	daa[1120+ 5] = 153.00; daa[1120+ 6] = 83.00; daa[1120+ 7] = 27.00; daa[1120+ 8] = 26.00;
292	daa[1120+ 9] = 46.00; daa[1120+10] = 18.00; daa[1220+ 0] = 72.00; daa[1220+ 1] = 90.00;
293	daa[1220+ 2] = 1.00; daa[1220+ 3] = 0.00; daa[1220+ 4] = 0.00; daa[1220+ 5] = 114.00;
294	daa[1220+ 6] = 30.00; daa[1220+ 7] = 17.00; daa[1220+ 8] = 0.00; daa[1220+ 9] = 336.00;
295	daa[1220+10] = 527.00; daa[1220+11] = 243.00; daa[1320+ 0] = 18.00; daa[1320+ 1] = 14.00;
296	daa[1320+ 2] = 14.00; daa[1320+ 3] = 0.00; daa[1320+ 4] = 0.00; daa[1320+ 5] = 0.00;
297	daa[1320+ 6] = 0.00; daa[1320+ 7] = 15.00; daa[1320+ 8] = 48.00; daa[1320+ 9] = 196.00;
298	daa[1320+10] = 157.00; daa[1320+11] = 0.00; daa[1320+12] = 92.00; daa[1420+ 0] = 250.00;
299	daa[1420+ 1] = 103.00; daa[1420+ 2] = 42.00; daa[1420+ 3] = 13.00; daa[1420+ 4] = 19.00;
300	daa[1420+ 5] = 153.00; daa[1420+ 6] = 51.00; daa[1420+ 7] = 34.00; daa[1420+ 8] = 94.00;
301	daa[1420+ 9] = 12.00; daa[1420+10] = 32.00; daa[1420+11] = 33.00; daa[1420+12] = 17.00;
302	daa[1420+13] = 11.00; daa[1520+ 0] = 409.00; daa[1520+ 1] = 154.00; daa[1520+ 2] = 495.00;
303	daa[1520+ 3] = 95.00; daa[1520+ 4] = 161.00; daa[1520+ 5] = 56.00; daa[1520+ 6] = 79.00;
304	daa[1520+ 7] = 234.00; daa[1520+ 8] = 35.00; daa[1520+ 9] = 24.00; daa[1520+10] = 17.00;
305	daa[1520+11] = 96.00; daa[1520+12] = 62.00; daa[1520+13] = 46.00; daa[1520+14] = 245.00;
306	daa[1620+ 0] = 371.00; daa[1620+ 1] = 26.00; daa[1620+ 2] = 229.00; daa[1620+ 3] = 66.00;
307	daa[1620+ 4] = 16.00; daa[1620+ 5] = 53.00; daa[1620+ 6] = 34.00; daa[1620+ 7] = 30.00;
308	daa[1620+ 8] = 22.00; daa[1620+ 9] = 192.00; daa[1620+10] = 33.00; daa[1620+11] = 136.00;
309	daa[1620+12] = 104.00; daa[1620+13] = 13.00; daa[1620+14] = 78.00; daa[1620+15] = 550.00;
310	daa[1720+ 0] = 0.00; daa[1720+ 1] = 201.00; daa[1720+ 2] = 23.00; daa[1720+ 3] = 0.00;
311	daa[1720+ 4] = 0.00; daa[1720+ 5] = 0.00; daa[1720+ 6] = 0.00; daa[1720+ 7] = 0.00;
312	daa[1720+ 8] = 27.00; daa[1720+ 9] = 0.00; daa[1720+10] = 46.00; daa[1720+11] = 0.00;
313	daa[1720+12] = 0.00; daa[1720+13] = 76.00; daa[1720+14] = 0.00; daa[1720+15] = 75.00;
314	daa[1720+16] = 0.00; daa[1820+ 0] = 24.00; daa[1820+ 1] = 8.00; daa[1820+ 2] = 95.00;
315	daa[1820+ 3] = 0.00; daa[1820+ 4] = 96.00; daa[1820+ 5] = 0.00; daa[1820+ 6] = 22.00;
316	daa[1820+ 7] = 0.00; daa[1820+ 8] = 127.00; daa[1820+ 9] = 37.00; daa[1820+10] = 28.00;
317	daa[1820+11] = 13.00; daa[1820+12] = 0.00; daa[1820+13] = 698.00; daa[1820+14] = 0.00;
318	daa[1820+15] = 34.00; daa[1820+16] = 42.00; daa[1820+17] = 61.00; daa[1920+ 0] = 208.00;
319	daa[1920+ 1] = 24.00; daa[1920+ 2] = 15.00; daa[1920+ 3] = 18.00; daa[1920+ 4] = 49.00;
320	daa[1920+ 5] = 35.00; daa[1920+ 6] = 37.00; daa[1920+ 7] = 54.00; daa[1920+ 8] = 44.00;
321	daa[1920+ 9] = 889.00; daa[1920+10] = 175.00; daa[1920+11] = 10.00; daa[1920+12] = 258.00;
322	daa[1920+13] = 12.00; daa[1920+14] = 48.00; daa[1920+15] = 30.00; daa[1920+16] = 157.00;
323	daa[1920+17] = 0.00; daa[1920+18] = 28.00;
324
325
326	f[ 0] = 0.087000; f[ 1] = 0.041000; f[ 2] = 0.040000; f[ 3] = 0.047000;
327	f[ 4] = 0.034000; f[ 5] = 0.038000; f[ 6] = 0.050000; f[ 7] = 0.089000;
328	f[ 8] = 0.034000; f[ 9] = 0.037000; f[10] = 0.085000; f[11] = 0.080000;
329	f[12] = 0.014000; f[13] = 0.040000; f[14] = 0.051000; f[15] = 0.070000;
330	f[16] = 0.058000; f[17] = 0.011000; f[18] = 0.030000; f[19] = 0.064000;
331	}
332	break;
333	case DCMUT:
334	{
335	daa[ 120+ 0] = 26.78280; daa[ 220+ 0] = 98.44740; daa[ 220+ 1] = 32.70590; daa[ 320+ 0] = 119.98050;
336	daa[ 320+ 1] = 0.00000; daa[ 320+ 2] = 893.15150; daa[ 420+ 0] = 36.00160; daa[ 420+ 1] = 23.23740;
337	daa[ 420+ 2] = 0.00000; daa[ 420+ 3] = 0.00000; daa[ 520+ 0] = 88.77530; daa[ 520+ 1] = 243.99390;
338	daa[ 520+ 2] = 102.85090; daa[ 520+ 3] = 134.85510; daa[ 520+ 4] = 0.00000; daa[ 620+ 0] = 196.11670;
339	daa[ 620+ 1] = 0.00000; daa[ 620+ 2] = 149.34090; daa[ 620+ 3] = 1138.86590; daa[ 620+ 4] = 0.00000;
340	daa[ 620+ 5] = 708.60220; daa[ 720+ 0] = 238.61110; daa[ 720+ 1] = 8.77910; daa[ 720+ 2] = 138.53520;
341	daa[ 720+ 3] = 124.09810; daa[ 720+ 4] = 10.72780; daa[ 720+ 5] = 28.15810; daa[ 720+ 6] = 81.19070;
342	daa[ 820+ 0] = 22.81160; daa[ 820+ 1] = 238.31480; daa[ 820+ 2] = 529.00240; daa[ 820+ 3] = 86.82410;
343	daa[ 820+ 4] = 28.27290; daa[ 820+ 5] = 601.16130; daa[ 820+ 6] = 43.94690; daa[ 820+ 7] = 10.68020;
344	daa[ 920+ 0] = 65.34160; daa[ 920+ 1] = 63.26290; daa[ 920+ 2] = 76.80240; daa[ 920+ 3] = 23.92480;
345	daa[ 920+ 4] = 43.80740; daa[ 920+ 5] = 18.03930; daa[ 920+ 6] = 60.95260; daa[ 920+ 7] = 0.00000;
346	daa[ 920+ 8] = 7.69810; daa[1020+ 0] = 40.64310; daa[1020+ 1] = 15.49240; daa[1020+ 2] = 34.11130;
347	daa[1020+ 3] = 0.00000; daa[1020+ 4] = 0.00000; daa[1020+ 5] = 73.07720; daa[1020+ 6] = 11.28800;
348	daa[1020+ 7] = 7.15140; daa[1020+ 8] = 44.35040; daa[1020+ 9] = 255.66850; daa[1120+ 0] = 25.86350;
349	daa[1120+ 1] = 461.01240; daa[1120+ 2] = 314.83710; daa[1120+ 3] = 71.69130; daa[1120+ 4] = 0.00000;
350	daa[1120+ 5] = 151.90780; daa[1120+ 6] = 83.00780; daa[1120+ 7] = 26.76830; daa[1120+ 8] = 27.04750;
351	daa[1120+ 9] = 46.08570; daa[1120+10] = 18.06290; daa[1220+ 0] = 71.78400; daa[1220+ 1] = 89.63210;
352	daa[1220+ 2] = 0.00000; daa[1220+ 3] = 0.00000; daa[1220+ 4] = 0.00000; daa[1220+ 5] = 112.74990;
353	daa[1220+ 6] = 30.48030; daa[1220+ 7] = 17.03720; daa[1220+ 8] = 0.00000; daa[1220+ 9] = 333.27320;
354	daa[1220+10] = 523.01150; daa[1220+11] = 241.17390; daa[1320+ 0] = 18.36410; daa[1320+ 1] = 13.69060;
355	daa[1320+ 2] = 13.85030; daa[1320+ 3] = 0.00000; daa[1320+ 4] = 0.00000; daa[1320+ 5] = 0.00000;
356	daa[1320+ 6] = 0.00000; daa[1320+ 7] = 15.34780; daa[1320+ 8] = 47.59270; daa[1320+ 9] = 195.19510;
357	daa[1320+10] = 156.51600; daa[1320+11] = 0.00000; daa[1320+12] = 92.18600; daa[1420+ 0] = 248.59200;
358	daa[1420+ 1] = 102.83130; daa[1420+ 2] = 41.92440; daa[1420+ 3] = 13.39400; daa[1420+ 4] = 18.75500;
359	daa[1420+ 5] = 152.61880; daa[1420+ 6] = 50.70030; daa[1420+ 7] = 34.71530; daa[1420+ 8] = 93.37090;
360	daa[1420+ 9] = 11.91520; daa[1420+10] = 31.62580; daa[1420+11] = 33.54190; daa[1420+12] = 17.02050;
361	daa[1420+13] = 11.05060; daa[1520+ 0] = 405.18700; daa[1520+ 1] = 153.15900; daa[1520+ 2] = 488.58920;
362	daa[1520+ 3] = 95.60970; daa[1520+ 4] = 159.83560; daa[1520+ 5] = 56.18280; daa[1520+ 6] = 79.39990;
363	daa[1520+ 7] = 232.22430; daa[1520+ 8] = 35.36430; daa[1520+ 9] = 24.79550; daa[1520+10] = 17.14320;
364	daa[1520+11] = 95.45570; daa[1520+12] = 61.99510; daa[1520+13] = 45.99010; daa[1520+14] = 242.72020;
365	daa[1620+ 0] = 368.03650; daa[1620+ 1] = 26.57450; daa[1620+ 2] = 227.16970; daa[1620+ 3] = 66.09300;
366	daa[1620+ 4] = 16.23660; daa[1620+ 5] = 52.56510; daa[1620+ 6] = 34.01560; daa[1620+ 7] = 30.66620;
367	daa[1620+ 8] = 22.63330; daa[1620+ 9] = 190.07390; daa[1620+10] = 33.10900; daa[1620+11] = 135.05990;
368	daa[1620+12] = 103.15340; daa[1620+13] = 13.66550; daa[1620+14] = 78.28570; daa[1620+15] = 543.66740;
369	daa[1720+ 0] = 0.00000; daa[1720+ 1] = 200.13750; daa[1720+ 2] = 22.49680; daa[1720+ 3] = 0.00000;
370	daa[1720+ 4] = 0.00000; daa[1720+ 5] = 0.00000; daa[1720+ 6] = 0.00000; daa[1720+ 7] = 0.00000;
371	daa[1720+ 8] = 27.05640; daa[1720+ 9] = 0.00000; daa[1720+10] = 46.17760; daa[1720+11] = 0.00000;
372	daa[1720+12] = 0.00000; daa[1720+13] = 76.23540; daa[1720+14] = 0.00000; daa[1720+15] = 74.08190;
373	daa[1720+16] = 0.00000; daa[1820+ 0] = 24.41390; daa[1820+ 1] = 7.80120; daa[1820+ 2] = 94.69400;
374	daa[1820+ 3] = 0.00000; daa[1820+ 4] = 95.31640; daa[1820+ 5] = 0.00000; daa[1820+ 6] = 21.47170;
375	daa[1820+ 7] = 0.00000; daa[1820+ 8] = 126.54000; daa[1820+ 9] = 37.48340; daa[1820+10] = 28.65720;
376	daa[1820+11] = 13.21420; daa[1820+12] = 0.00000; daa[1820+13] = 695.26290; daa[1820+14] = 0.00000;
377	daa[1820+15] = 33.62890; daa[1820+16] = 41.78390; daa[1820+17] = 60.80700; daa[1920+ 0] = 205.95640;
378	daa[1920+ 1] = 24.03680; daa[1920+ 2] = 15.80670; daa[1920+ 3] = 17.83160; daa[1920+ 4] = 48.46780;
379	daa[1920+ 5] = 34.69830; daa[1920+ 6] = 36.72500; daa[1920+ 7] = 53.81650; daa[1920+ 8] = 43.87150;
380	daa[1920+ 9] = 881.00380; daa[1920+10] = 174.51560; daa[1920+11] = 10.38500; daa[1920+12] = 256.59550;
381	daa[1920+13] = 12.36060; daa[1920+14] = 48.50260; daa[1920+15] = 30.38360; daa[1920+16] = 156.19970;
382	daa[1920+17] = 0.00000; daa[1920+18] = 27.93790;
383
384	/*
385	f[ 0] = 0.087127; f[ 1] = 0.040904; f[ 2] = 0.040432; f[ 3] = 0.046872;
386	f[ 4] = 0.033474; f[ 5] = 0.038255; f[ 6] = 0.049530; f[ 7] = 0.088612;
387	f[ 8] = 0.033619; f[ 9] = 0.036886; f[10] = 0.085357; f[11] = 0.080481;
388	f[12] = 0.014753; f[13] = 0.039772; f[14] = 0.050680; f[15] = 0.069577;
389	f[16] = 0.058542; f[17] = 0.010494; f[18] = 0.029916; f[19] = 0.064718;
390	*/
391
392	f[ 0] = 0.08700; f[ 1] = 0.04100; f[ 2] = 0.04000; f[ 3] = 0.04700;
393	f[ 4] = 0.03300; f[ 5] = 0.03800; f[ 6] = 0.04900; f[ 7] = 0.08900;
394	f[ 8] = 0.03400; f[ 9] = 0.03700; f[10] = 0.08500; f[11] = 0.08000;
395	f[12] = 0.01500; f[13] = 0.04000; f[14] = 0.05200; f[15] = 0.06900;
396	f[16] = 0.05900; f[17] = 0.01000; f[18] = 0.03000; f[19] = 0.06500;
397
398	}
399	break;
400	case JTT:
401	{
402	daa[ 120+ 0] = 58.00; daa[ 220+ 0] = 54.00; daa[ 220+ 1] = 45.00; daa[ 320+ 0] = 81.00;
403	daa[ 320+ 1] = 16.00; daa[ 320+ 2] = 528.00; daa[ 420+ 0] = 56.00; daa[ 420+ 1] = 113.00;
404	daa[ 420+ 2] = 34.00; daa[ 420+ 3] = 10.00; daa[ 520+ 0] = 57.00; daa[ 520+ 1] = 310.00;
405	daa[ 520+ 2] = 86.00; daa[ 520+ 3] = 49.00; daa[ 520+ 4] = 9.00; daa[ 620+ 0] = 105.00;
406	daa[ 620+ 1] = 29.00; daa[ 620+ 2] = 58.00; daa[ 620+ 3] = 767.00; daa[ 620+ 4] = 5.00;
407	daa[ 620+ 5] = 323.00; daa[ 720+ 0] = 179.00; daa[ 720+ 1] = 137.00; daa[ 720+ 2] = 81.00;
408	daa[ 720+ 3] = 130.00; daa[ 720+ 4] = 59.00; daa[ 720+ 5] = 26.00; daa[ 720+ 6] = 119.00;
409	daa[ 820+ 0] = 27.00; daa[ 820+ 1] = 328.00; daa[ 820+ 2] = 391.00; daa[ 820+ 3] = 112.00;
410	daa[ 820+ 4] = 69.00; daa[ 820+ 5] = 597.00; daa[ 820+ 6] = 26.00; daa[ 820+ 7] = 23.00;
411	daa[ 920+ 0] = 36.00; daa[ 920+ 1] = 22.00; daa[ 920+ 2] = 47.00; daa[ 920+ 3] = 11.00;
412	daa[ 920+ 4] = 17.00; daa[ 920+ 5] = 9.00; daa[ 920+ 6] = 12.00; daa[ 920+ 7] = 6.00;
413	daa[ 920+ 8] = 16.00; daa[1020+ 0] = 30.00; daa[1020+ 1] = 38.00; daa[1020+ 2] = 12.00;
414	daa[1020+ 3] = 7.00; daa[1020+ 4] = 23.00; daa[1020+ 5] = 72.00; daa[1020+ 6] = 9.00;
415	daa[1020+ 7] = 6.00; daa[1020+ 8] = 56.00; daa[1020+ 9] = 229.00; daa[1120+ 0] = 35.00;
416	daa[1120+ 1] = 646.00; daa[1120+ 2] = 263.00; daa[1120+ 3] = 26.00; daa[1120+ 4] = 7.00;
417	daa[1120+ 5] = 292.00; daa[1120+ 6] = 181.00; daa[1120+ 7] = 27.00; daa[1120+ 8] = 45.00;
418	daa[1120+ 9] = 21.00; daa[1120+10] = 14.00; daa[1220+ 0] = 54.00; daa[1220+ 1] = 44.00;
419	daa[1220+ 2] = 30.00; daa[1220+ 3] = 15.00; daa[1220+ 4] = 31.00; daa[1220+ 5] = 43.00;
420	daa[1220+ 6] = 18.00; daa[1220+ 7] = 14.00; daa[1220+ 8] = 33.00; daa[1220+ 9] = 479.00;
421	daa[1220+10] = 388.00; daa[1220+11] = 65.00; daa[1320+ 0] = 15.00; daa[1320+ 1] = 5.00;
422	daa[1320+ 2] = 10.00; daa[1320+ 3] = 4.00; daa[1320+ 4] = 78.00; daa[1320+ 5] = 4.00;
423	daa[1320+ 6] = 5.00; daa[1320+ 7] = 5.00; daa[1320+ 8] = 40.00; daa[1320+ 9] = 89.00;
424	daa[1320+10] = 248.00; daa[1320+11] = 4.00; daa[1320+12] = 43.00; daa[1420+ 0] = 194.00;
425	daa[1420+ 1] = 74.00; daa[1420+ 2] = 15.00; daa[1420+ 3] = 15.00; daa[1420+ 4] = 14.00;
426	daa[1420+ 5] = 164.00; daa[1420+ 6] = 18.00; daa[1420+ 7] = 24.00; daa[1420+ 8] = 115.00;
427	daa[1420+ 9] = 10.00; daa[1420+10] = 102.00; daa[1420+11] = 21.00; daa[1420+12] = 16.00;
428	daa[1420+13] = 17.00; daa[1520+ 0] = 378.00; daa[1520+ 1] = 101.00; daa[1520+ 2] = 503.00;
429	daa[1520+ 3] = 59.00; daa[1520+ 4] = 223.00; daa[1520+ 5] = 53.00; daa[1520+ 6] = 30.00;
430	daa[1520+ 7] = 201.00; daa[1520+ 8] = 73.00; daa[1520+ 9] = 40.00; daa[1520+10] = 59.00;
431	daa[1520+11] = 47.00; daa[1520+12] = 29.00; daa[1520+13] = 92.00; daa[1520+14] = 285.00;
432	daa[1620+ 0] = 475.00; daa[1620+ 1] = 64.00; daa[1620+ 2] = 232.00; daa[1620+ 3] = 38.00;
433	daa[1620+ 4] = 42.00; daa[1620+ 5] = 51.00; daa[1620+ 6] = 32.00; daa[1620+ 7] = 33.00;
434	daa[1620+ 8] = 46.00; daa[1620+ 9] = 245.00; daa[1620+10] = 25.00; daa[1620+11] = 103.00;
435	daa[1620+12] = 226.00; daa[1620+13] = 12.00; daa[1620+14] = 118.00; daa[1620+15] = 477.00;
436	daa[1720+ 0] = 9.00; daa[1720+ 1] = 126.00; daa[1720+ 2] = 8.00; daa[1720+ 3] = 4.00;
437	daa[1720+ 4] = 115.00; daa[1720+ 5] = 18.00; daa[1720+ 6] = 10.00; daa[1720+ 7] = 55.00;
438	daa[1720+ 8] = 8.00; daa[1720+ 9] = 9.00; daa[1720+10] = 52.00; daa[1720+11] = 10.00;
439	daa[1720+12] = 24.00; daa[1720+13] = 53.00; daa[1720+14] = 6.00; daa[1720+15] = 35.00;
440	daa[1720+16] = 12.00; daa[1820+ 0] = 11.00; daa[1820+ 1] = 20.00; daa[1820+ 2] = 70.00;
441	daa[1820+ 3] = 46.00; daa[1820+ 4] = 209.00; daa[1820+ 5] = 24.00; daa[1820+ 6] = 7.00;
442	daa[1820+ 7] = 8.00; daa[1820+ 8] = 573.00; daa[1820+ 9] = 32.00; daa[1820+10] = 24.00;
443	daa[1820+11] = 8.00; daa[1820+12] = 18.00; daa[1820+13] = 536.00; daa[1820+14] = 10.00;
444	daa[1820+15] = 63.00; daa[1820+16] = 21.00; daa[1820+17] = 71.00; daa[1920+ 0] = 298.00;
445	daa[1920+ 1] = 17.00; daa[1920+ 2] = 16.00; daa[1920+ 3] = 31.00; daa[1920+ 4] = 62.00;
446	daa[1920+ 5] = 20.00; daa[1920+ 6] = 45.00; daa[1920+ 7] = 47.00; daa[1920+ 8] = 11.00;
447	daa[1920+ 9] = 961.00; daa[1920+10] = 180.00; daa[1920+11] = 14.00; daa[1920+12] = 323.00;
448	daa[1920+13] = 62.00; daa[1920+14] = 23.00; daa[1920+15] = 38.00; daa[1920+16] = 112.00;
449	daa[1920+17] = 25.00; daa[1920+18] = 16.00;
450
451	/* f[ 0] = 0.076748; f[ 1] = 0.051691; f[ 2] = 0.042645; f[ 3] = 0.051544;
452	f[ 4] = 0.019803; f[ 5] = 0.040752; f[ 6] = 0.061830; f[ 7] = 0.073152;
453	f[ 8] = 0.022944; f[ 9] = 0.053761; f[10] = 0.091904; f[11] = 0.058676;
454	f[12] = 0.023826; f[13] = 0.040126; f[14] = 0.050901; f[15] = 0.068765;
455	f[16] = 0.058565; f[17] = 0.014261; f[18] = 0.032102; f[19] = 0.066005;
456	*/
457
458	f[ 0] = 0.07700; f[ 1] = 0.05200; f[ 2] = 0.04200; f[ 3] = 0.05100;
459	f[ 4] = 0.02000; f[ 5] = 0.04100; f[ 6] = 0.06200; f[ 7] = 0.07300;
460	f[ 8] = 0.02300; f[ 9] = 0.05400; f[10] = 0.09200; f[11] = 0.05900;
461	f[12] = 0.02400; f[13] = 0.04000; f[14] = 0.05100; f[15] = 0.06900;
462	f[16] = 0.05800; f[17] = 0.01400; f[18] = 0.03200; f[19] = 0.06600;
463
464	}
465	break;
466	case MTREV:
467	{
468	daa[ 120+ 0] = 23.18; daa[ 220+ 0] = 26.95; daa[ 220+ 1] = 13.24; daa[ 320+ 0] = 17.67;
469	daa[ 320+ 1] = 1.90; daa[ 320+ 2] = 794.38; daa[ 420+ 0] = 59.93; daa[ 420+ 1] = 103.33;
470	daa[ 420+ 2] = 58.94; daa[ 420+ 3] = 1.90; daa[ 520+ 0] = 1.90; daa[ 520+ 1] = 220.99;
471	daa[ 520+ 2] = 173.56; daa[ 520+ 3] = 55.28; daa[ 520+ 4] = 75.24; daa[ 620+ 0] = 9.77;
472	daa[ 620+ 1] = 1.90; daa[ 620+ 2] = 63.05; daa[ 620+ 3] = 583.55; daa[ 620+ 4] = 1.90;
473	daa[ 620+ 5] = 313.56; daa[ 720+ 0] = 120.71; daa[ 720+ 1] = 23.03; daa[ 720+ 2] = 53.30;
474	daa[ 720+ 3] = 56.77; daa[ 720+ 4] = 30.71; daa[ 720+ 5] = 6.75; daa[ 720+ 6] = 28.28;
475	daa[ 820+ 0] = 13.90; daa[ 820+ 1] = 165.23; daa[ 820+ 2] = 496.13; daa[ 820+ 3] = 113.99;
476	daa[ 820+ 4] = 141.49; daa[ 820+ 5] = 582.40; daa[ 820+ 6] = 49.12; daa[ 820+ 7] = 1.90;
477	daa[ 920+ 0] = 96.49; daa[ 920+ 1] = 1.90; daa[ 920+ 2] = 27.10; daa[ 920+ 3] = 4.34;
478	daa[ 920+ 4] = 62.73; daa[ 920+ 5] = 8.34; daa[ 920+ 6] = 3.31; daa[ 920+ 7] = 5.98;
479	daa[ 920+ 8] = 12.26; daa[1020+ 0] = 25.46; daa[1020+ 1] = 15.58; daa[1020+ 2] = 15.16;
480	daa[1020+ 3] = 1.90; daa[1020+ 4] = 25.65; daa[1020+ 5] = 39.70; daa[1020+ 6] = 1.90;
481	daa[1020+ 7] = 2.41; daa[1020+ 8] = 11.49; daa[1020+ 9] = 329.09; daa[1120+ 0] = 8.36;
482	daa[1120+ 1] = 141.40; daa[1120+ 2] = 608.70; daa[1120+ 3] = 2.31; daa[1120+ 4] = 1.90;
483	daa[1120+ 5] = 465.58; daa[1120+ 6] = 313.86; daa[1120+ 7] = 22.73; daa[1120+ 8] = 127.67;
484	daa[1120+ 9] = 19.57; daa[1120+10] = 14.88; daa[1220+ 0] = 141.88; daa[1220+ 1] = 1.90;
485	daa[1220+ 2] = 65.41; daa[1220+ 3] = 1.90; daa[1220+ 4] = 6.18; daa[1220+ 5] = 47.37;
486	daa[1220+ 6] = 1.90; daa[1220+ 7] = 1.90; daa[1220+ 8] = 11.97; daa[1220+ 9] = 517.98;
487	daa[1220+10] = 537.53; daa[1220+11] = 91.37; daa[1320+ 0] = 6.37; daa[1320+ 1] = 4.69;
488	daa[1320+ 2] = 15.20; daa[1320+ 3] = 4.98; daa[1320+ 4] = 70.80; daa[1320+ 5] = 19.11;
489	daa[1320+ 6] = 2.67; daa[1320+ 7] = 1.90; daa[1320+ 8] = 48.16; daa[1320+ 9] = 84.67;
490	daa[1320+10] = 216.06; daa[1320+11] = 6.44; daa[1320+12] = 90.82; daa[1420+ 0] = 54.31;
491	daa[1420+ 1] = 23.64; daa[1420+ 2] = 73.31; daa[1420+ 3] = 13.43; daa[1420+ 4] = 31.26;
492	daa[1420+ 5] = 137.29; daa[1420+ 6] = 12.83; daa[1420+ 7] = 1.90; daa[1420+ 8] = 60.97;
493	daa[1420+ 9] = 20.63; daa[1420+10] = 40.10; daa[1420+11] = 50.10; daa[1420+12] = 18.84;
494	daa[1420+13] = 17.31; daa[1520+ 0] = 387.86; daa[1520+ 1] = 6.04; daa[1520+ 2] = 494.39;
495	daa[1520+ 3] = 69.02; daa[1520+ 4] = 277.05; daa[1520+ 5] = 54.11; daa[1520+ 6] = 54.71;
496	daa[1520+ 7] = 125.93; daa[1520+ 8] = 77.46; daa[1520+ 9] = 47.70; daa[1520+10] = 73.61;
497	daa[1520+11] = 105.79; daa[1520+12] = 111.16; daa[1520+13] = 64.29; daa[1520+14] = 169.90;
498	daa[1620+ 0] = 480.72; daa[1620+ 1] = 2.08; daa[1620+ 2] = 238.46; daa[1620+ 3] = 28.01;
499	daa[1620+ 4] = 179.97; daa[1620+ 5] = 94.93; daa[1620+ 6] = 14.82; daa[1620+ 7] = 11.17;
500	daa[1620+ 8] = 44.78; daa[1620+ 9] = 368.43; daa[1620+10] = 126.40; daa[1620+11] = 136.33;
501	daa[1620+12] = 528.17; daa[1620+13] = 33.85; daa[1620+14] = 128.22; daa[1620+15] = 597.21;
502	daa[1720+ 0] = 1.90; daa[1720+ 1] = 21.95; daa[1720+ 2] = 10.68; daa[1720+ 3] = 19.86;
503	daa[1720+ 4] = 33.60; daa[1720+ 5] = 1.90; daa[1720+ 6] = 1.90; daa[1720+ 7] = 10.92;
504	daa[1720+ 8] = 7.08; daa[1720+ 9] = 1.90; daa[1720+10] = 32.44; daa[1720+11] = 24.00;
505	daa[1720+12] = 21.71; daa[1720+13] = 7.84; daa[1720+14] = 4.21; daa[1720+15] = 38.58;
506	daa[1720+16] = 9.99; daa[1820+ 0] = 6.48; daa[1820+ 1] = 1.90; daa[1820+ 2] = 191.36;
507	daa[1820+ 3] = 21.21; daa[1820+ 4] = 254.77; daa[1820+ 5] = 38.82; daa[1820+ 6] = 13.12;
508	daa[1820+ 7] = 3.21; daa[1820+ 8] = 670.14; daa[1820+ 9] = 25.01; daa[1820+10] = 44.15;
509	daa[1820+11] = 51.17; daa[1820+12] = 39.96; daa[1820+13] = 465.58; daa[1820+14] = 16.21;
510	daa[1820+15] = 64.92; daa[1820+16] = 38.73; daa[1820+17] = 26.25; daa[1920+ 0] = 195.06;
511	daa[1920+ 1] = 7.64; daa[1920+ 2] = 1.90; daa[1920+ 3] = 1.90; daa[1920+ 4] = 1.90;
512	daa[1920+ 5] = 19.00; daa[1920+ 6] = 21.14; daa[1920+ 7] = 2.53; daa[1920+ 8] = 1.90;
513	daa[1920+ 9] = 1222.94; daa[1920+10] = 91.67; daa[1920+11] = 1.90; daa[1920+12] = 387.54;
514	daa[1920+13] = 6.35; daa[1920+14] = 8.23; daa[1920+15] = 1.90; daa[1920+16] = 204.54;
515	daa[1920+17] = 5.37; daa[1920+18] = 1.90;
516
517
518	f[ 0] = 0.072000; f[ 1] = 0.019000; f[ 2] = 0.039000; f[ 3] = 0.019000;
519	f[ 4] = 0.006000; f[ 5] = 0.025000; f[ 6] = 0.024000; f[ 7] = 0.056000;
520	f[ 8] = 0.028000; f[ 9] = 0.088000; f[10] = 0.169000; f[11] = 0.023000;
521	f[12] = 0.054000; f[13] = 0.061000; f[14] = 0.054000; f[15] = 0.072000;
522	f[16] = 0.086000; f[17] = 0.029000; f[18] = 0.033000; f[19] = 0.043000;
523	}
524	break;
525	case WAG:
526	{
527	daa[ 120+ 0] = 55.15710; daa[ 220+ 0] = 50.98480; daa[ 2*20+ 1] = 63.53460;
528	daa[ 320+ 0] = 73.89980; daa[ 320+ 1] = 14.73040; daa[ 3*20+ 2] = 542.94200;
529	daa[ 420+ 0] = 102.70400; daa[ 420+ 1] = 52.81910; daa[ 4*20+ 2] = 26.52560;
530	daa[ 420+ 3] = 3.02949; daa[ 520+ 0] = 90.85980; daa[ 5*20+ 1] = 303.55000;
531	daa[ 520+ 2] = 154.36400; daa[ 520+ 3] = 61.67830; daa[ 5*20+ 4] = 9.88179;
532	daa[ 620+ 0] = 158.28500; daa[ 620+ 1] = 43.91570; daa[ 6*20+ 2] = 94.71980;
533	daa[ 620+ 3] = 617.41600; daa[ 620+ 4] = 2.13520; daa[ 6*20+ 5] = 546.94700;
534	daa[ 720+ 0] = 141.67200; daa[ 720+ 1] = 58.46650; daa[ 7*20+ 2] = 112.55600;
535	daa[ 720+ 3] = 86.55840; daa[ 720+ 4] = 30.66740; daa[ 7*20+ 5] = 33.00520;
536	daa[ 720+ 6] = 56.77170; daa[ 820+ 0] = 31.69540; daa[ 8*20+ 1] = 213.71500;
537	daa[ 820+ 2] = 395.62900; daa[ 820+ 3] = 93.06760; daa[ 8*20+ 4] = 24.89720;
538	daa[ 820+ 5] = 429.41100; daa[ 820+ 6] = 57.00250; daa[ 8*20+ 7] = 24.94100;
539	daa[ 920+ 0] = 19.33350; daa[ 920+ 1] = 18.69790; daa[ 9*20+ 2] = 55.42360;
540	daa[ 920+ 3] = 3.94370; daa[ 920+ 4] = 17.01350; daa[ 9*20+ 5] = 11.39170;
541	daa[ 920+ 6] = 12.73950; daa[ 920+ 7] = 3.04501; daa[ 9*20+ 8] = 13.81900;
542	daa[1020+ 0] = 39.79150; daa[1020+ 1] = 49.76710; daa[10*20+ 2] = 13.15280;
543	daa[1020+ 3] = 8.48047; daa[1020+ 4] = 38.42870; daa[10*20+ 5] = 86.94890;
544	daa[1020+ 6] = 15.42630; daa[1020+ 7] = 6.13037; daa[10*20+ 8] = 49.94620;
545	daa[1020+ 9] = 317.09700; daa[1120+ 0] = 90.62650; daa[11*20+ 1] = 535.14200;
546	daa[1120+ 2] = 301.20100; daa[1120+ 3] = 47.98550; daa[11*20+ 4] = 7.40339;
547	daa[1120+ 5] = 389.49000; daa[1120+ 6] = 258.44300; daa[11*20+ 7] = 37.35580;
548	daa[1120+ 8] = 89.04320; daa[1120+ 9] = 32.38320; daa[11*20+10] = 25.75550;
549	daa[1220+ 0] = 89.34960; daa[1220+ 1] = 68.31620; daa[12*20+ 2] = 19.82210;
550	daa[1220+ 3] = 10.37540; daa[1220+ 4] = 39.04820; daa[12*20+ 5] = 154.52600;
551	daa[1220+ 6] = 31.51240; daa[1220+ 7] = 17.41000; daa[12*20+ 8] = 40.41410;
552	daa[1220+ 9] = 425.74600; daa[1220+10] = 485.40200; daa[12*20+11] = 93.42760;
553	daa[1320+ 0] = 21.04940; daa[1320+ 1] = 10.27110; daa[13*20+ 2] = 9.61621;
554	daa[1320+ 3] = 4.67304; daa[1320+ 4] = 39.80200; daa[13*20+ 5] = 9.99208;
555	daa[1320+ 6] = 8.11339; daa[1320+ 7] = 4.99310; daa[13*20+ 8] = 67.93710;
556	daa[1320+ 9] = 105.94700; daa[1320+10] = 211.51700; daa[13*20+11] = 8.88360;
557	daa[1320+12] = 119.06300; daa[1420+ 0] = 143.85500; daa[14*20+ 1] = 67.94890;
558	daa[1420+ 2] = 19.50810; daa[1420+ 3] = 42.39840; daa[14*20+ 4] = 10.94040;
559	daa[1420+ 5] = 93.33720; daa[1420+ 6] = 68.23550; daa[14*20+ 7] = 24.35700;
560	daa[1420+ 8] = 69.61980; daa[1420+ 9] = 9.99288; daa[14*20+10] = 41.58440;
561	daa[1420+11] = 55.68960; daa[1420+12] = 17.13290; daa[14*20+13] = 16.14440;
562	daa[1520+ 0] = 337.07900; daa[1520+ 1] = 122.41900; daa[15*20+ 2] = 397.42300;
563	daa[1520+ 3] = 107.17600; daa[1520+ 4] = 140.76600; daa[15*20+ 5] = 102.88700;
564	daa[1520+ 6] = 70.49390; daa[1520+ 7] = 134.18200; daa[15*20+ 8] = 74.01690;
565	daa[1520+ 9] = 31.94400; daa[1520+10] = 34.47390; daa[15*20+11] = 96.71300;
566	daa[1520+12] = 49.39050; daa[1520+13] = 54.59310; daa[15*20+14] = 161.32800;
567	daa[1620+ 0] = 212.11100; daa[1620+ 1] = 55.44130; daa[16*20+ 2] = 203.00600;
568	daa[1620+ 3] = 37.48660; daa[1620+ 4] = 51.29840; daa[16*20+ 5] = 85.79280;
569	daa[1620+ 6] = 82.27650; daa[1620+ 7] = 22.58330; daa[16*20+ 8] = 47.33070;
570	daa[1620+ 9] = 145.81600; daa[1620+10] = 32.66220; daa[16*20+11] = 138.69800;
571	daa[1620+12] = 151.61200; daa[1620+13] = 17.19030; daa[16*20+14] = 79.53840;
572	daa[1620+15] = 437.80200; daa[1720+ 0] = 11.31330; daa[17*20+ 1] = 116.39200;
573	daa[1720+ 2] = 7.19167; daa[1720+ 3] = 12.97670; daa[17*20+ 4] = 71.70700;
574	daa[1720+ 5] = 21.57370; daa[1720+ 6] = 15.65570; daa[17*20+ 7] = 33.69830;
575	daa[1720+ 8] = 26.25690; daa[1720+ 9] = 21.24830; daa[17*20+10] = 66.53090;
576	daa[1720+11] = 13.75050; daa[1720+12] = 51.57060; daa[17*20+13] = 152.96400;
577	daa[1720+14] = 13.94050; daa[1720+15] = 52.37420; daa[17*20+16] = 11.08640;
578	daa[1820+ 0] = 24.07350; daa[1820+ 1] = 38.15330; daa[18*20+ 2] = 108.60000;
579	daa[1820+ 3] = 32.57110; daa[1820+ 4] = 54.38330; daa[18*20+ 5] = 22.77100;
580	daa[1820+ 6] = 19.63030; daa[1820+ 7] = 10.36040; daa[18*20+ 8] = 387.34400;
581	daa[1820+ 9] = 42.01700; daa[1820+10] = 39.86180; daa[18*20+11] = 13.32640;
582	daa[1820+12] = 42.84370; daa[1820+13] = 645.42800; daa[18*20+14] = 21.60460;
583	daa[1820+15] = 78.69930; daa[1820+16] = 29.11480; daa[18*20+17] = 248.53900;
584	daa[1920+ 0] = 200.60100; daa[1920+ 1] = 25.18490; daa[19*20+ 2] = 19.62460;
585	daa[1920+ 3] = 15.23350; daa[1920+ 4] = 100.21400; daa[19*20+ 5] = 30.12810;
586	daa[1920+ 6] = 58.87310; daa[1920+ 7] = 18.72470; daa[19*20+ 8] = 11.83580;
587	daa[1920+ 9] = 782.13000; daa[1920+10] = 180.03400; daa[19*20+11] = 30.54340;
588	daa[1920+12] = 205.84500; daa[1920+13] = 64.98920; daa[19*20+14] = 31.48870;
589	daa[1920+15] = 23.27390; daa[1920+16] = 138.82300; daa[19*20+17] = 36.53690;
590	daa[19*20+18] = 31.47300;
591
592	/* f[0] = 0.0866279; f[1] = 0.043972; f[2] = 0.0390894; f[3] = 0.0570451;
593	f[4] = 0.0193078; f[5] = 0.0367281; f[6] = 0.0580589; f[7] = 0.0832518;
594	f[8] = 0.0244313; f[9] = 0.048466; f[10] = 0.086209; f[11] = 0.0620286;
595	f[12] = 0.0195027; f[13] = 0.0384319; f[14] = 0.0457631; f[15] = 0.0695179;
596	f[16] = 0.0610127; f[17] = 0.0143859; f[18] = 0.0352742; f[19] = 0.0708956;
597	*/
598	f[0] = 0.08700; f[1] = 0.04400; f[2] = 0.03900; f[3] = 0.05700;
599	f[4] = 0.01900; f[5] = 0.03700; f[6] = 0.05800; f[7] = 0.08300;
600	f[8] = 0.02400; f[9] = 0.04900; f[10] = 0.08600; f[11] = 0.06200;
601	f[12] = 0.02000; f[13] = 0.03800; f[14] = 0.04600; f[15] = 0.07000;
602	f[16] = 0.06100; f[17] = 0.01400; f[18] = 0.03500; f[19] = 0.07100;
603	}
604	break;
605	case RTREV:
606	{
607	daa[120+0]= 34; daa[220+0]= 51; daa[220+1]= 35; daa[320+0]= 10;
608	daa[320+1]= 30; daa[320+2]= 384; daa[420+0]= 439; daa[420+1]= 92;
609	daa[420+2]= 128; daa[420+3]= 1; daa[520+0]= 32; daa[520+1]= 221;
610	daa[520+2]= 236; daa[520+3]= 78; daa[520+4]= 70; daa[620+0]= 81;
611	daa[620+1]= 10; daa[620+2]= 79; daa[620+3]= 542; daa[620+4]= 1;
612	daa[620+5]= 372; daa[720+0]= 135; daa[720+1]= 41; daa[720+2]= 94;
613	daa[720+3]= 61; daa[720+4]= 48; daa[720+5]= 18; daa[720+6]= 70;
614	daa[820+0]= 30; daa[820+1]= 90; daa[820+2]= 320; daa[820+3]= 91;
615	daa[820+4]= 124; daa[820+5]= 387; daa[820+6]= 34; daa[820+7]= 68;
616	daa[920+0]= 1; daa[920+1]= 24; daa[920+2]= 35; daa[920+3]= 1;
617	daa[920+4]= 104; daa[920+5]= 33; daa[920+6]= 1; daa[920+7]= 1;
618	daa[920+8]= 34; daa[1020+0]= 45; daa[1020+1]= 18; daa[1020+2]= 15;
619	daa[1020+3]= 5; daa[1020+4]= 110; daa[1020+5]= 54; daa[1020+6]= 21;
620	daa[1020+7]= 3; daa[1020+8]= 51; daa[1020+9]= 385; daa[1120+0]= 38;
621	daa[1120+1]= 593; daa[1120+2]= 123; daa[1120+3]= 20; daa[1120+4]= 16;
622	daa[1120+5]= 309; daa[1120+6]= 141; daa[1120+7]= 30; daa[1120+8]= 76;
623	daa[1120+9]= 34; daa[1120+10]= 23; daa[1220+0]= 235; daa[1220+1]= 57;
624	daa[1220+2]= 1; daa[1220+3]= 1; daa[1220+4]= 156; daa[1220+5]= 158;
625	daa[1220+6]= 1; daa[1220+7]= 37; daa[1220+8]= 116; daa[1220+9]= 375;
626	daa[1220+10]= 581; daa[1220+11]= 134; daa[1320+0]= 1; daa[1320+1]= 7;
627	daa[1320+2]= 49; daa[1320+3]= 1; daa[1320+4]= 70; daa[1320+5]= 1;
628	daa[1320+6]= 1; daa[1320+7]= 7; daa[1320+8]= 141; daa[1320+9]= 64;
629	daa[1320+10]= 179; daa[1320+11]= 14; daa[1320+12]= 247; daa[1420+0]= 97;
630	daa[1420+1]= 24; daa[1420+2]= 33; daa[1420+3]= 55; daa[1420+4]= 1;
631	daa[1420+5]= 68; daa[1420+6]= 52; daa[1420+7]= 17; daa[1420+8]= 44;
632	daa[1420+9]= 10; daa[1420+10]= 22; daa[1420+11]= 43; daa[1420+12]= 1;
633	daa[1420+13]= 11; daa[1520+0]= 460; daa[1520+1]= 102; daa[1520+2]= 294;
634	daa[1520+3]= 136; daa[1520+4]= 75; daa[1520+5]= 225; daa[1520+6]= 95;
635	daa[1520+7]= 152; daa[1520+8]= 183; daa[1520+9]= 4; daa[1520+10]= 24;
636	daa[1520+11]= 77; daa[1520+12]= 1; daa[1520+13]= 20; daa[1520+14]= 134;
637	daa[1620+0]= 258; daa[1620+1]= 64; daa[1620+2]= 148; daa[1620+3]= 55;
638	daa[1620+4]= 117; daa[1620+5]= 146; daa[1620+6]= 82; daa[1620+7]= 7;
639	daa[1620+8]= 49; daa[1620+9]= 72; daa[1620+10]= 25; daa[1620+11]= 110;
640	daa[1620+12]= 131; daa[1620+13]= 69; daa[1620+14]= 62; daa[1620+15]= 671;
641	daa[1720+0]= 5; daa[1720+1]= 13; daa[1720+2]= 16; daa[1720+3]= 1;
642	daa[1720+4]= 55; daa[1720+5]= 10; daa[1720+6]= 17; daa[1720+7]= 23;
643	daa[1720+8]= 48; daa[1720+9]= 39; daa[1720+10]= 47; daa[1720+11]= 6;
644	daa[1720+12]= 111; daa[1720+13]= 182; daa[1720+14]= 9; daa[1720+15]= 14;
645	daa[1720+16]= 1; daa[1820+0]= 55; daa[1820+1]= 47; daa[1820+2]= 28;
646	daa[1820+3]= 1; daa[1820+4]= 131; daa[1820+5]= 45; daa[1820+6]= 1;
647	daa[1820+7]= 21; daa[1820+8]= 307; daa[1820+9]= 26; daa[1820+10]= 64;
648	daa[1820+11]= 1; daa[1820+12]= 74; daa[1820+13]= 1017; daa[1820+14]= 14;
649	daa[1820+15]= 31; daa[1820+16]= 34; daa[1820+17]= 176; daa[1920+0]= 197;
650	daa[1920+1]= 29; daa[1920+2]= 21; daa[1920+3]= 6; daa[1920+4]= 295;
651	daa[1920+5]= 36; daa[1920+6]= 35; daa[1920+7]= 3; daa[1920+8]= 1;
652	daa[1920+9]= 1048; daa[1920+10]= 112; daa[1920+11]= 19; daa[1920+12]= 236;
653	daa[1920+13]= 92; daa[1920+14]= 25; daa[1920+15]= 39; daa[1920+16]= 196;
654	daa[1920+17]= 26; daa[1920+18]= 59;
655
656	f[0]= 0.0646; f[1]= 0.0453; f[2]= 0.0376; f[3]= 0.0422;
657	f[4]= 0.0114; f[5]= 0.0606; f[6]= 0.0607; f[7]= 0.0639;
658	f[8]= 0.0273; f[9]= 0.0679; f[10]= 0.1018; f[11]= 0.0751;
659	f[12]= 0.015; f[13]= 0.0287; f[14]= 0.0681; f[15]= 0.0488;
660	f[16]= 0.0622; f[17]= 0.0251; f[18]= 0.0318; f[19]= 0.0619;
661	}
662	break;
663	case CPREV:
664	{
665	daa[120+0]= 105; daa[220+0]= 227; daa[220+1]= 357; daa[320+0]= 175;
666	daa[320+1]= 43; daa[320+2]= 4435; daa[420+0]= 669; daa[420+1]= 823;
667	daa[420+2]= 538; daa[420+3]= 10; daa[520+0]= 157; daa[520+1]= 1745;
668	daa[520+2]= 768; daa[520+3]= 400; daa[520+4]= 10; daa[620+0]= 499;
669	daa[620+1]= 152; daa[620+2]= 1055; daa[620+3]= 3691; daa[620+4]= 10;
670	daa[620+5]= 3122; daa[720+0]= 665; daa[720+1]= 243; daa[720+2]= 653;
671	daa[720+3]= 431; daa[720+4]= 303; daa[720+5]= 133; daa[720+6]= 379;
672	daa[820+0]= 66; daa[820+1]= 715; daa[820+2]= 1405; daa[820+3]= 331;
673	daa[820+4]= 441; daa[820+5]= 1269; daa[820+6]= 162; daa[820+7]= 19;
674	daa[920+0]= 145; daa[920+1]= 136; daa[920+2]= 168; daa[920+3]= 10;
675	daa[920+4]= 280; daa[920+5]= 92; daa[920+6]= 148; daa[920+7]= 40;
676	daa[920+8]= 29; daa[1020+0]= 197; daa[1020+1]= 203; daa[1020+2]= 113;
677	daa[1020+3]= 10; daa[1020+4]= 396; daa[1020+5]= 286; daa[1020+6]= 82;
678	daa[1020+7]= 20; daa[1020+8]= 66; daa[1020+9]= 1745; daa[1120+0]= 236;
679	daa[1120+1]= 4482; daa[1120+2]= 2430; daa[1120+3]= 412; daa[1120+4]= 48;
680	daa[1120+5]= 3313; daa[1120+6]= 2629; daa[1120+7]= 263; daa[1120+8]= 305;
681	daa[1120+9]= 345; daa[1120+10]= 218; daa[1220+0]= 185; daa[1220+1]= 125;
682	daa[1220+2]= 61; daa[1220+3]= 47; daa[1220+4]= 159; daa[1220+5]= 202;
683	daa[1220+6]= 113; daa[1220+7]= 21; daa[1220+8]= 10; daa[1220+9]= 1772;
684	daa[1220+10]= 1351; daa[1220+11]= 193; daa[1320+0]= 68; daa[1320+1]= 53;
685	daa[1320+2]= 97; daa[1320+3]= 22; daa[1320+4]= 726; daa[1320+5]= 10;
686	daa[1320+6]= 145; daa[1320+7]= 25; daa[1320+8]= 127; daa[1320+9]= 454;
687	daa[1320+10]= 1268; daa[1320+11]= 72; daa[1320+12]= 327; daa[1420+0]= 490;
688	daa[1420+1]= 87; daa[1420+2]= 173; daa[1420+3]= 170; daa[1420+4]= 285;
689	daa[1420+5]= 323; daa[1420+6]= 185; daa[1420+7]= 28; daa[1420+8]= 152;
690	daa[1420+9]= 117; daa[1420+10]= 219; daa[1420+11]= 302; daa[1420+12]= 100;
691	daa[1420+13]= 43; daa[1520+0]= 2440; daa[1520+1]= 385; daa[1520+2]= 2085;
692	daa[1520+3]= 590; daa[1520+4]= 2331; daa[1520+5]= 396; daa[1520+6]= 568;
693	daa[1520+7]= 691; daa[1520+8]= 303; daa[1520+9]= 216; daa[1520+10]= 516;
694	daa[1520+11]= 868; daa[1520+12]= 93; daa[1520+13]= 487; daa[1520+14]= 1202;
695	daa[1620+0]= 1340; daa[1620+1]= 314; daa[1620+2]= 1393; daa[1620+3]= 266;
696	daa[1620+4]= 576; daa[1620+5]= 241; daa[1620+6]= 369; daa[1620+7]= 92;
697	daa[1620+8]= 32; daa[1620+9]= 1040; daa[1620+10]= 156; daa[1620+11]= 918;
698	daa[1620+12]= 645; daa[1620+13]= 148; daa[1620+14]= 260; daa[1620+15]= 2151;
699	daa[1720+0]= 14; daa[1720+1]= 230; daa[1720+2]= 40; daa[1720+3]= 18;
700	daa[1720+4]= 435; daa[1720+5]= 53; daa[1720+6]= 63; daa[1720+7]= 82;
701	daa[1720+8]= 69; daa[1720+9]= 42; daa[1720+10]= 159; daa[1720+11]= 10;
702	daa[1720+12]= 86; daa[1720+13]= 468; daa[1720+14]= 49; daa[1720+15]= 73;
703	daa[1720+16]= 29; daa[1820+0]= 56; daa[1820+1]= 323; daa[1820+2]= 754;
704	daa[1820+3]= 281; daa[1820+4]= 1466; daa[1820+5]= 391; daa[1820+6]= 142;
705	daa[1820+7]= 10; daa[1820+8]= 1971; daa[1820+9]= 89; daa[1820+10]= 189;
706	daa[1820+11]= 247; daa[1820+12]= 215; daa[1820+13]= 2370; daa[1820+14]= 97;
707	daa[1820+15]= 522; daa[1820+16]= 71; daa[1820+17]= 346; daa[1920+0]= 968;
708	daa[1920+1]= 92; daa[1920+2]= 83; daa[1920+3]= 75; daa[1920+4]= 592;
709	daa[1920+5]= 54; daa[1920+6]= 200; daa[1920+7]= 91; daa[1920+8]= 25;
710	daa[1920+9]= 4797; daa[1920+10]= 865; daa[1920+11]= 249; daa[1920+12]= 475;
711	daa[1920+13]= 317; daa[1920+14]= 122; daa[1920+15]= 167; daa[1920+16]= 760;
712	daa[1920+17]= 10; daa[1920+18]= 119;
713
714	f[0]= 0.076; f[1]= 0.062; f[2]= 0.041; f[3]= 0.037;
715	f[4]= 0.009; f[5]= 0.038; f[6]= 0.049; f[7]= 0.084;
716	f[8]= 0.025; f[9]= 0.081; f[10]= 0.101; f[11]= 0.05;
717	f[12]= 0.022; f[13]= 0.051; f[14]= 0.043; f[15]= 0.062;
718	f[16]= 0.054; f[17]= 0.018; f[18]= 0.031; f[19]= 0.066;
719	}
720	break;
721	case VT:
722	{
723	daa[120+0]= 0.233108; daa[220+0]= 0.199097; daa[220+1]= 0.210797; daa[320+0]= 0.265145;
724	daa[320+1]= 0.105191; daa[320+2]= 0.883422; daa[420+0]= 0.227333; daa[420+1]= 0.031726;
725	daa[420+2]= 0.027495; daa[420+3]= 0.010313; daa[520+0]= 0.310084; daa[520+1]= 0.493763;
726	daa[520+2]= 0.2757; daa[520+3]= 0.205842; daa[520+4]= 0.004315; daa[620+0]= 0.567957;
727	daa[620+1]= 0.25524; daa[620+2]= 0.270417; daa[620+3]= 1.599461; daa[620+4]= 0.005321;
728	daa[620+5]= 0.960976; daa[720+0]= 0.876213; daa[720+1]= 0.156945; daa[720+2]= 0.362028;
729	daa[720+3]= 0.311718; daa[720+4]= 0.050876; daa[720+5]= 0.12866; daa[720+6]= 0.250447;
730	daa[820+0]= 0.078692; daa[820+1]= 0.213164; daa[820+2]= 0.290006; daa[820+3]= 0.134252;
731	daa[820+4]= 0.016695; daa[820+5]= 0.315521; daa[820+6]= 0.104458; daa[820+7]= 0.058131;
732	daa[920+0]= 0.222972; daa[920+1]= 0.08151; daa[920+2]= 0.087225; daa[920+3]= 0.01172;
733	daa[920+4]= 0.046398; daa[920+5]= 0.054602; daa[920+6]= 0.046589; daa[920+7]= 0.051089;
734	daa[920+8]= 0.020039; daa[1020+0]= 0.42463; daa[1020+1]= 0.192364; daa[1020+2]= 0.069245;
735	daa[1020+3]= 0.060863; daa[1020+4]= 0.091709; daa[1020+5]= 0.24353; daa[1020+6]= 0.151924;
736	daa[1020+7]= 0.087056; daa[1020+8]= 0.103552; daa[1020+9]= 2.08989; daa[1120+0]= 0.393245;
737	daa[1120+1]= 1.755838; daa[1120+2]= 0.50306; daa[1120+3]= 0.261101; daa[1120+4]= 0.004067;
738	daa[1120+5]= 0.738208; daa[1120+6]= 0.88863; daa[1120+7]= 0.193243; daa[1120+8]= 0.153323;
739	daa[1120+9]= 0.093181; daa[1120+10]= 0.201204; daa[1220+0]= 0.21155; daa[1220+1]= 0.08793;
740	daa[1220+2]= 0.05742; daa[1220+3]= 0.012182; daa[1220+4]= 0.02369; daa[1220+5]= 0.120801;
741	daa[1220+6]= 0.058643; daa[1220+7]= 0.04656; daa[1220+8]= 0.021157; daa[1220+9]= 0.493845;
742	daa[1220+10]= 1.105667; daa[1220+11]= 0.096474; daa[1320+0]= 0.116646; daa[1320+1]= 0.042569;
743	daa[1320+2]= 0.039769; daa[1320+3]= 0.016577; daa[1320+4]= 0.051127; daa[1320+5]= 0.026235;
744	daa[1320+6]= 0.028168; daa[1320+7]= 0.050143; daa[1320+8]= 0.079807; daa[1320+9]= 0.32102;
745	daa[1320+10]= 0.946499; daa[1320+11]= 0.038261; daa[1320+12]= 0.173052; daa[1420+0]= 0.399143;
746	daa[1420+1]= 0.12848; daa[1420+2]= 0.083956; daa[1420+3]= 0.160063; daa[1420+4]= 0.011137;
747	daa[1420+5]= 0.15657; daa[1420+6]= 0.205134; daa[1420+7]= 0.124492; daa[1420+8]= 0.078892;
748	daa[1420+9]= 0.054797; daa[1420+10]= 0.169784; daa[1420+11]= 0.212302; daa[1420+12]= 0.010363;
749	daa[1420+13]= 0.042564; daa[1520+0]= 1.817198; daa[1520+1]= 0.292327; daa[1520+2]= 0.847049;
750	daa[1520+3]= 0.461519; daa[1520+4]= 0.17527; daa[1520+5]= 0.358017; daa[1520+6]= 0.406035;
751	daa[1520+7]= 0.612843; daa[1520+8]= 0.167406; daa[1520+9]= 0.081567; daa[1520+10]= 0.214977;
752	daa[1520+11]= 0.400072; daa[1520+12]= 0.090515; daa[1520+13]= 0.138119; daa[1520+14]= 0.430431;
753	daa[1620+0]= 0.877877; daa[1620+1]= 0.204109; daa[1620+2]= 0.471268; daa[1620+3]= 0.178197;
754	daa[1620+4]= 0.079511; daa[1620+5]= 0.248992; daa[1620+6]= 0.321028; daa[1620+7]= 0.136266;
755	daa[1620+8]= 0.101117; daa[1620+9]= 0.376588; daa[1620+10]= 0.243227; daa[1620+11]= 0.446646;
756	daa[1620+12]= 0.184609; daa[1620+13]= 0.08587; daa[1620+14]= 0.207143; daa[1620+15]= 1.767766;
757	daa[1720+0]= 0.030309; daa[1720+1]= 0.046417; daa[1720+2]= 0.010459; daa[1720+3]= 0.011393;
758	daa[1720+4]= 0.007732; daa[1720+5]= 0.021248; daa[1720+6]= 0.018844; daa[1720+7]= 0.02399;
759	daa[1720+8]= 0.020009; daa[1720+9]= 0.034954; daa[1720+10]= 0.083439; daa[1720+11]= 0.023321;
760	daa[1720+12]= 0.022019; daa[1720+13]= 0.12805; daa[1720+14]= 0.014584; daa[1720+15]= 0.035933;
761	daa[1720+16]= 0.020437; daa[1820+0]= 0.087061; daa[1820+1]= 0.09701; daa[1820+2]= 0.093268;
762	daa[1820+3]= 0.051664; daa[1820+4]= 0.042823; daa[1820+5]= 0.062544; daa[1820+6]= 0.0552;
763	daa[1820+7]= 0.037568; daa[1820+8]= 0.286027; daa[1820+9]= 0.086237; daa[1820+10]= 0.189842;
764	daa[1820+11]= 0.068689; daa[1820+12]= 0.073223; daa[1820+13]= 0.898663; daa[1820+14]= 0.032043;
765	daa[1820+15]= 0.121979; daa[1820+16]= 0.094617; daa[1820+17]= 0.124746; daa[1920+0]= 1.230985;
766	daa[1920+1]= 0.113146; daa[1920+2]= 0.049824; daa[1920+3]= 0.048769; daa[1920+4]= 0.163831;
767	daa[1920+5]= 0.112027; daa[1920+6]= 0.205868; daa[1920+7]= 0.082579; daa[1920+8]= 0.068575;
768	daa[1920+9]= 3.65443; daa[1920+10]= 1.337571; daa[1920+11]= 0.144587; daa[1920+12]= 0.307309;
769	daa[1920+13]= 0.247329; daa[1920+14]= 0.129315; daa[1920+15]= 0.1277; daa[1920+16]= 0.740372;
770	daa[1920+17]= 0.022134; daa[1920+18]= 0.125733;
771
772	/*f[0]= 0.078837; f[1]= 0.051238; f[2]= 0.042313; f[3]= 0.053066;
773	f[4]= 0.015175; f[5]= 0.036713; f[6]= 0.061924; f[7]= 0.070852;
774	f[8]= 0.023082; f[9]= 0.062056; f[10]= 0.096371; f[11]= 0.057324;
775	f[12]= 0.023771; f[13]= 0.043296; f[14]= 0.043911; f[15]= 0.063403;
776	f[16]= 0.055897; f[17]= 0.013272; f[18]= 0.034399; f[19]= 0.073101; */
777
778	f[0] = 0.07900; f[1]= 0.05100; f[2] = 0.04200; f[3]= 0.05300;
779	f[4] = 0.01500; f[5]= 0.03700; f[6] = 0.06200; f[7]= 0.07100;
780	f[8] = 0.02300; f[9]= 0.06200; f[10] = 0.09600; f[11]= 0.05700;
781	f[12] = 0.02400; f[13]= 0.04300; f[14] = 0.04400; f[15]= 0.06400;
782	f[16] = 0.05600; f[17]= 0.01300; f[18] = 0.03500; f[19]= 0.07300;
783
784	}
785	break;
786	case BLOSUM62:
787	{
788	daa[120+0]= 0.735790389698; daa[220+0]= 0.485391055466; daa[2*20+1]= 1.297446705134;
789	daa[3*20+0]= 0.543161820899;
790	daa[320+1]= 0.500964408555; daa[320+2]= 3.180100048216; daa[4*20+0]= 1.45999531047;
791	daa[4*20+1]= 0.227826574209;
792	daa[420+2]= 0.397358949897; daa[420+3]= 0.240836614802; daa[5*20+0]= 1.199705704602;
793	daa[5*20+1]= 3.020833610064;
794	daa[520+2]= 1.839216146992; daa[520+3]= 1.190945703396; daa[5*20+4]= 0.32980150463;
795	daa[6*20+0]= 1.1709490428;
796	daa[620+1]= 1.36057419042; daa[620+2]= 1.24048850864; daa[6*20+3]= 3.761625208368;
797	daa[6*20+4]= 0.140748891814;
798	daa[620+5]= 5.528919177928; daa[720+0]= 1.95588357496; daa[7*20+1]= 0.418763308518;
799	daa[7*20+2]= 1.355872344485;
800	daa[720+3]= 0.798473248968; daa[720+4]= 0.418203192284; daa[7*20+5]= 0.609846305383;
801	daa[7*20+6]= 0.423579992176;
802	daa[820+0]= 0.716241444998; daa[820+1]= 1.456141166336; daa[8*20+2]= 2.414501434208;
803	daa[8*20+3]= 0.778142664022;
804	daa[820+4]= 0.354058109831; daa[820+5]= 2.43534113114; daa[8*20+6]= 1.626891056982;
805	daa[8*20+7]= 0.539859124954;
806	daa[920+0]= 0.605899003687; daa[920+1]= 0.232036445142; daa[9*20+2]= 0.283017326278;
807	daa[9*20+3]= 0.418555732462;
808	daa[920+4]= 0.774894022794; daa[920+5]= 0.236202451204; daa[9*20+6]= 0.186848046932;
809	daa[9*20+7]= 0.189296292376;
810	daa[920+8]= 0.252718447885; daa[1020+0]= 0.800016530518; daa[10*20+1]= 0.622711669692;
811	daa[10*20+2]= 0.211888159615;
812	daa[1020+3]= 0.218131577594; daa[1020+4]= 0.831842640142; daa[10*20+5]= 0.580737093181;
813	daa[10*20+6]= 0.372625175087;
814	daa[1020+7]= 0.217721159236; daa[1020+8]= 0.348072209797; daa[10*20+9]= 3.890963773304;
815	daa[11*20+0]= 1.295201266783;
816	daa[1120+1]= 5.411115141489; daa[1120+2]= 1.593137043457; daa[11*20+3]= 1.032447924952;
817	daa[11*20+4]= 0.285078800906;
818	daa[1120+5]= 3.945277674515; daa[1120+6]= 2.802427151679; daa[11*20+7]= 0.752042440303;
819	daa[11*20+8]= 1.022507035889;
820	daa[1120+9]= 0.406193586642; daa[1120+10]= 0.445570274261;daa[12*20+0]= 1.253758266664;
821	daa[12*20+1]= 0.983692987457;
822	daa[1220+2]= 0.648441278787; daa[1220+3]= 0.222621897958; daa[12*20+4]= 0.76768882348;
823	daa[12*20+5]= 2.494896077113;
824	daa[1220+6]= 0.55541539747; daa[1220+7]= 0.459436173579; daa[12*20+8]= 0.984311525359;
825	daa[12*20+9]= 3.364797763104;
826	daa[1220+10]= 6.030559379572;daa[1220+11]= 1.073061184332;daa[13*20+0]= 0.492964679748;
827	daa[13*20+1]= 0.371644693209;
828	daa[1320+2]= 0.354861249223; daa[1320+3]= 0.281730694207; daa[13*20+4]= 0.441337471187;
829	daa[13*20+5]= 0.14435695975;
830	daa[1320+6]= 0.291409084165; daa[1320+7]= 0.368166464453; daa[13*20+8]= 0.714533703928;
831	daa[13*20+9]= 1.517359325954;
832	daa[1320+10]= 2.064839703237;daa[1320+11]= 0.266924750511;daa[13*20+12]= 1.77385516883;
833	daa[14*20+0]= 1.173275900924;
834	daa[1420+1]= 0.448133661718; daa[1420+2]= 0.494887043702; daa[14*20+3]= 0.730628272998;
835	daa[14*20+4]= 0.356008498769;
836	daa[1420+5]= 0.858570575674; daa[1420+6]= 0.926563934846; daa[1420+7]= 0.504086599527; daa[1420+8]= 0.527007339151;
837	daa[1420+9]= 0.388355409206; daa[1420+10]= 0.374555687471;daa[1420+11]= 1.047383450722;daa[1420+12]= 0.454123625103;
838	daa[1420+13]= 0.233597909629;daa[1520+0]= 4.325092687057; daa[1520+1]= 1.12278310421; daa[1520+2]= 2.904101656456;
839	daa[1520+3]= 1.582754142065; daa[1520+4]= 1.197188415094; daa[1520+5]= 1.934870924596; daa[1520+6]= 1.769893238937;
840	daa[1520+7]= 1.509326253224; daa[1520+8]= 1.11702976291; daa[1520+9]= 0.35754441246; daa[1520+10]= 0.352969184527;
841	daa[1520+11]= 1.752165917819;daa[1520+12]= 0.918723415746;daa[1520+13]= 0.540027644824;daa[1520+14]= 1.169129577716;
842	daa[1620+0]= 1.729178019485; daa[1620+1]= 0.914665954563; daa[1620+2]= 1.898173634533; daa[1620+3]= 0.934187509431;
843	daa[1620+4]= 1.119831358516; daa[1620+5]= 1.277480294596; daa[1620+6]= 1.071097236007; daa[1620+7]= 0.641436011405;
844	daa[1620+8]= 0.585407090225; daa[1620+9]= 1.17909119726; daa[1620+10]= 0.915259857694;daa[1620+11]= 1.303875200799;
845	daa[1620+12]= 1.488548053722;daa[1620+13]= 0.488206118793;daa[1620+14]= 1.005451683149;daa[1620+15]= 5.15155629227;
846	daa[1720+0]= 0.465839367725; daa[1720+1]= 0.426382310122; daa[1720+2]= 0.191482046247; daa[1720+3]= 0.145345046279;
847	daa[1720+4]= 0.527664418872; daa[1720+5]= 0.758653808642; daa[1720+6]= 0.407635648938; daa[1720+7]= 0.508358924638;
848	daa[1720+8]= 0.30124860078; daa[1720+9]= 0.34198578754; daa[1720+10]= 0.6914746346; daa[1720+11]= 0.332243040634;
849	daa[1720+12]= 0.888101098152;daa[1720+13]= 2.074324893497;daa[1720+14]= 0.252214830027;daa[1720+15]= 0.387925622098;
850	daa[1720+16]= 0.513128126891;daa[1820+0]= 0.718206697586; daa[1820+1]= 0.720517441216; daa[1820+2]= 0.538222519037;
851	daa[1820+3]= 0.261422208965; daa[1820+4]= 0.470237733696; daa[1820+5]= 0.95898974285; daa[1820+6]= 0.596719300346;
852	daa[1820+7]= 0.308055737035; daa[1820+8]= 4.218953969389; daa[1820+9]= 0.674617093228; daa[1820+10]= 0.811245856323;
853	daa[1820+11]= 0.7179934869; daa[1820+12]= 0.951682162246;daa[1820+13]= 6.747260430801;daa[1820+14]= 0.369405319355;
854	daa[1820+15]= 0.796751520761;daa[1820+16]= 0.801010243199;daa[1820+17]= 4.054419006558;daa[1920+0]= 2.187774522005;
855	daa[1920+1]= 0.438388343772; daa[1920+2]= 0.312858797993; daa[1920+3]= 0.258129289418; daa[1920+4]= 1.116352478606;
856	daa[1920+5]= 0.530785790125; daa[1920+6]= 0.524253846338; daa[1920+7]= 0.25334079019; daa[1920+8]= 0.20155597175;
857	daa[1920+9]= 8.311839405458; daa[1920+10]= 2.231405688913;daa[1920+11]= 0.498138475304;daa[1920+12]= 2.575850755315;
858	daa[1920+13]= 0.838119610178;daa[1920+14]= 0.496908410676;daa[1920+15]= 0.561925457442;daa[1920+16]= 2.253074051176;
859	daa[1920+17]= 0.266508731426;daa[1920+18]= 1;
860
861	f[0]= 0.074; f[1]= 0.052; f[2]= 0.045; f[3]= 0.054;
862	f[4]= 0.025; f[5]= 0.034; f[6]= 0.054; f[7]= 0.074;
863	f[8]= 0.026; f[9]= 0.068; f[10]= 0.099; f[11]= 0.058;
864	f[12]= 0.025; f[13]= 0.047; f[14]= 0.039; f[15]= 0.057;
865	f[16]= 0.051; f[17]= 0.013; f[18]= 0.032; f[19]= 0.073;
866	}
867	break;
868	case MTMAM:
869	{
870	daa[120+0]= 32; daa[220+0]= 2; daa[220+1]= 4; daa[320+0]= 11;
871	daa[320+1]= 0; daa[320+2]= 864; daa[420+0]= 0; daa[420+1]= 186;
872	daa[420+2]= 0; daa[420+3]= 0; daa[520+0]= 0; daa[520+1]= 246;
873	daa[520+2]= 8; daa[520+3]= 49; daa[520+4]= 0; daa[620+0]= 0;
874	daa[620+1]= 0; daa[620+2]= 0; daa[620+3]= 569; daa[620+4]= 0;
875	daa[620+5]= 274; daa[720+0]= 78; daa[720+1]= 18; daa[720+2]= 47;
876	daa[720+3]= 79; daa[720+4]= 0; daa[720+5]= 0; daa[720+6]= 22;
877	daa[820+0]= 8; daa[820+1]= 232; daa[820+2]= 458; daa[820+3]= 11;
878	daa[820+4]= 305; daa[820+5]= 550; daa[820+6]= 22; daa[820+7]= 0;
879	daa[920+0]= 75; daa[920+1]= 0; daa[920+2]= 19; daa[920+3]= 0;
880	daa[920+4]= 41; daa[920+5]= 0; daa[920+6]= 0; daa[920+7]= 0;
881	daa[920+8]= 0; daa[1020+0]= 21; daa[1020+1]= 6; daa[1020+2]= 0;
882	daa[1020+3]= 0; daa[1020+4]= 27; daa[1020+5]= 20; daa[1020+6]= 0;
883	daa[1020+7]= 0; daa[1020+8]= 26; daa[1020+9]= 232; daa[1120+0]= 0;
884	daa[1120+1]= 50; daa[1120+2]= 408; daa[1120+3]= 0; daa[1120+4]= 0;
885	daa[1120+5]= 242; daa[1120+6]= 215; daa[1120+7]= 0; daa[1120+8]= 0;
886	daa[1120+9]= 6; daa[1120+10]= 4; daa[1220+0]= 76; daa[1220+1]= 0;
887	daa[1220+2]= 21; daa[1220+3]= 0; daa[1220+4]= 0; daa[1220+5]= 22;
888	daa[1220+6]= 0; daa[1220+7]= 0; daa[1220+8]= 0; daa[1220+9]= 378;
889	daa[1220+10]= 609; daa[1220+11]= 59; daa[1320+0]= 0; daa[1320+1]= 0;
890	daa[1320+2]= 6; daa[1320+3]= 5; daa[1320+4]= 7; daa[1320+5]= 0;
891	daa[1320+6]= 0; daa[1320+7]= 0; daa[1320+8]= 0; daa[1320+9]= 57;
892	daa[1320+10]= 246; daa[1320+11]= 0; daa[1320+12]= 11; daa[1420+0]= 53;
893	daa[1420+1]= 9; daa[1420+2]= 33; daa[1420+3]= 2; daa[1420+4]= 0;
894	daa[1420+5]= 51; daa[1420+6]= 0; daa[1420+7]= 0; daa[1420+8]= 53;
895	daa[1420+9]= 5; daa[1420+10]= 43; daa[1420+11]= 18; daa[1420+12]= 0;
896	daa[1420+13]= 17; daa[1520+0]= 342; daa[1520+1]= 3; daa[1520+2]= 446;
897	daa[1520+3]= 16; daa[1520+4]= 347; daa[1520+5]= 30; daa[1520+6]= 21;
898	daa[1520+7]= 112; daa[1520+8]= 20; daa[1520+9]= 0; daa[1520+10]= 74;
899	daa[1520+11]= 65; daa[1520+12]= 47; daa[1520+13]= 90; daa[1520+14]= 202;
900	daa[1620+0]= 681; daa[1620+1]= 0; daa[1620+2]= 110; daa[1620+3]= 0;
901	daa[1620+4]= 114; daa[1620+5]= 0; daa[1620+6]= 4; daa[1620+7]= 0;
902	daa[1620+8]= 1; daa[1620+9]= 360; daa[1620+10]= 34; daa[1620+11]= 50;
903	daa[1620+12]= 691; daa[1620+13]= 8; daa[1620+14]= 78; daa[1620+15]= 614;
904	daa[1720+0]= 5; daa[1720+1]= 16; daa[1720+2]= 6; daa[1720+3]= 0;
905	daa[1720+4]= 65; daa[1720+5]= 0; daa[1720+6]= 0; daa[1720+7]= 0;
906	daa[1720+8]= 0; daa[1720+9]= 0; daa[1720+10]= 12; daa[1720+11]= 0;
907	daa[1720+12]= 13; daa[1720+13]= 0; daa[1720+14]= 7; daa[1720+15]= 17;
908	daa[1720+16]= 0; daa[1820+0]= 0; daa[1820+1]= 0; daa[1820+2]= 156;
909	daa[1820+3]= 0; daa[1820+4]= 530; daa[1820+5]= 54; daa[1820+6]= 0;
910	daa[1820+7]= 1; daa[1820+8]= 1525;daa[1820+9]= 16; daa[1820+10]= 25;
911	daa[1820+11]= 67; daa[1820+12]= 0; daa[1820+13]= 682; daa[1820+14]= 8;
912	daa[1820+15]= 107; daa[1820+16]= 0; daa[1820+17]= 14; daa[1920+0]= 398;
913	daa[1920+1]= 0; daa[1920+2]= 0; daa[1920+3]= 10; daa[1920+4]= 0;
914	daa[1920+5]= 33; daa[1920+6]= 20; daa[1920+7]= 5; daa[1920+8]= 0;
915	daa[1920+9]= 2220; daa[1920+10]= 100;daa[1920+11]= 0; daa[1920+12]= 832;
916	daa[1920+13]= 6; daa[1920+14]= 0; daa[1920+15]= 0; daa[1920+16]= 237;
917	daa[1920+17]= 0; daa[1920+18]= 0;
918
919	f[0]= 0.06920; f[1]= 0.01840; f[2]= 0.04000; f[3]= 0.018600;
920	f[4]= 0.00650; f[5]= 0.02380; f[6]= 0.02360; f[7]= 0.055700;
921	f[8]= 0.02770; f[9]= 0.09050; f[10]=0.16750; f[11]= 0.02210;
922	f[12]=0.05610; f[13]= 0.06110; f[14]=0.05360; f[15]= 0.07250;
923	f[16]=0.08700; f[17]= 0.02930; f[18]=0.03400; f[19]= 0.04280;
924	}
925	break;
926	case GTR:
927	printf("this function should not be called by GTR model for proteins\n");
928	exit(-1);
929	default:
930	printf("FATAL ERROR not defined\n");
931	exit(-1);
932	}
933	}
934
935
936	/*
937
938	TODO review frequency sums for fixed as well as empirical base frequencies !
939
940	NUMERICAL BUG fix, rounded AA freqs in some models, such that
941	they actually really sum to 1.0 +/- epsilon
942	{
943	double acc = 0.0;
944
945	for(i = 0; i < 20; i++)
946	acc += f[i];
947
948	printf("%1.80f\n", acc);
949	assert(acc == 1.0);
950	}
951	*/
952
953
954	for (i=0; i<20; i++)
955	for (j=0; j<i; j++)
956	daa[j20+i] = daa[i20+j];
957
958	/*
959	for (i=0; i<20; i++)
960	{
961	for (j=0; j<20; j++)
962	{
963	if(i == j)
964	printf("0.0 ");
965	else
966	printf("%f ", daa[i * 20 + j]);
967	}
968	printf("\n");
969	}
970
971	for (i=0; i<20; i++)
972	printf("%f ", f[i]);
973	printf("\n");
974	*/
975
976	max = 0;
977
978	for(i = 0; i < 19; i++)
979	for(j = i + 1; j < 20; j++)
980	{
981	q[i][j] = temp = daa[i * 20 + j];
982	if(temp > max)
983	max = temp;
984	}
985
986	scaler = AA_SCALE / max;
987
988	/* SCALING HAS BEEN RE-INTRODUCED TO RESOLVE NUMERICAL PROBLEMS */
989
990	r = 0;
991	for(i = 0; i < 19; i++)
992	{
993	for(j = i + 1; j < 20; j++)
994	{
995	q[i][j] *= scaler;
996
997	assert(q[i][j] <= AA_SCALE_PLUS_EPSILON);
998
999	initialRates[r++] = q[i][j];
1000	}
1001	}
1002	}
1003
1004	static void updateFracChange(tree *tr)
1005	{
1006	if(tr->NumberOfModels == 1)
1007	{
1008	assert(tr->fracchanges[0] != -1.0);
1009	tr->fracchange = tr->fracchanges[0];
1010	tr->fracchanges[0] = -1.0;
1011	}
1012	else
1013	{
1014	int model, i;
1015	double modelWeights = (double )calloc(tr->NumberOfModels, sizeof(double));
1016	double wgtsum = 0.0;
1017
1018	/assert(tr->rateHetModel == GAMMA \|\| tr->rateHetModel == GAMMA_I);/
1019	assert(tr->NumberOfModels > 1);
1020
1021	tr->fracchange = 0.0;
1022
1023	for(i = 0; i < tr->cdta->endsite; i++)
1024	{
1025	modelWeights[tr->model[i]] += (double)tr->cdta->aliaswgt[i];
1026	wgtsum += (double)tr->cdta->aliaswgt[i];
1027	}
1028
1029	for(model = 0; model < tr->NumberOfModels; model++)
1030	{
1031	/assert(tr->fracchanges[model] != -1.0);/
1032	tr->partitionContributions[model] = modelWeights[model] / wgtsum;
1033	tr->fracchange += tr->partitionContributions[model] * tr->fracchanges[model];
1034	}
1035
1036	free(modelWeights);
1037	}
1038	}
1039
1040
1041	void initReversibleGTR(tree tr, analdef adef, int model)
1042	{
1043	double
1044	*ext_EIGN,
1045	*EV,
1046	*EI,
1047	*frequencies,
1048	*ext_initialRates,
1049	*fracchanges = tr->fracchanges,
1050	*tipVector;
1051
1052	switch(tr->partitionData[model].dataType)
1053	{
1054	case AA_DATA:
1055	ext_EIGN = tr->EIGN_AA;
1056	EV = tr->EV_AA;
1057	EI = tr->EI_AA;
1058	frequencies = tr->frequencies_AA;
1059	ext_initialRates = tr->initialRates_AA;
1060	tipVector = tr->tipVectorAA;
1061
1062	{
1063	double r[20][20], a[20][20], initialRates = &(ext_initialRates[model 190]),
1064	f[20], e[20], d[20];
1065	int i, j, k, m, l;
1066	double invfreq[20], EIGN[20], EIGV[20][20];
1067	double *eptr;
1068
1069	if(adef->useMultipleModel)
1070	{
1071	if(adef->proteinMatrix == GTR)
1072	{
1073	printf("FATAL ERROR no GTR for AA with multiple models\n");
1074	exit(-1);
1075	}
1076	else
1077	{
1078	initProtMat(model, f, tr->partitionData[model].protModels, ext_initialRates, adef->userProteinModel,
1079	adef->externalAAMatrix);
1080
1081	if(tr->partitionData[model].protFreqs)
1082	{
1083	for(l = 0; l < 20; l++)
1084	f[l] = frequencies[model * 20 + l];
1085	}
1086	else
1087	{
1088	for(l = 0; l < 20; l++)
1089	frequencies[model * 20 + l] = f[l];
1090	}
1091	}
1092	}
1093	else
1094	{
1095	if(adef->proteinMatrix == GTR)
1096	{
1097	for(l = 0; l < 20; l++)
1098	f[l] = frequencies[model * 20 + l];
1099	}
1100	else
1101	{
1102	initProtMat(model, f, adef->proteinMatrix, ext_initialRates, adef->userProteinModel,
1103	adef->externalAAMatrix);
1104	if(adef->protEmpiricalFreqs)
1105	{
1106	for(l = 0; l < 20; l++)
1107	f[l] = frequencies[model * 20 + l];
1108	}
1109	else
1110	{
1111	for(l = 0; l < 20; l++)
1112	frequencies[model * 20 + l] = f[l];
1113	}
1114	}
1115	}
1116
1117	i = 0;
1118
1119	for(j = 0; j < 19; j++)
1120	for (k = j+1; k < 20; k++)
1121	r[j][k] = initialRates[i++];
1122
1123	for (j = 0; j < 20; j++)
1124	{
1125	r[j][j] = 0.0;
1126	for (k = 0; k < j; k++)
1127	r[j][k] = r[k][j];
1128	}
1129
1130	fracchanges[model] = 0.0;
1131
1132	for (j = 0; j< 20; j++)
1133	for (k = 0; k< 20; k++)
1134	fracchanges[model] += f[j] * r[j][k] * f[k];
1135
1136
1137	m = 0;
1138
1139	for(i=0; i< 20; i++)
1140	a[i][i] = 0;
1141
1142	for(i=0; i < 20; i++)
1143	{
1144	for(j=i+1; j < 20 ; j++)
1145	{
1146	double factor = initialRates[m++];
1147	a[i][j] = a[j][i] = factor * sqrt( f[i] * f[j]);
1148	a[i][i] -= factor * f[j];
1149	a[j][j] -= factor * f[i];
1150	}
1151	}
1152
1153	tred2((double *)a,20,20,d,e);
1154	tqli(d, e, 20 , 20, (double *)a);
1155
1156	for(i=0; i<20; i++)
1157	for(j=0; j<20; j++)
1158	a[i][j] *= sqrt(f[j]);
1159
1160	for (i=0; i<20; i++)
1161	{
1162	if (d[i] > -1e-8)
1163	{
1164	if (i != 0)
1165	{
1166	double tmp = d[i], sum=0;
1167	d[i] = d[0];
1168	d[0] = tmp;
1169	for (j=0; j < 20; j++)
1170	{
1171	tmp = a[i][j];
1172	a[i][j] = a[0][j];
1173	sum += (a[0][j] = tmp);
1174	}
1175	for (j=0; j < 20; j++)
1176	a[0][j] /= sum;
1177	}
1178	break;
1179	}
1180	}
1181
1182	for (i=0; i< 20; i++)
1183	{
1184	EIGN[i] = -d[i];
1185
1186	for (j=0; j<20; j++)
1187	EIGV[i][j] = a[j][i];
1188	invfreq[i] = 1 / EIGV[i][0];
1189	}
1190
1191	for(l = 1; l < 20; l++)
1192	ext_EIGN[model * 19 + (l - 1)] = EIGN[l];
1193
1194	eptr = &(EV[model * 400]);
1195
1196	for(i = 0; i < 20; i++)
1197	for(j = 0; j < 20; j++)
1198	*eptr++ = EIGV[i][j];
1199
1200	for(i = 0; i < 20; i++)
1201	for(j = 1; j < 20; j++)
1202	EI[model * 380 + i * 19 + (j - 1)] = EV[model * 400 + i * 20 + j] * invfreq[i];
1203
1204	for(i=0; i < 23; i++)
1205	{
1206	for(j = 0; j < 20; j++)
1207	{
1208	tipVector[model * 460 + 20 * i + j] = 0.0;
1209	}
1210
1211	if(i < 20)
1212	{
1213	for (j = 0; j < 20; j++)
1214	{
1215	tipVector[model * 460 + 20 * i + j] += EIGV[i][j];
1216	}
1217	}
1218	else
1219	{
1220	if(i == 20)
1221	{
1222	for (j = 0; j < 20; j++)
1223	{
1224	tipVector[model * 460 + 20 * i + j] += EIGV[2][j] + EIGV[3][j];
1225	}
1226	}
1227	else
1228	{
1229	if(i == 21)
1230	{
1231	for (j = 0; j < 20; j++)
1232	{
1233	tipVector[model * 460 + 20 * i + j] += EIGV[5][j] + EIGV[6][j];
1234	}
1235	}
1236	else
1237	{
1238	if(i == 22)
1239	{
1240	for (j = 0; j < 20; j++)
1241	{
1242	for(k = 0; k < 20; k++)
1243	{
1244	tipVector[model * 460 + 20 * i + j] += EIGV[k][j];
1245	}
1246	}
1247	}
1248	}
1249	}
1250	}
1251	}
1252	}
1253	break;
1254	case DNA_DATA:
1255	ext_EIGN = tr->EIGN_DNA;
1256	EV = tr->EV_DNA;
1257	EI = tr->EI_DNA;
1258	frequencies = tr->frequencies_DNA;
1259	ext_initialRates = tr->initialRates_DNA;
1260	tipVector = tr->tipVectorDNA;
1261
1262	{
1263	double r[4][4], a[4][4], initialRates = &(ext_initialRates[model 5]),
1264	f[4], e[4], d[4];
1265	int i, j, k, m, code;
1266	double invfreq[4], EIGN[4], EIGV[4][4];
1267	double *eptr;
1268
1269	f[0] = frequencies[model * 4];
1270	f[1] = frequencies[model * 4 + 1];
1271	f[2] = frequencies[model * 4 + 2];
1272	f[3] = frequencies[model * 4 + 3];
1273
1274	i = 0;
1275	for (j = 0; j < 2; j++)
1276	for (k = j+1; k< 4; k++)
1277	r[j][k] = initialRates[i++];
1278	r[2][3] = 1.0;
1279	for (j = 0; j < 4; j++)
1280	{
1281	r[j][j] = 0;
1282	for (k = 0; k < j; k++)
1283	r[j][k] = r[k][j];
1284	}
1285
1286	fracchanges[model] = 0.0;
1287
1288	for (j = 0; j<4; j++)
1289	for (k = 0; k<4; k++)
1290	fracchanges[model] += f[j] * r[j][k] * f[k];
1291
1292	m = 0;
1293
1294	for(i=0; i<4; i++)
1295	a[i][i] = 0;
1296
1297	for(i=0; i<4; i++)
1298	{
1299	for(j=i+1; j<4; j++)
1300	{
1301	double factor = ((m>=5) ? 1.0 : initialRates[m++]);
1302	a[i][j] = a[j][i] = factor * sqrt(f[i] * f[j]);
1303	a[i][i] -= factor * f[j];
1304	a[j][j] -= factor * f[i];
1305	}
1306	}
1307
1308	tred2((double *)a,4,4,d,e);
1309	tqli(d, e, 4 , 4, (double *)a);
1310
1311	for(i=0; i<4; i++)
1312	for(j=0; j<4; j++)
1313	a[i][j] *= sqrt(f[j]);
1314
1315	for (i=0; i<4; i++)
1316	{
1317	if (d[i] > -1e-8)
1318	{
1319	if (i != 0)
1320	{
1321	double tmp = d[i], sum=0;
1322	d[i] = d[0];
1323	d[0] = tmp;
1324	for (j=0; j < 4; j++)
1325	{
1326	tmp = a[i][j];
1327	a[i][j] = a[0][j];
1328	sum += (a[0][j] = tmp);
1329	}
1330	for (j=0; j < 4; j++)
1331	a[0][j] /= sum;
1332	}
1333	break;
1334	}
1335	}
1336
1337	for (i=0; i< 4; i++)
1338	{
1339	EIGN[i] = -d[i];
1340
1341	for (j=0; j<4; j++)
1342	EIGV[i][j] = a[j][i];
1343	invfreq[i] = 1 / EIGV[i][0];
1344	}
1345
1346	ext_EIGN[model * 3] = EIGN[1];
1347	ext_EIGN[model * 3 + 1] = EIGN[2];
1348	ext_EIGN[model * 3 + 2] = EIGN[3];
1349
1350	eptr = &(EV[model * 16]);
1351
1352	for(i = 0; i < 4; i++)
1353	for(j = 0; j < 4; j++)
1354	*eptr++ = EIGV[i][j];
1355
1356	EI[model * 12] = EV[model * 16 + 1] * invfreq[0];
1357	EI[model * 12 + 1] = EV[model * 16 + 2] * invfreq[0];
1358	EI[model * 12 + 2] = EV[model * 16 + 3] * invfreq[0];
1359
1360	EI[model * 12 + 3] = EV[model * 16 + 5] * invfreq[1];
1361	EI[model * 12 + 4] = EV[model * 16 + 6] * invfreq[1];
1362	EI[model * 12 + 5] = EV[model * 16 + 7] * invfreq[1];
1363
1364	EI[model * 12 + 6] = EV[model * 16 + 9] * invfreq[2];
1365	EI[model * 12 + 7] = EV[model * 16 + 10] * invfreq[2];
1366	EI[model * 12 + 8] = EV[model * 16 + 11] * invfreq[2];
1367
1368	EI[model * 12 + 9] = EV[model * 16 + 13] * invfreq[3];
1369	EI[model * 12 + 10] = EV[model * 16 + 14] * invfreq[3];
1370	EI[model * 12 + 11] = EV[model * 16 + 15] * invfreq[3];
1371
1372	for(i=0; i < 16; i++)
1373	{
1374	code = i;
1375
1376	tipVector[model * 64 + i * 4] = 0;
1377	tipVector[model * 64 + i * 4 + 1] = 0;
1378	tipVector[model * 64 + i * 4 + 2] = 0;
1379	tipVector[model * 64 + i * 4 + 3] = 0;
1380
1381	if(i > 0)
1382	{
1383	for (j = 0; j < 4; j++)
1384	{
1385	if ((code >> j) & 1)
1386	{
1387	int jj = "0123"[j] - '0';
1388	tipVector[model * 64 + i * 4] += EIGV[jj][0];
1389	tipVector[model * 64 + i * 4 + 1] += EIGV[jj][1];
1390	tipVector[model * 64 + i * 4 + 2] += EIGV[jj][2];
1391	tipVector[model * 64 + i * 4 + 3] += EIGV[jj][3];
1392	}
1393	}
1394	}
1395	}
1396	}
1397	break;
1398	default:
1399	assert(0);
1400	}
1401
1402
1403	updateFracChange(tr);
1404	}
1405
1406
1407	double LnGamma (double alpha)
1408	{
1409	/* returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.
1410	Stirling's formula is used for the central polynomial part of the procedure.
1411	Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
1412	Communications of the Association for Computing Machinery, 9:684
1413	*/
1414	double x, f, z, result;
1415
1416	x = alpha;
1417	f = 0.0;
1418
1419	if ( x < 7.0)
1420	{
1421	f = 1.0;
1422	z = alpha - 1.0;
1423
1424	while ((z = z + 1.0) < 7.0)
1425	{
1426	f *= z;
1427	}
1428	x = z;
1429
1430	assert(f != 0.0);
1431
1432	f=-log(f);
1433	}
1434
1435	z = 1/(x*x);
1436
1437	result = f + (x-0.5)*log(x) - x + .918938533204673
1438	+ (((-.000595238095238z+.000793650793651)z-.002777777777778)*z
1439	+.083333333333333)/x;
1440
1441	return result;
1442	}
1443
1444
1445
1446	double IncompleteGamma (double x, double alpha, double ln_gamma_alpha)
1447	{
1448	/* returns the incomplete gamma ratio I(x,alpha) where x is the upper
1449	limit of the integration and alpha is the shape parameter.
1450	returns (-1) if in error
1451	ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
1452	(1) series expansion if (alpha>x \|\| x<=1)
1453	(2) continued fraction otherwise
1454	RATNEST FORTRAN by
1455	Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics,
1456	19: 285-287 (AS32)
1457	*/
1458	int i;
1459	double p=alpha, g=ln_gamma_alpha;
1460	double accurate=1e-8, overflow=1e30;
1461	double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, pn[6];
1462
1463
1464	if (x==0) return (0);
1465	if (x<0 \|\| p<=0) return (-1);
1466
1467
1468	factor=exp(p*log(x)-x-g);
1469	if (x>1 && x>=p) goto l30;
1470	/* (1) series expansion */
1471	gin=1; term=1; rn=p;
1472	l20:
1473	rn++;
1474	term*=x/rn; gin+=term;
1475
1476	if (term > accurate) goto l20;
1477	gin*=factor/p;
1478	goto l50;
1479	l30:
1480	/* (2) continued fraction */
1481	a=1-p; b=a+x+1; term=0;
1482	pn[0]=1; pn[1]=x; pn[2]=x+1; pn[3]=x*b;
1483	gin=pn[2]/pn[3];
1484	l32:
1485	a++;
1486	b+=2;
1487	term++;
1488	an=a*term;
1489	for (i=0; i<2; i++)
1490	pn[i+4]=bpn[i+2]-anpn[i];
1491	if (pn[5] == 0) goto l35;
1492	rn=pn[4]/pn[5];
1493	dif=fabs(gin-rn);
1494	if (dif>accurate) goto l34;
1495	if (dif<=accurate*rn) goto l42;
1496	l34:
1497	gin=rn;
1498	l35:
1499	for (i=0; i<4; i++)
1500	pn[i]=pn[i+2];
1501	if (fabs(pn[4]) < overflow)
1502	goto l32;
1503
1504	for (i=0; i<4; i++)
1505	pn[i]/=overflow;
1506
1507
1508	goto l32;
1509	l42:
1510	gin=1-factor*gin;
1511
1512	l50:
1513	return (gin);
1514	}
1515
1516
1517
1518
1519	double PointNormal (double prob)
1520	{
1521	/* returns z so that Prob{x<z}=prob where x ~ N(0,1) and (1e-12)<prob<1-(1e-12)
1522	returns (-9999) if in error
1523	Odeh RE & Evans JO (1974) The percentage points of the normal distribution.
1524	Applied Statistics 22: 96-97 (AS70)
1525
1526	Newer methods:
1527	Wichura MJ (1988) Algorithm AS 241: the percentage points of the
1528	normal distribution. 37: 477-484.
1529	Beasley JD & Springer SG (1977). Algorithm AS 111: the percentage
1530	points of the normal distribution. 26: 118-121.
1531
1532	*/
1533	double a0=-.322232431088, a1=-1, a2=-.342242088547, a3=-.0204231210245;
1534	double a4=-.453642210148e-4, b0=.0993484626060, b1=.588581570495;
1535	double b2=.531103462366, b3=.103537752850, b4=.0038560700634;
1536	double y, z=0, p=prob, p1;
1537
1538	p1 = (p<0.5 ? p : 1-p);
1539	if (p1<1e-20) return (-9999);
1540
1541	y = sqrt (log(1/(p1*p1)));
1542	z = y + ((((ya4+a3)y+a2)y+a1)y+a0) / ((((yb4+b3)y+b2)y+b1)y+b0);
1543	return (p<0.5 ? -z : z);
1544	}
1545
1546
1547	double PointChi2 (double prob, double v)
1548	{
1549	/* returns z so that Prob{x<z}=prob where x is Chi2 distributed with df=v
1550	returns -1 if in error. 0.000002<prob<0.999998
1551	RATNEST FORTRAN by
1552	Best DJ & Roberts DE (1975) The percentage points of the
1553	Chi2 distribution. Applied Statistics 24: 385-388. (AS91)
1554	Converted into C by Ziheng Yang, Oct. 1993.
1555	*/
1556	double e=.5e-6, aa=.6931471805, p=prob, g;
1557	double xx, c, ch, a=0,q=0,p1=0,p2=0,t=0,x=0,b=0,s1,s2,s3,s4,s5,s6;
1558
1559	if (p<.000002 \|\| p>.999998 \|\| v<=0) return (-1);
1560
1561	g = LnGamma(v/2);
1562
1563	xx=v/2; c=xx-1;
1564	if (v >= -1.24*log(p)) goto l1;
1565
1566	ch=pow((pxxexp(g+xx*aa)), 1/xx);
1567	if (ch-e<0) return (ch);
1568	goto l4;
1569	l1:
1570	if (v>.32) goto l3;
1571	ch=0.4; a=log(1-p);
1572	l2:
1573	q=ch; p1=1+ch(4.67+ch); p2=ch(6.73+ch*(6.66+ch));
1574	t=-0.5+(4.67+2ch)/p1 - (6.73+ch(13.32+3*ch))/p2;
1575	ch-=(1-exp(a+g+.5ch+caa)*p2/p1)/t;
1576	if (fabs(q/ch-1)-.01 <= 0) goto l4;
1577	else goto l2;
1578
1579	l3:
1580	x=PointNormal (p);
1581	p1=0.222222/v; ch=vpow((xsqrt(p1)+1-p1), 3.0);
1582	if (ch>2.2v+6) ch=-2(log(1-p)-clog(.5ch)+g);
1583	l4:
1584	q=ch; p1=.5*ch;
1585	if ((t=IncompleteGamma (p1, xx, g))< 0.0)
1586	{
1587	printf ("IncompleteGamma \n");
1588	return (-1);
1589	}
1590
1591	p2=p-t;
1592	t=p2exp(xxaa+g+p1-c*log(ch));
1593	b=t/ch; a=0.5t-bc;
1594
1595	s1=(210+a(140+a(105+a(84+a(70+60*a))))) / 420;
1596	s2=(420+a(735+a(966+a(1141+1278a))))/2520;
1597	s3=(210+a(462+a(707+932*a)))/2520;
1598	s4=(252+a(672+1182a)+c(294+a(889+1740*a)))/5040;
1599	s5=(84+264a+c(175+606*a))/2520;
1600	s6=(120+c(346+127c))/5040;
1601	ch+=t(1+0.5ts1-bc(s1-b(s2-b(s3-b(s4-b(s5-bs6))))));
1602	if (fabs(q/ch-1) > e) goto l4;
1603
1604	return (ch);
1605	}
1606
1607
1608
1609
1610
1611
1612	void makeGammaCats(int model, double alphas, double gammaRates)
1613	{
1614	int i, K = 4;
1615	double factor, lnga1, alfa, beta;
1616	double gammaProbs = (double )malloc(K * sizeof(double));
1617
1618	alfa = beta = alphas[model];
1619
1620	/* Note that ALPHA_MIN setting is somewhat critical due to */
1621	/* numerical instability caused by very small rate[0] values */
1622	/* induced by low alpha values around 0.01 */
1623
1624	assert(alfa >= ALPHA_MIN);
1625
1626	factor = alfa / beta * K;
1627
1628	lnga1=LnGamma(alfa+1);
1629
1630	for (i=0; i<K-1; i++)
1631	gammaProbs[i]=PointGamma((i+1.0)/K, alfa, beta);
1632
1633	for (i=0; i<K-1; i++)
1634	gammaProbs[i]=IncompleteGamma(gammaProbs[i]*beta, alfa+1, lnga1);
1635
1636	gammaRates[model * K] = gammaProbs[0]*factor;
1637
1638	gammaRates[model * K + (K-1)] = (1 - gammaProbs[K-2])*factor;
1639
1640	for (i=1; i<K-1; i++)
1641	gammaRates[model * K + i] = (gammaProbs[i] - gammaProbs[i-1])*factor;
1642
1643	free(gammaProbs);
1644
1645	return;
1646	}
1647
1648
1649	void assignLikelihoodFunctions(tree tr, analdef adef)
1650	{
1651	switch(adef->model)
1652	{
1653	case M_PROTGAMMA:
1654	if(adef->useMultipleModel)
1655	{
1656	if(adef->useInvariant)
1657	tr->likelihoodFunction = PROTGAMMAMULTI;
1658	else
1659	tr->likelihoodFunction = PROTGAMMAMULT;
1660	}
1661	else
1662	{
1663	if(adef->useInvariant)
1664	tr->likelihoodFunction = PROTGAMMAI;
1665	else
1666	tr->likelihoodFunction = PROTGAMMA;
1667	}
1668	break;
1669	case M_PROTCAT:
1670	if(adef->useMultipleModel)
1671	tr->likelihoodFunction = PROTCATMULT;
1672	else
1673	tr->likelihoodFunction = PROTCAT;
1674	break;
1675	case M_GTRGAMMA:
1676	if(adef->useMultipleModel)
1677	{
1678	if(adef->useInvariant)
1679	tr->likelihoodFunction = GTRGAMMAMULTI;
1680	else
1681	tr->likelihoodFunction = GTRGAMMAMULT;
1682	}
1683	else
1684	{
1685	if(adef->useInvariant)
1686	tr->likelihoodFunction = GTRGAMMAI;
1687	else
1688	tr->likelihoodFunction = GTRGAMMA;
1689	}
1690	break;
1691	case M_GTRCAT:
1692	if(adef->useMultipleModel)
1693	tr->likelihoodFunction = GTRCATMULT;
1694	else
1695	tr->likelihoodFunction = GTRCAT;
1696	break;
1697	default:
1698	printf("FATAL ERROR: assignLikelihoodFunctions\n");
1699	exit(1);
1700	break;
1701	}
1702	}
1703
1704
1705	static void initInvariant(tree *tr, int lower, int upper,
1706	int numberOfInvariableColumns, int weightOfInvariableColumns,
1707	int dataType)
1708	{
1709	int count = 0, sum = 0, i, j;
1710
1711
1712	switch(dataType)
1713	{
1714	case AA_DATA:
1715	for(i = lower; i < upper; i++)
1716	{
1717	char c[23];
1718	int aaSum;
1719
1720	for(j = 0; j < 23; j++)
1721	c[j] = 0;
1722
1723	for(j = 1; j <= tr->mxtips; j++)
1724	c[tr->yVector[tr->nodep[j]->number][i]] = 1;
1725
1726	aaSum = 0;
1727
1728	for(j = 0; j < 20; j++)
1729	aaSum += c[j];
1730
1731	if(aaSum == 1)
1732	{
1733	if(((c[20] == 1) && ((c[2] == 0) \|\| (c[3] == 0))) \|\| ((c[21] == 1) && ((c[5] == 0) \|\| (c[6] == 0))))
1734	tr->invariant[i] = 20;
1735	else
1736	{
1737	int which = -1;
1738
1739	for(j = 0; (j < 20) && (which < 0); j++)
1740	if(c[j] == 1)
1741	which = j;
1742
1743	tr->invariant[i] = which;
1744	count++;
1745	sum += tr->cdta->aliaswgt[i];
1746	}
1747	}
1748	else
1749	tr->invariant[i] = 20;
1750	}
1751	break;
1752	case DNA_DATA:
1753	for(i = lower; i < upper; i++)
1754	{
1755	char c[16];
1756
1757	for(j = 0; j < 16; j++)
1758	c[j] = 0;
1759
1760	for(j = 1; j <= tr->mxtips; j++)
1761	c[tr->yVector[tr->nodep[j]->number][i]] = 1;
1762
1763
1764	if(c[1] + c[2] + c[4] + c[8] == 1)
1765	{
1766	if((c[1] && !(c[6] \|\| c[10] \|\| c[12] \|\| c[14])) \|\|
1767	(c[2] && !(c[5] \|\| c[9] \|\| c[12] \|\| c[13])) \|\|
1768	(c[4] && !(c[3] \|\| c[9] \|\| c[10] \|\| c[11])) \|\|
1769	(c[8] && !(c[3] \|\| c[5] \|\| c[6] \|\| c[7])))
1770	{
1771	count++;
1772	sum += tr->cdta->aliaswgt[i];
1773	if(c[1]) tr->invariant[i] = 0;
1774	if(c[2]) tr->invariant[i] = 1;
1775	if(c[4]) tr->invariant[i] = 2;
1776	if(c[8]) tr->invariant[i] = 3;
1777	}
1778	else
1779	tr->invariant[i] = 4;
1780	}
1781	else
1782	tr->invariant[i] = 4;
1783	}
1784	break;
1785	default:
1786	assert(0);
1787	}
1788
1789	*numberOfInvariableColumns += count;
1790	*weightOfInvariableColumns += sum;
1791	}
1792
1793	void initModel(tree tr, rawdata rdta, cruncheddata cdta, analdef adef)
1794	{
1795	int model, i, j;
1796	double temp, wtemp;
1797
1798	optimizeRatesInvocations = 1;
1799	optimizeRateCategoryInvocations = 1;
1800	optimizeAlphaInvocations = 1;
1801	optimizeInvarInvocations = 1;
1802
1803	tr->numberOfInvariableColumns = 0;
1804	tr->weightOfInvariableColumns = 0;
1805
1806
1807	if(!adef->useMultipleModel)
1808	{
1809	assert(tr->partitionData[0].lower == 0);
1810	assert(tr->partitionData[0].upper == tr->cdta->endsite);
1811	}
1812
1813	if(adef->useInvariant)
1814	{
1815	for(model = 0; model < tr->NumberOfModels; model++)
1816	{
1817	int lower = tr->partitionData[model].lower;
1818	int upper = tr->partitionData[model].upper;
1819
1820	initInvariant(tr, lower, upper,
1821	&(tr->numberOfInvariableColumns),
1822	&(tr->weightOfInvariableColumns),
1823	tr->partitionData[model].dataType);
1824	}
1825	}
1826
1827	for(model = 0; model < tr->NumberOfModels; model++)
1828	{
1829	for(i = 0; i < DNA_RATES; i++)
1830	tr->initialRates_DNA[model * DNA_RATES + i] = 0.5;
1831	for(i = 0; i < AA_RATES; i++)
1832	tr->initialRates_AA[model * AA_RATES + i] = 0.5;
1833
1834	if(adef->useInvariant)
1835	{
1836	int lower, upper;
1837	int count = 0;
1838	int total = 0;
1839
1840	lower = tr->partitionData[model].lower;
1841	upper = tr->partitionData[model].upper;
1842
1843	for(i = lower; i < upper; i++)
1844	{
1845	if(tr->invariant[i] < 4)
1846	count += tr->cdta->aliaswgt[i];
1847	total += tr->cdta->aliaswgt[i];
1848	}
1849	tr->invariants[model] = ((double)count)/((double) total);
1850	}
1851	}
1852
1853	assignLikelihoodFunctions(tr, adef);
1854
1855	tr->NumberOfCategories = 1;
1856
1857	switch(adef->model)
1858	{
1859	case M_PROTGAMMA:
1860	case M_GTRGAMMA:
1861	for(j = 0; j < tr->cdta->endsite; j++)
1862	tr->cdta->wr[j] = tr->cdta->aliaswgt[j];
1863	break;
1864	case M_PROTCAT:
1865	case M_GTRCAT:
1866	for (j = 0; j < tr->cdta->endsite; j++)
1867	{
1868	tr->cdta->patrat[j] = temp = 1.0;
1869	tr->cdta->patratStored[j] = 1.0;
1870	tr->cdta->rateCategory[j] = 0;
1871	tr->cdta->wr[j] = wtemp = temp * tr->cdta->aliaswgt[j];
1872	tr->cdta->wr2[j] = temp * wtemp;
1873	}
1874	break;
1875	default:
1876	assert(0);
1877	}
1878
1879	baseFrequenciesGTR(rdta, cdta, tr, adef);
1880
1881	for(model = 0; model < tr->NumberOfModels; model++)
1882	{
1883	tr->alphas[model] = 1.0;
1884	initReversibleGTR(tr, adef, model);
1885	makeGammaCats(model, tr->alphas, tr->gammaRates);
1886	}
1887
1888
1889	if(tr->NumberOfModels > 1)
1890	{
1891	tr->fracchange = 0;
1892	for(model = 0; model < tr->NumberOfModels; model++)
1893	{
1894	tr->fracchange += tr->fracchanges[model];
1895	}
1896	tr->fracchange /= ((double)tr->NumberOfModels);
1897	}
1898
1899	#ifdef _LOCAL_DATA
1900	masterBarrier(THREAD_COPY_INIT_MODEL, tr);
1901	#endif
1902
1903	}
1904
1905
1906
1907

Note: See TracBrowser for help on using the repository browser.

Download in other formats: