← Back to branch summary

~ubuntu-branches/ubuntu/raring/libbpp-phyl/raring

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Model/HKY85.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Julien Dutheil
  • Date: 2011-06-09 11:00:00 UTC
  • Revision ID: james.westby@ubuntu.com-20110609110000-yvx78svv6w7xxgph
Tags: upstream-2.0.2
Import upstream version 2.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Bpp/Phyl/Model/HKY85.cpp
 
1
//
 
2
// File: HKY85.cpp
 
3
// Created by: Julien Dutheil
 
4
// Created on: Thu Jan 22 16:17:39 2004
 
5
//
 
6
 
 
7
/*
 
8
Copyright or © or Copr. CNRS, (November 16, 2004)
 
9
 
 
10
This software is a computer program whose purpose is to provide classes
 
11
for phylogenetic data analysis.
 
12
 
 
13
This software is governed by the CeCILL  license under French law and
 
14
abiding by the rules of distribution of free software.  You can  use, 
 
15
modify and/ or redistribute the software under the terms of the CeCILL
 
16
license as circulated by CEA, CNRS and INRIA at the following URL
 
17
"http://www.cecill.info". 
 
18
 
 
19
As a counterpart to the access to the source code and  rights to copy,
 
20
modify and redistribute granted by the license, users are provided only
 
21
with a limited warranty  and the software's author,  the holder of the
 
22
economic rights,  and the successive licensors  have only  limited
 
23
liability. 
 
24
 
 
25
In this respect, the user's attention is drawn to the risks associated
 
26
with loading,  using,  modifying and/or developing or reproducing the
 
27
software by the user in light of its specific status of free software,
 
28
that may mean  that it is complicated to manipulate,  and  that  also
 
29
therefore means  that it is reserved for developers  and  experienced
 
30
professionals having in-depth computer knowledge. Users are therefore
 
31
encouraged to load and test the software's suitability as regards their
 
32
requirements in conditions enabling the security of their systems and/or 
 
33
data to be ensured and,  more generally, to use and operate it in the 
 
34
same conditions as regards security. 
 
35
 
 
36
The fact that you are presently reading this means that you have had
 
37
knowledge of the CeCILL license and that you accept its terms.
 
38
*/
 
39
 
 
40
#include "HKY85.h"
 
41
 
 
42
#include <Bpp/Numeric/Matrix/MatrixTools.h>
 
43
 
 
44
// From SeqLib:
 
45
#include <Bpp/Seq/Container/SequenceContainerTools.h>
 
46
 
 
47
using namespace bpp;
 
48
 
 
49
// From the STL:
 
50
#include <cmath>
 
51
 
 
52
using namespace std;
 
53
 
 
54
/******************************************************************************/
 
55
 
 
56
HKY85::HKY85(
 
57
        const NucleicAlphabet* alpha,
 
58
        double kappa,
 
59
        double piA,
 
60
        double piC,
 
61
        double piG,
 
62
        double piT):
 
63
        AbstractReversibleSubstitutionModel(alpha, "HKY85."),
 
64
  kappa_(kappa), k1_(), k2_(), r_(),
 
65
  piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
 
66
  theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
 
67
  exp1_(), exp21_(), exp22_(), l_(), p_(size_, size_)
 
68
{
 
69
        Parameter kappaP("HKY85.kappa", kappa, &Parameter::R_PLUS_STAR);
 
70
        addParameter_(kappaP);
 
71
        Parameter thetaP("HKY85.theta" , theta_, new IncludingInterval(0.001, 0.999), true);
 
72
        addParameter_(thetaP);
 
73
        Parameter theta1P("HKY85.theta1", theta1_, new IncludingInterval(0.001, 0.999), true);
 
74
        addParameter_(theta1P);
 
75
        Parameter theta2P("HKY85.theta2", theta2_, new IncludingInterval(0.001, 0.999), true);
 
76
        addParameter_(theta2P);
 
77
        updateMatrices();
 
78
}
 
79
 
 
80
/******************************************************************************/
 
81
 
 
82
void HKY85::updateMatrices()
 
83
{
 
84
        kappa_  = getParameterValue("kappa");
 
85
        theta_  = getParameterValue("theta");
 
86
        theta1_ = getParameterValue("theta1");
 
87
        theta2_ = getParameterValue("theta2");
 
88
  piA_ = theta1_ * (1. - theta_);
 
89
  piC_ = (1. - theta2_) * theta_;
 
90
  piG_ = theta2_ * theta_;
 
91
  piT_ = (1. - theta1_) * (1. - theta_);
 
92
        piR_   = piA_ + piG_;
 
93
        piY_   = piT_ + piC_;
 
94
        k1_    = kappa_ * piY_ + piR_;
 
95
        k2_    = kappa_ * piR_ + piY_;
 
96
 
 
97
  freq_[0] = piA_;
 
98
  freq_[1] = piC_;
 
99
  freq_[2] = piG_;
 
100
  freq_[3] = piT_;
 
101
        
 
102
        generator_(0, 0) = -(                     piC_ + kappa_*piG_ +        piT_);
 
103
        generator_(1, 1) = -(       piA_ +                      piG_ + kappa_*piT_); 
 
104
        generator_(2, 2) = -(kappa_*piA_ +        piC_               +       piT_);
 
105
        generator_(3, 3) = -(       piA_ + kappa_*piC_ +        piG_             );
 
106
 
 
107
        generator_(1, 0) = piA_;
 
108
        generator_(3, 0) = piA_;
 
109
        generator_(0, 1) = piC_;
 
110
        generator_(2, 1) = piC_;
 
111
        generator_(1, 2) = piG_;
 
112
        generator_(3, 2) = piG_;
 
113
        generator_(0, 3) = piT_;
 
114
        generator_(2, 3) = piT_;
 
115
        
 
116
        generator_(2, 0) = kappa_ * piA_;
 
117
        generator_(3, 1) = kappa_ * piC_;
 
118
        generator_(0, 2) = kappa_ * piG_;
 
119
        generator_(1, 3) = kappa_ * piT_;
 
120
        
 
121
        // Normalization:
 
122
        r_ = 1. / (2. * (piA_ * piC_ + piC_ * piG_ + piA_ * piT_ + piG_ * piT_ + kappa_ * (piC_ * piT_ + piA_ * piG_)));
 
123
        MatrixTools::scale(generator_, r_);
 
124
        
 
125
        // Exchangeability:
 
126
        exchangeability_(0,0) = generator_(0,0) / piA_;
 
127
        exchangeability_(0,1) = generator_(0,1) / piC_; 
 
128
        exchangeability_(0,2) = generator_(0,2) / piG_; 
 
129
        exchangeability_(0,3) = generator_(0,3) / piT_;
 
130
 
 
131
        exchangeability_(1,0) = generator_(1,0) / piA_; 
 
132
        exchangeability_(1,1) = generator_(1,1) / piC_; 
 
133
        exchangeability_(1,2) = generator_(1,2) / piG_; 
 
134
        exchangeability_(1,3) = generator_(1,3) / piT_; 
 
135
        
 
136
        exchangeability_(2,0) = generator_(2,0) / piA_; 
 
137
        exchangeability_(2,1) = generator_(2,1) / piC_; 
 
138
        exchangeability_(2,2) = generator_(2,2) / piG_; 
 
139
        exchangeability_(2,3) = generator_(2,3) / piT_; 
 
140
        
 
141
        exchangeability_(3,0) = generator_(3,0) / piA_;
 
142
        exchangeability_(3,1) = generator_(3,1) / piC_; 
 
143
        exchangeability_(3,2) = generator_(3,2) / piG_; 
 
144
        exchangeability_(3,3) = generator_(3,3) / piT_;
 
145
 
 
146
        // Eigen values:
 
147
        eigenValues_[0] = 0;
 
148
        eigenValues_[1] = -r_ * (kappa_ * piY_ + piR_);
 
149
        eigenValues_[2] = -r_ * (kappa_ * piR_ + piY_); 
 
150
        eigenValues_[3] = -r_;
 
151
        
 
152
        // Eigen vectors:
 
153
        leftEigenVectors_(0,0) = piA_;
 
154
        leftEigenVectors_(0,1) = piC_;
 
155
        leftEigenVectors_(0,2) = piG_;
 
156
        leftEigenVectors_(0,3) = piT_;
 
157
 
 
158
        leftEigenVectors_(1,0) = 0.;
 
159
        leftEigenVectors_(1,1) = piT_ / piY_;
 
160
        leftEigenVectors_(1,2) = 0.;
 
161
        leftEigenVectors_(1,3) = -piT_ / piY_;
 
162
 
 
163
        leftEigenVectors_(2,0) = piG_ / piR_;
 
164
        leftEigenVectors_(2,1) = 0.;
 
165
        leftEigenVectors_(2,2) = -piG_ / piR_;
 
166
        leftEigenVectors_(2,3) = 0.;
 
167
 
 
168
        leftEigenVectors_(3,0) = piA_*piY_ / piR_;
 
169
        leftEigenVectors_(3,1) = -piC_;
 
170
        leftEigenVectors_(3,2) = piG_*piY_ / piR_;
 
171
        leftEigenVectors_(3,3) = -piT_;
 
172
 
 
173
        rightEigenVectors_(0,0) = 1.;
 
174
        rightEigenVectors_(0,1) = 0.;
 
175
        rightEigenVectors_(0,2) = 1.;
 
176
        rightEigenVectors_(0,3) = 1.;
 
177
        
 
178
        rightEigenVectors_(1,0) = 1.;
 
179
        rightEigenVectors_(1,1) = 1.;
 
180
        rightEigenVectors_(1,2) = 0.;;
 
181
        rightEigenVectors_(1,3) = -piR_ / piY_;
 
182
 
 
183
        rightEigenVectors_(2,0) = 1.;
 
184
        rightEigenVectors_(2,1) = 0.;
 
185
        rightEigenVectors_(2,2) = -piA_ / piG_;
 
186
        rightEigenVectors_(2,3) = 1.;
 
187
 
 
188
        rightEigenVectors_(3,0) = 1.;
 
189
        rightEigenVectors_(3,1) = -piC_ / piT_;
 
190
        rightEigenVectors_(3,2) = 0.;
 
191
        rightEigenVectors_(3,3) = -piR_ / piY_;
 
192
}
 
193
        
 
194
/******************************************************************************/
 
195
 
 
196
double HKY85::Pij_t(int i, int j, double d) const
 
197
{
 
198
        l_     = rate_ * r_ * d;
 
199
        exp1_  = exp(-l_);
 
200
        exp22_ = exp(-k2_ * l_);
 
201
        exp21_ = exp(-k1_ * l_);
 
202
        
 
203
        switch(i)
 
204
  {
 
205
                //A
 
206
                case 0 : {
 
207
                        switch(j) {
 
208
                                case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp22_; //A
 
209
                                case 1 : return piC_ * (1. -               exp1_);                        //C
 
210
                                case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp22_; //G
 
211
                                case 3 : return piT_ * (1. -               exp1_);                        //T, U
 
212
                        }
 
213
                } 
 
214
                //C
 
215
                case 1 : {
 
216
                        switch(j) {
 
217
                                case 0 : return piA_ * (1. -               exp1_);                        //A
 
218
                                case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp21_; //C
 
219
                                case 2 : return piG_ * (1. -               exp1_);                        //G
 
220
                                case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp21_; //T, U
 
221
                        }
 
222
                }
 
223
                //G
 
224
                case 2 : {
 
225
                        switch(j) {
 
226
                                case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp22_; //A
 
227
                                case 1 : return piC_ * (1. -               exp1_);                        //C
 
228
                                case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp22_; //G
 
229
                                case 3 : return piT_ * (1. -               exp1_);                        //T, U
 
230
                        }
 
231
                }
 
232
                //T, U
 
233
                case 3 : {
 
234
                        switch(j) {
 
235
                                case 0 : return piA_ * (1. -               exp1_);                        //A
 
236
                                case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp21_; //C
 
237
                                case 2 : return piG_ * (1. -               exp1_);                        //G
 
238
                                case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp21_; //T, U
 
239
                        }
 
240
                }
 
241
        }
 
242
        return 0;
 
243
}
 
244
 
 
245
/******************************************************************************/
 
246
 
 
247
double HKY85::dPij_dt(int i, int j, double d) const
 
248
{
 
249
        l_     = rate_ * r_ * d;
 
250
        exp1_  = exp(-l_);
 
251
        exp22_ = exp(-k2_ * l_);
 
252
        exp21_ = exp(-k1_ * l_);
 
253
        
 
254
        switch(i)
 
255
  {
 
256
                //A
 
257
                case 0 : {
 
258
                        switch(j) {
 
259
                                case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp22_); //A
 
260
                                case 1 : return rate_ * r_ * (piC_ *                exp1_);                              //C
 
261
                                case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp22_); //G
 
262
                                case 3 : return rate_ * r_ * (piT_ *                exp1_);                              //T, U
 
263
                        }
 
264
                } 
 
265
                //C
 
266
                case 1 : {
 
267
                        switch(j) {
 
268
                                case 0 : return rate_ * r_ * (piA_ *                exp1_);                              //A
 
269
                                case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k1_ * exp21_); //C
 
270
                                case 2 : return rate_ * r_ * (piG_ *                exp1_);                              //G
 
271
                                case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k1_ * exp21_); //T, U
 
272
                        }
 
273
                }
 
274
                //G
 
275
                case 2 : {
 
276
                        switch(j) {
 
277
                                case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp22_); //A
 
278
                                case 1 : return rate_ * r_ * (piC_ *                exp1_);                              //C
 
279
                                case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp22_); //G
 
280
                                case 3 : return rate_ * r_ * (piT_ *                exp1_);                              //T, U
 
281
                        }
 
282
                }
 
283
                //T, U
 
284
                case 3 : {
 
285
                        switch(j) {
 
286
                                case 0 : return rate_ * r_ * (piA_ *                exp1_);                              //A
 
287
                                case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k1_ * exp21_); //C
 
288
                                case 2 : return rate_ * r_ * (piG_ *                exp1_);                              //G
 
289
                                case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k1_ * exp21_); //T, U
 
290
                        }
 
291
                }
 
292
        }
 
293
        return 0;
 
294
}
 
295
 
 
296
/******************************************************************************/
 
297
 
 
298
double HKY85::d2Pij_dt2(int i, int j, double d) const
 
299
{
 
300
        double r_2 = rate_ * rate_ * r_ * r_;
 
301
        l_ = rate_ * r_ * d;
 
302
        double k1_2 = k1_ * k1_;
 
303
        double k2_2 = k2_ * k2_;
 
304
        exp1_ = exp(-l_);
 
305
        exp22_ = exp(-k2_ * l_);
 
306
        exp21_ = exp(-k1_ * l_);
 
307
        
 
308
        switch(i)
 
309
  {
 
310
                //A
 
311
                case 0 : {
 
312
                        switch(j) {
 
313
                                case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp22_); //A
 
314
                                case 1 : return r_2 * (piC_ *             - exp1_);                               //C
 
315
                                case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp22_); //G
 
316
                                case 3 : return r_2 * (piT_ *             - exp1_);                               //T, U
 
317
                        }
 
318
                } 
 
319
                //C
 
320
                case 1 : {
 
321
                        switch(j) {
 
322
                                case 0 : return r_2 * (piA_ *             - exp1_);                               //A
 
323
                                case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k1_2 * exp21_); //C
 
324
                                case 2 : return r_2 * (piG_ *             - exp1_);                               //G
 
325
                                case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k1_2 * exp21_); //T, U
 
326
                        }
 
327
                }
 
328
                //G
 
329
                case 2 : {
 
330
                        switch(j) {
 
331
                                case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp22_); //A
 
332
                                case 1 : return r_2 * (piC_ *             - exp1_);                               //C
 
333
                                case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp22_); //G
 
334
                                case 3 : return r_2 * (piT_ *             - exp1_);                               //T, U
 
335
                        }
 
336
                }
 
337
                //T, U
 
338
                case 3 : {
 
339
                        switch(j) {
 
340
                                case 0 : return r_2 * (piA_ *             - exp1_);                              //A
 
341
                                case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k1_2 * exp21_); //C
 
342
                                case 2 : return r_2 * (piG_ *             - exp1_);                              //G
 
343
                                case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k1_2 * exp21_); //T, U
 
344
                        }
 
345
                }
 
346
        }
 
347
        return 0;
 
348
}
 
349
 
 
350
/******************************************************************************/
 
351
 
 
352
const Matrix<double> & HKY85::getPij_t(double d) const
 
353
{
 
354
  l_ = rate_ * r_ * d;
 
355
        exp1_ = exp(-l_);
 
356
        exp22_ = exp(-k2_ * l_);
 
357
        exp21_ = exp(-k1_ * l_);
 
358
 
 
359
        //A
 
360
        p_(0, 0) = piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp22_; //A
 
361
        p_(0, 1) = piC_ * (1. -               exp1_);                        //C
 
362
        p_(0, 2) = piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp22_; //G
 
363
        p_(0, 3) = piT_ * (1. -               exp1_);                        //T, U
 
364
 
 
365
        //C
 
366
        p_(1, 0) = piA_ * (1. -               exp1_);                        //A
 
367
        p_(1, 1) = piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp21_; //C
 
368
        p_(1, 2) = piG_ * (1. -               exp1_);                        //G
 
369
        p_(1, 3) = piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp21_; //T, U
 
370
 
 
371
        //G
 
372
        p_(2, 0) = piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp22_; //A
 
373
        p_(2, 1) = piC_ * (1. -               exp1_);                        //C
 
374
        p_(2, 2) = piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp22_; //G
 
375
        p_(2, 3) = piT_ * (1. -               exp1_);                        //T, U
 
376
 
 
377
        //T, U
 
378
        p_(3, 0) = piA_ * (1. -               exp1_);                        //A
 
379
        p_(3, 1) = piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp21_; //C
 
380
        p_(3, 2) = piG_ * (1. -               exp1_);                        //G
 
381
        p_(3, 3) = piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp21_; //T, U
 
382
 
 
383
        return p_;
 
384
}
 
385
 
 
386
const Matrix<double> & HKY85::getdPij_dt(double d) const
 
387
{
 
388
        l_ = rate_ * r_ * d;
 
389
        exp1_ = exp(-l_);
 
390
        exp22_ = exp(-k2_ * l_);
 
391
        exp21_ = exp(-k1_ * l_);
 
392
 
 
393
        //A
 
394
        p_(0, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp22_); //A
 
395
        p_(0, 1) = rate_ * r_ * (piC_ *                exp1_);                              //C
 
396
        p_(0, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp22_); //G
 
397
        p_(0, 3) = rate_ * r_ * (piT_ *                exp1_);                              //T, U
 
398
 
 
399
        //C
 
400
        p_(1, 0) = rate_ * r_ * (piA_ *                exp1_);                              //A
 
401
        p_(1, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k1_ * exp21_); //C
 
402
        p_(1, 2) = rate_ * r_ * (piG_ *                exp1_);                              //G
 
403
        p_(1, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k1_ * exp21_); //T, U
 
404
 
 
405
        //G
 
406
        p_(2, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp22_); //A
 
407
        p_(2, 1) = rate_ * r_ * (piC_ *                exp1_);                              //C
 
408
        p_(2, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp22_); //G
 
409
        p_(2, 3) = rate_ * r_ * (piT_ *                exp1_);                              //T, U
 
410
 
 
411
        //T, U
 
412
        p_(3, 0) = rate_ * r_ * (piA_ *                exp1_);                              //A
 
413
        p_(3, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k1_ * exp21_); //C
 
414
        p_(3, 2) = rate_ * r_ * (piG_ *                exp1_);                              //G
 
415
        p_(3, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k1_ * exp21_); //T, U
 
416
 
 
417
        return p_;
 
418
}
 
419
 
 
420
const Matrix<double> & HKY85::getd2Pij_dt2(double d) const
 
421
{
 
422
        double r_2 = rate_ * rate_ * r_ * r_;
 
423
        l_ = rate_ * r_ * d;
 
424
        double k1_2 = k1_ * k1_;
 
425
        double k2_2 = k2_ * k2_;
 
426
        exp1_ = exp(-l_);
 
427
        exp22_ = exp(-k2_ * l_);
 
428
        exp21_ = exp(-k1_ * l_);
 
429
 
 
430
        //A
 
431
        p_(0, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp22_); //A
 
432
        p_(0, 1) = r_2 * (piC_ *             - exp1_);                               //C
 
433
        p_(0, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp22_); //G
 
434
        p_(0, 3) = r_2 * (piT_ *             - exp1_);                               //T, U
 
435
 
 
436
        //C
 
437
        p_(1, 0) = r_2 * (piA_ *             - exp1_);                               //A
 
438
        p_(1, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k1_2 * exp21_); //C
 
439
        p_(1, 2) = r_2 * (piG_ *             - exp1_);                               //G
 
440
        p_(1, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k1_2 * exp21_); //T, U
 
441
 
 
442
        //G
 
443
        p_(2, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp22_); //A
 
444
        p_(2, 1) = r_2 * (piC_ *             - exp1_);                               //C
 
445
        p_(2, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp22_); //G
 
446
        p_(2, 3) = r_2 * (piT_ *             - exp1_);                               //T, U
 
447
 
 
448
        //T, U
 
449
        p_(3, 0) = r_2 * (piA_ *             - exp1_);                               //A
 
450
        p_(3, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k1_2 * exp21_); //C
 
451
        p_(3, 2) = r_2 * (piG_ *             - exp1_);                               //G
 
452
        p_(3, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k1_2 * exp21_); //T, U
 
453
 
 
454
        return p_;
 
455
}
 
456
 
 
457
/******************************************************************************/
 
458
 
 
459
void HKY85::setFreq(std::map<int, double>& freqs)
 
460
{
 
461
  piA_ = freqs[0];
 
462
  piC_ = freqs[1];
 
463
  piG_ = freqs[2];
 
464
  piT_ = freqs[3];
 
465
  vector<string> thetas(3);
 
466
  thetas[0] = getNamespace() + "theta";
 
467
  thetas[1] = getNamespace() + "theta1";
 
468
  thetas[2] = getNamespace() + "theta2";
 
469
  ParameterList pl = getParameters().subList(thetas);
 
470
  pl[0].setValue(piC_ + piG_);
 
471
  pl[1].setValue(piA_ / (piA_ + piT_));
 
472
  pl[2].setValue(piG_ / (piC_ + piG_));
 
473
  setParametersValues(pl);
 
474
}
 
475
 
 
476
/******************************************************************************/
 
477