1 //
2 // C++ Implementation: substmodel
3 //
4 // Description:
5 //
6 //
7 // Author: BUI Quang Minh, Steffen Klaere, Arndt von Haeseler <minh.bui@univie.ac.at>, (C) 2008
8 //
9 // Copyright: See COPYING file that comes with this distribution
10 //
11 //
12 #include "modelsubst.h"
13 #include "utils/tools.h"
14
ModelSubst(int nstates)15 ModelSubst::ModelSubst(int nstates) : Optimization(), CheckpointFactory()
16 {
17 num_states = nstates;
18 name = "JC";
19 full_name = "JC (Juke and Cantor, 1969)";
20 state_freq = new double[num_states];
21 for (int i = 0; i < num_states; i++)
22 state_freq[i] = 1.0 / num_states;
23 freq_type = FREQ_EQUAL;
24 fixed_parameters = false;
25 // linked_model = NULL;
26 }
27
startCheckpoint()28 void ModelSubst::startCheckpoint() {
29 checkpoint->startStruct("ModelSubst");
30 }
31
saveCheckpoint()32 void ModelSubst::saveCheckpoint() {
33 startCheckpoint();
34 // CKP_SAVE(num_states);
35 // CKP_SAVE(name);
36 // CKP_SAVE(full_name);
37 // CKP_SAVE(freq_type);
38 if (freq_type == FREQ_ESTIMATE && !fixed_parameters)
39 CKP_ARRAY_SAVE(num_states, state_freq);
40 endCheckpoint();
41 CheckpointFactory::saveCheckpoint();
42 }
43
restoreCheckpoint()44 void ModelSubst::restoreCheckpoint() {
45 CheckpointFactory::restoreCheckpoint();
46 startCheckpoint();
47 // CKP_RESTORE(num_states);
48 // CKP_RESTORE(name);
49 // CKP_RESTORE(full_name);
50 // int freq_type = this->freq_type;
51 // CKP_RESTORE(freq_type);
52 // this->freq_type = (StateFreqType)freq_type;
53 if (freq_type == FREQ_ESTIMATE && !fixed_parameters)
54 CKP_ARRAY_RESTORE(num_states, state_freq);
55 endCheckpoint();
56
57 decomposeRateMatrix();
58 }
59
60 // here the simplest Juke-Cantor model is implemented, valid for all kind of data (DNA, AA,...)
computeTransMatrix(double time,double * trans_matrix,int mixture)61 void ModelSubst::computeTransMatrix(double time, double *trans_matrix, int mixture) {
62 double non_diagonal = (1.0 - exp(-time*num_states/(num_states - 1))) / num_states;
63 double diagonal = 1.0 - non_diagonal * (num_states - 1);
64 int nstates_sqr = num_states * num_states;
65
66 for (int i = 0; i < nstates_sqr; i++)
67 if (i % (num_states+1) == 0)
68 trans_matrix[i] = diagonal;
69 else
70 trans_matrix[i] = non_diagonal;
71 }
72
73
computeTrans(double time,int state1,int state2)74 double ModelSubst::computeTrans(double time, int state1, int state2) {
75 double expt = exp(-time * num_states / (num_states-1));
76 if (state1 != state2) {
77 return (1.0 - expt) / num_states;
78 }
79 return (1.0 + (num_states-1)*expt) / num_states;
80
81 /* double non_diagonal = (1.0 - exp(-time*num_states/(num_states - 1))) / num_states;
82 if (state1 != state2)
83 return non_diagonal;
84 return 1.0 - non_diagonal * (num_states - 1);*/
85 }
86
computeTrans(double time,int model_id,int state1,int state2)87 double ModelSubst::computeTrans(double time, int model_id, int state1, int state2) {
88 return computeTrans(time, state1, state2);
89 }
90
computeTrans(double time,int state1,int state2,double & derv1,double & derv2)91 double ModelSubst::computeTrans(double time, int state1, int state2, double &derv1, double &derv2) {
92 double coef = -double(num_states) / (num_states-1);
93 double expt = exp(time * coef);
94 if (state1 != state2) {
95 derv1 = expt / (num_states-1);
96 derv2 = derv1 * coef;
97 return (1.0 - expt) / num_states;
98 }
99
100 derv1 = -expt;
101 derv2 = derv1 * coef;
102 return (1.0 + (num_states-1)*expt) / num_states;
103 }
104
computeTrans(double time,int model_id,int state1,int state2,double & derv1,double & derv2)105 double ModelSubst::computeTrans(double time, int model_id, int state1, int state2, double &derv1, double &derv2) {
106 return computeTrans(time, state1, state2, derv1, derv2);
107 }
108
getRateMatrix(double * rate_mat)109 void ModelSubst::getRateMatrix(double *rate_mat) {
110 int nrate = getNumRateEntries();
111 for (int i = 0; i < nrate; i++)
112 rate_mat[i] = 1.0;
113 }
114
getQMatrix(double * q_mat)115 void ModelSubst::getQMatrix(double *q_mat) {
116 int i, j, k;
117 for (i = 0, k = 0; i < num_states; i++)
118 for (j = 0; j < num_states; j++, k++)
119 if (i == j) q_mat[k] = -1.0; else q_mat[k] = 1.0/3;
120 }
121
getStateFrequency(double * state_freq,int mixture)122 void ModelSubst::getStateFrequency(double *state_freq, int mixture) {
123 double freq = 1.0 / num_states;
124 for (int i = 0; i < num_states; i++)
125 state_freq[i] = freq;
126 }
127
setStateFrequency(double * state_freq)128 void ModelSubst::setStateFrequency(double *state_freq) {
129 memcpy(this->state_freq, state_freq, sizeof(double)*num_states);
130 }
131
computeTransDerv(double time,double * trans_matrix,double * trans_derv1,double * trans_derv2,int mixture)132 void ModelSubst::computeTransDerv(double time, double *trans_matrix,
133 double *trans_derv1, double *trans_derv2, int mixture)
134 {
135 double expf = exp(-time*num_states/(num_states - 1));
136 double non_diag = (1.0 - expf) / num_states;
137 double diag = 1.0 - non_diag * (num_states - 1);
138 double derv1_non_diag = expf / (num_states-1);
139 double derv1_diag = -expf;
140 double derv2_non_diag = -derv1_non_diag*num_states/(num_states-1);
141 double derv2_diag = -derv1_diag*num_states/(num_states-1);
142
143 int nstates_sqr = num_states * num_states;
144 int i;
145 for (i = 0; i < nstates_sqr; i++)
146 if (i % (num_states+1) == 0) {
147 trans_matrix[i] = diag;
148 trans_derv1[i] = derv1_diag;
149 trans_derv2[i] = derv2_diag;
150 } else {
151 trans_matrix[i] = non_diag;
152 trans_derv1[i] = derv1_non_diag;
153 trans_derv2[i] = derv2_non_diag;
154 }
155
156 // DEBUG
157 /*int j;
158 if (verbose_mode == VB_DEBUG) {
159 cout.precision(4);
160 cout << "time = " << time << endl;
161 for (i = 0; i < num_states; i++, cout << endl) {
162 for (j = 0; j < num_states; j++) {
163 cout.width(8);
164 cout << right << trans_matrix[i*num_states+j] << " ";
165 }
166 cout << "| ";
167 for (j = 0; j < num_states; j++) {
168 cout << right << trans_derv1[i*num_states+j] << " ";
169 cout.width(8);
170 }
171 cout << "| ";
172 for (j = 0; j < num_states; j++) {
173 cout.width(8);
174 cout << right << trans_derv2[i*num_states+j] << " ";
175 }
176 }
177 cout.precision(10);
178 }*/
179
180 }
181
newTransMatrix()182 double *ModelSubst::newTransMatrix() {
183 return new double[num_states * num_states];
184 }
185
~ModelSubst()186 ModelSubst::~ModelSubst()
187 {
188 // mem space pointing to target model and thus avoid double free here
189 // if (linked_model && linked_model != this)
190 // return;
191
192 if (state_freq) delete [] state_freq;
193 }
194
195
196
197