← Back to branch summary

~ubuntu-branches/ubuntu/precise/libbpp-phyl/precise

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Model/K80.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/K80.cpp
 
1
//
 
2
// File: K80.cpp
 
3
// Created by: Julien Dutheil
 
4
// Created on: Tue May 27 15:24:30 2003
 
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 "K80.h"
 
41
 
 
42
// From the STL:
 
43
#include <cmath>
 
44
 
 
45
using namespace bpp;
 
46
 
 
47
#include <Bpp/Numeric/Matrix/MatrixTools.h>
 
48
 
 
49
using namespace std;
 
50
 
 
51
/******************************************************************************/
 
52
 
 
53
K80::K80(const NucleicAlphabet* alpha, double kappa) :
 
54
  AbstractReversibleSubstitutionModel(alpha, "K80."),
 
55
  kappa_(kappa), r_(), l_(), k_(), exp1_(), exp2_(), p_(size_, size_)
 
56
{
 
57
        Parameter kappaP(getNamespace() + "kappa", kappa, &Parameter::R_PLUS_STAR);
 
58
        addParameter_(kappaP);
 
59
        updateMatrices();
 
60
}
 
61
 
 
62
/******************************************************************************/
 
63
 
 
64
void K80::updateMatrices()
 
65
{
 
66
  kappa_ = getParameterValue("kappa");
 
67
        k_ = (kappa_ + 1.) / 2.;
 
68
        r_ = 4. / (kappa_ + 2.);
 
69
        
 
70
  // Frequences:
 
71
        freq_[0] = freq_[1] = freq_[2] = freq_[3] = 1. / 4.;
 
72
 
 
73
        // Generator:
 
74
        generator_(0, 0) = -2. - kappa_;
 
75
        generator_(1, 1) = -2. - kappa_;
 
76
        generator_(2, 2) = -2. - kappa_;
 
77
        generator_(3, 3) = -2. - kappa_;
 
78
 
 
79
        generator_(0, 1) = 1.;
 
80
        generator_(0, 3) = 1.;
 
81
        generator_(1, 0) = 1.;
 
82
        generator_(1, 2) = 1.;
 
83
        generator_(2, 1) = 1.;
 
84
        generator_(2, 3) = 1.;
 
85
        generator_(3, 0) = 1.;
 
86
        generator_(3, 2) = 1.;
 
87
        
 
88
        generator_(0, 2) = kappa_;
 
89
        generator_(1, 3) = kappa_;
 
90
        generator_(2, 0) = kappa_;
 
91
        generator_(3, 1) = kappa_;
 
92
 
 
93
        // Normalization:
 
94
        MatrixTools::scale(generator_, r_/4);
 
95
 
 
96
        // Exchangeability:
 
97
        exchangeability_ = generator_;
 
98
        MatrixTools::scale(exchangeability_, 4.);
 
99
 
 
100
        // Eigen values:
 
101
        eigenValues_[0] = 0;
 
102
        eigenValues_[1] = -r_ * (1. + kappa_)/2;
 
103
        eigenValues_[2] = -r_ * (1. + kappa_)/2;
 
104
        eigenValues_[3] = -r_;
 
105
        
 
106
        // Eigen vectors:
 
107
        leftEigenVectors_(0,0) = 1. / 4.;
 
108
        leftEigenVectors_(0,1) = 1. / 4.;
 
109
        leftEigenVectors_(0,2) = 1. / 4.;
 
110
        leftEigenVectors_(0,3) = 1. / 4.;
 
111
        leftEigenVectors_(1,0) = 0.;
 
112
        leftEigenVectors_(1,1) = 1. / 2.;
 
113
        leftEigenVectors_(1,2) = 0.;
 
114
        leftEigenVectors_(1,3) = -1. / 2.;
 
115
        leftEigenVectors_(2,0) = 1. / 2.;
 
116
        leftEigenVectors_(2,1) = 0.;
 
117
        leftEigenVectors_(2,2) = -1. / 2.;
 
118
        leftEigenVectors_(2,3) = 0.;
 
119
        leftEigenVectors_(3,0) = 1. / 4.;
 
120
        leftEigenVectors_(3,1) = -1. / 4.;
 
121
        leftEigenVectors_(3,2) = 1. / 4.;
 
122
        leftEigenVectors_(3,3) = -1. / 4.;
 
123
 
 
124
        rightEigenVectors_(0,0) = 1.;
 
125
        rightEigenVectors_(0,1) = 0.;
 
126
        rightEigenVectors_(0,2) = 1.;
 
127
        rightEigenVectors_(0,3) = 1.;
 
128
        rightEigenVectors_(1,0) = 1.;
 
129
        rightEigenVectors_(1,1) = 1.;
 
130
        rightEigenVectors_(1,2) = 0.;
 
131
        rightEigenVectors_(1,3) = -1.;
 
132
        rightEigenVectors_(2,0) = 1.;
 
133
        rightEigenVectors_(2,1) = 0.;
 
134
        rightEigenVectors_(2,2) = -1.;
 
135
        rightEigenVectors_(2,3) = 1.;
 
136
        rightEigenVectors_(3,0) = 1.;
 
137
        rightEigenVectors_(3,1) = -1.;
 
138
        rightEigenVectors_(3,2) = 0;
 
139
        rightEigenVectors_(3,3) = -1.;
 
140
}
 
141
        
 
142
/******************************************************************************/
 
143
 
 
144
double K80::Pij_t(int i, int j, double d) const
 
145
{
 
146
        l_ = rate_ * r_ * d;
 
147
        exp1_ = exp(-l_);
 
148
        exp2_ = exp(-k_ * l_);
 
149
        
 
150
        switch(i) {
 
151
                //A
 
152
                case 0 : {
 
153
                        switch(j) {
 
154
                                case 0 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //A
 
155
                                case 1 : return 0.25 * (1. - exp1_);               //C
 
156
                                case 2 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //G
 
157
                                case 3 : return 0.25 * (1. - exp1_);               //T, U
 
158
                        }
 
159
                } 
 
160
                //C
 
161
                case 1 : {
 
162
                        switch(j) {
 
163
                                case 0 : return 0.25 * (1. - exp1_);               //A
 
164
                                case 1 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //C
 
165
                                case 2 : return 0.25 * (1. - exp1_);               //G
 
166
                                case 3 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //T, U
 
167
                        }
 
168
                }
 
169
                //G
 
170
                case 2 : {
 
171
                        switch(j) {
 
172
                                case 0 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //A
 
173
                                case 1 : return 0.25 * (1. - exp1_);               //C
 
174
                                case 2 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //G
 
175
                                case 3 : return 0.25 * (1. - exp1_);               //T, U
 
176
                        }
 
177
                }
 
178
                //T, U
 
179
                case 3 : {
 
180
                        switch(j) {
 
181
                                case 0 : return 0.25 * (1. - exp1_);               //A
 
182
                                case 1 : return 0.25 * (1. + exp1_) - 0.5 * exp2_; //C
 
183
                                case 2 : return 0.25 * (1. - exp1_);               //G
 
184
                                case 3 : return 0.25 * (1. + exp1_) + 0.5 * exp2_; //T, U
 
185
                        }
 
186
                }
 
187
        }
 
188
        return 0;
 
189
}
 
190
 
 
191
/******************************************************************************/
 
192
 
 
193
double K80::dPij_dt(int i, int j, double d) const
 
194
{
 
195
  l_ = rate_ * r_ * d;
 
196
        exp1_ = exp(-l_);
 
197
        exp2_ = exp(-k_ * l_);
 
198
 
 
199
        switch(i) {
 
200
                //A
 
201
                case 0 : {
 
202
                        switch(j) {
 
203
                        case 0 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //A
 
204
                        case 1 : return rate_ * r_/4. * (  exp1_);                   //C
 
205
                        case 2 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //G
 
206
                        case 3 : return rate_ * r_/4. * (  exp1_);                   //T, U
 
207
                        }
 
208
                } 
 
209
                //C
 
210
                case 1 : {
 
211
                        switch(j) {
 
212
                        case 0 : return rate_ * r_/4. * (  exp1_);                   //A
 
213
                        case 1 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //C
 
214
                        case 2 : return rate_ * r_/4. * (  exp1_);                   //G
 
215
                        case 3 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //T, U
 
216
                        }
 
217
                }
 
218
                //G
 
219
                case 2 : {
 
220
                        switch(j) {
 
221
                        case 0 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //A
 
222
                        case 1 : return rate_ * r_/4. * (  exp1_);                   //C
 
223
                        case 2 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //G
 
224
                        case 3 : return rate_ * r_/4. * (  exp1_);                   //T, U
 
225
                        }
 
226
                }
 
227
                //T, U
 
228
                case 3 : {
 
229
                        switch(j) {
 
230
                        case 0 : return rate_ * r_/4. * (  exp1_);                   //A
 
231
                        case 1 : return rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //C
 
232
                        case 2 : return rate_ * r_/4. * (  exp1_);                   //G
 
233
                        case 3 : return rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //T, U
 
234
                        }
 
235
                }
 
236
        }
 
237
        return 0;
 
238
}
 
239
 
 
240
/******************************************************************************/
 
241
 
 
242
double K80::d2Pij_dt2(int i, int j, double d) const
 
243
{
 
244
        double k_2 = k_ * k_;
 
245
        double r_2 = rate_ * rate_ * r_ * r_;
 
246
        l_ = rate_ * r_ * d;
 
247
        exp1_ = exp(-l_);
 
248
        exp2_ = exp(-k_ * l_);
 
249
 
 
250
        switch(i) {
 
251
                //A
 
252
                case 0 : {
 
253
                        switch(j) {
 
254
                                case 0 : return r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //A
 
255
                                case 1 : return r_2/4. * (- exp1_);                    //C
 
256
                                case 2 : return r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //G
 
257
                                case 3 : return r_2/4. * (- exp1_);                    //T, U
 
258
                        }
 
259
                } 
 
260
                //C
 
261
                case 1 : {
 
262
                        switch(j) {
 
263
                                case 0 : return r_2/4. * (- exp1_);                    //A
 
264
                                case 1 : return r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //C
 
265
                                case 2 : return r_2/4. * (- exp1_);                    //G
 
266
                                case 3 : return r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //T, U
 
267
                        }
 
268
                }
 
269
                //G
 
270
                case 2 : {
 
271
                        switch(j) {
 
272
                                case 0 : return r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //A
 
273
                                case 1 : return r_2/4. * (- exp1_);                    //C
 
274
                                case 2 : return r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //G
 
275
                                case 3 : return r_2/4. * (- exp1_);                    //T, U
 
276
                        }
 
277
                }
 
278
                //T, U
 
279
                case 3 : {
 
280
                        switch(j) {
 
281
                                case 0 : return r_2/4. * (- exp1_);                    //A
 
282
                                case 1 : return r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //C
 
283
                                case 2 : return r_2/4. * (- exp1_);                    //G
 
284
                                case 3 : return r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //T, U
 
285
                        }
 
286
                }
 
287
        }
 
288
        return 0;
 
289
}
 
290
 
 
291
/******************************************************************************/
 
292
 
 
293
const Matrix<double> & K80::getPij_t(double d) const
 
294
{
 
295
  l_ = rate_ * r_ * d;
 
296
        exp1_ = exp(-l_);
 
297
        exp2_ = exp(-k_ * l_);
 
298
 
 
299
        //A
 
300
        p_(0, 0) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //A
 
301
        p_(0, 1) = 0.25 * (1. - exp1_);               //C
 
302
        p_(0, 2) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //G
 
303
        p_(0, 3) = 0.25 * (1. - exp1_);               //T, U
 
304
 
 
305
        //C
 
306
        p_(1, 0) = 0.25 * (1. - exp1_);               //A
 
307
        p_(1, 1) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //C
 
308
        p_(1, 2) = 0.25 * (1. - exp1_);               //G
 
309
        p_(1, 3) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //T, U
 
310
 
 
311
        //G
 
312
        p_(2, 0) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //A
 
313
        p_(2, 1) = 0.25 * (1. - exp1_);               //C
 
314
        p_(2, 2) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //G
 
315
        p_(2, 3) = 0.25 * (1. - exp1_);               //T, U
 
316
 
 
317
        //T, U
 
318
        p_(3, 0) = 0.25 * (1. - exp1_);               //A
 
319
        p_(3, 1) = 0.25 * (1. + exp1_) - 0.5 * exp2_; //C
 
320
        p_(3, 2) = 0.25 * (1. - exp1_);               //G
 
321
        p_(3, 3) = 0.25 * (1. + exp1_) + 0.5 * exp2_; //T, U
 
322
 
 
323
        return p_;
 
324
}
 
325
 
 
326
const Matrix<double> & K80::getdPij_dt(double d) const
 
327
{
 
328
  l_ = rate_ * r_ * d;
 
329
        exp1_ = exp(-l_);
 
330
        exp2_ = exp(-k_ * l_);
 
331
 
 
332
        p_(0, 0) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //A
 
333
        p_(0, 1) = rate_ * r_/4. * (  exp1_);                   //C
 
334
        p_(0, 2) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //G
 
335
        p_(0, 3) = rate_ * r_/4. * (  exp1_);                   //T, U
 
336
 
 
337
        //C
 
338
        p_(1, 0) = rate_ * r_/4. * (  exp1_);                   //A
 
339
        p_(1, 1) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //C
 
340
        p_(1, 2) = rate_ * r_/4. * (  exp1_);                   //G
 
341
        p_(1, 3) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //T, U
 
342
 
 
343
        //G
 
344
        p_(2, 0) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //A
 
345
        p_(2, 1) = rate_ * r_/4. * (  exp1_);                   //C
 
346
        p_(2, 2) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //G
 
347
        p_(2, 3) = rate_ * r_/4. * (  exp1_);                   //T, U
 
348
 
 
349
        //T, U
 
350
        p_(3, 0) = rate_ * r_/4. * (  exp1_);                   //A
 
351
        p_(3, 1) = rate_ * r_/4. * (- exp1_ + 2. * k_ * exp2_); //C
 
352
        p_(3, 2) = rate_ * r_/4. * (  exp1_);                   //G
 
353
        p_(3, 3) = rate_ * r_/4. * (- exp1_ - 2. * k_ * exp2_); //T, U
 
354
 
 
355
        return p_;
 
356
}
 
357
 
 
358
const Matrix<double> & K80::getd2Pij_dt2(double d) const
 
359
{
 
360
        double k_2 = k_ * k_;
 
361
        double r_2 = rate_ * rate_ * r_ * r_;
 
362
        l_ = rate_ * r_ * d;
 
363
        exp1_ = exp(-l_);
 
364
        exp2_ = exp(-k_ * l_);
 
365
 
 
366
        p_(0, 0) = r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //A
 
367
        p_(0, 1) = r_2/4. * (- exp1_);                    //C
 
368
        p_(0, 2) = r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //G
 
369
        p_(0, 3) = r_2/4. * (- exp1_);                    //T, U
 
370
 
 
371
        //C
 
372
        p_(1, 0) = r_2/4. * (- exp1_);                    //A
 
373
        p_(1, 1) = r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //C
 
374
        p_(1, 2) = r_2/4. * (- exp1_);                    //G
 
375
        p_(1, 3) = r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //T, U
 
376
 
 
377
        //G
 
378
        p_(2, 0) = r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //A
 
379
        p_(2, 1) = r_2/4. * (- exp1_);                    //C
 
380
        p_(2, 2) = r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //G
 
381
        p_(2, 3) = r_2/4. * (- exp1_);                    //T, U
 
382
 
 
383
        //T, U
 
384
        p_(3, 0) = r_2/4. * (- exp1_);                    //A
 
385
        p_(3, 1) = r_2/4. * (  exp1_ - 2. * k_2 * exp2_); //C
 
386
        p_(3, 2) = r_2/4. * (- exp1_);                    //G
 
387
        p_(3, 3) = r_2/4. * (  exp1_ + 2. * k_2 * exp2_); //T, U
 
388
 
 
389
        return p_;
 
390
}
 
391
 
 
392
/******************************************************************************/
 
393