1codingMapperVector =
2{
3{                 4,                 9,                 7,                10,                17,                15,                12,                16,                 0,                19,                 6,                13,                11,                 8,                 2,                 3,                 1,                18,                14,                 5}
4};
5
6function checkPair (c1,c2)
7{
8	if (c1!=c2)
9	{
10		c1 = _Genetic_Code[c1];
11		c2 = _Genetic_Code[c2];
12		if (c1 == 10 || c2 == 10)
13		{
14			return 0;
15		}
16
17		if (c1>=10)
18		{
19			c1 = c1-1;
20		}
21		if (c2>=10)
22		{
23			c2 = c2-1;
24		}
25
26		c1 = codingMapperVector[c1];
27		c2 = codingMapperVector[c2];
28
29		allowedMoves[c1][c2] = 1;
30		allowedMoves[c2][c1] = 1;
31	}
32	return 0;
33}
34
35/*---------------------------------------------------------------------------------------------------------------------------------------------*/
36
37#include "chooseGeneticCode.def";
38
39allowedMoves = {20,20};
40
41for (h=0; h<64; h=h+1)
42{
43	if (_Genetic_Code[h] != 10)
44	{
45		pos12 = h-h%4;
46		pos13 = h-4*(h%16$4);
47		pos23 = h-16*(h$16);
48
49		checkPair (h,pos12);
50		checkPair (h,pos12+1);
51		checkPair (h,pos12+2);
52		checkPair (h,pos12+3);
53		checkPair (h,pos13);
54		checkPair (h,pos13+4);
55		checkPair (h,pos13+8);
56		checkPair (h,pos13+12);
57		checkPair (h,pos23);
58		checkPair (h,pos23+16);
59		checkPair (h,pos23+32);
60		checkPair (h,pos23+48);
61	}
62}
63
64#include "modelParameters2.mdl";
65
66aaNames = {{"Alanine",
67"Cysteine",
68"Aspartic_Acid",
69"Glutamic_Acid",
70"Phenylalanine",
71"Glycine",
72"Histidine",
73"Isoleucine",
74"Lysine",
75"Leucine",
76"Methionine",
77"Asparagine",
78"Proline",
79"Glutamine",
80"Arginine",
81"Serine",
82"Threonine",
83"Valine",
84"Tryptophan",
85"Tyrosine"}};
86
87
88if (!SKIP_MODEL_PARAMETER_LIST)
89{
90	if (modelType == 1)
91	{
92		#include "defineGamma.mdl";
93	}
94	if (modelType == 2)
95	{
96		#include "defineHM.mdl";
97	}
98}
99
100function PopulateModelMatrix (ModelMatrixName&, EFV)
101{
102
103	ModelMatrixName = {20,20};
104	if (categoriesUsed)
105	{
106		for (k=0; k<19; k=k+1)
107		{
108			for (k2 = k+1; k2 < 20; k2=k2+1)
109			{
110				if (allowedMoves[k][k2])
111				{
112					subString = aaNames[k]+aaNames[k2];
113				}
114				else
115				{
116					subString = "MultiStep";
117				}
118				ExecuteCommands("global "+subString +"=1; ModelMatrixName[k][k2]:= c*t*"+subString+
119								";ModelMatrixName[k2][k]:= c*t*"+subString+";");
120			}
121		}
122	}
123	else
124	{
125		for (k=0; k<19; k=k+1)
126		{
127			for (k2 = k+1; k2 < 20; k2=k2+1)
128			{
129				if (allowedMoves[k][k2])
130				{
131					subString = aaNames[k]+aaNames[k2];
132				}
133				else
134				{
135					subString = "MultiStep";
136				}
137				ExecuteCommands("global "+subString +"=1; ModelMatrixName[k][k2]:= t*"+subString+
138								";ModelMatrixName[k2][k]:= t*"+subString+";");
139			}
140		}
141	}
142
143	subString = aaNames[7]+aaNames[9];
144	ExecuteCommands(subString +":=1;");
145	return 1;
146}
147
148
149
150redREVMatrix = 0;
151
152HarvestFrequencies (vectorOfFrequencies,filteredData,1,1,0);
153
154MULTIPLY_BY_FREQS = PopulateModelMatrix ("redREVMatrix",vectorOfFrequencies);
155
156Model redREVModel = (redREVMatrix, vectorOfFrequencies, MULTIPLY_BY_FREQS);
157
158FREQUENCY_SENSITIVE = 1;
159