← Back to branch summary

~ubuntu-branches/ubuntu/saucy/libbpp-phyl/saucy

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Julien Dutheil
  • Date: 2013-02-09 16:00:00 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20130209160000-5v65ba68z8032nzj
Tags: 2.0.3-1
* Reorganized model hierarchy
* New pairwise models
* Several bugs fixed

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Bpp/Phyl/Model/YpR.cpp
 
1
//
 
2
// File: YpR.cpp
 
3
// Created by: Laurent Gueguen
 
4
// Created on: Thu August 2 2007
 
5
//
 
6
 
 
7
/*
 
8
   Copyright or � or Copr. CNRS, (November 16, 2004)
 
9
   This software is a computer program whose purpose is to provide
 
10
   classes for phylogenetic data analysis.
 
11
 
 
12
   This software is governed by the CeCILL license under French law and
 
13
   abiding by the rules of distribution of free software. You can use,
 
14
   modify and/ or redistribute the software under the terms of the CeCILL
 
15
   license as circulated by CEA, CNRS and INRIA at the following URL
 
16
   "http://www.cecill.info".
 
17
 
 
18
   As a counterpart to the access to the source code and rights to copy,
 
19
   modify and redistribute granted by the license, users are provided
 
20
   only with a limited warranty and the software's author, the holder of
 
21
   the economic rights, and the successive licensors have only limited
 
22
   liability.
 
23
 
 
24
   In this respect, the user's attention is drawn to the risks associated
 
25
   with loading, using, modifying and/or developing or reproducing the
 
26
   software by the user in light of its specific status of free software,
 
27
   that may mean that it is complicated to manipulate, and that also
 
28
   therefore means that it is reserved for developers and experienced
 
29
   professionals having in-depth computer knowledge. Users are therefore
 
30
   encouraged to load and test the software's suitability as regards
 
31
   their requirements in conditions enabling the security of their
 
32
   systems and/or data to be ensured and, more generally, to use and
 
33
   operate it in the same conditions as regards security.
 
34
 
 
35
   The fact that you are presently reading this means that you have had
 
36
   knowledge of the CeCILL license and that you accept its terms.
 
37
 */
 
38
 
 
39
#include "YpR.h"
 
40
 
 
41
// From the STL:
 
42
#include <cmath>
 
43
 
 
44
using namespace bpp;
 
45
 
 
46
#include <Bpp/Numeric/Matrix/MatrixTools.h>
 
47
#include <Bpp/Numeric/VectorTools.h>
 
48
#include <Bpp/Numeric/Matrix/EigenValue.h>
 
49
 
 
50
/******************************************************************************/
 
51
 
 
52
YpR::YpR(const RNY* alph, SubstitutionModel* const pm, const std::string& prefix) :
 
53
  AbstractParameterAliasable(prefix),
 
54
  AbstractSubstitutionModel(alph,prefix),
 
55
  _pmodel(pm->clone()),
 
56
  _nestedPrefix(pm->getNamespace())
 
57
{
 
58
  _pmodel->setNamespace(prefix + _nestedPrefix);
 
59
  _pmodel->enableEigenDecomposition(0);
 
60
  addParameters_(_pmodel->getParameters());
 
61
}
 
62
 
 
63
YpR::YpR(const YpR& ypr, const std::string& prefix) :
 
64
  AbstractParameterAliasable(ypr),
 
65
  AbstractSubstitutionModel(ypr),
 
66
  _pmodel(ypr._pmodel->clone()),
 
67
  _nestedPrefix(ypr.getNestedPrefix())
 
68
 
 
69
{
 
70
  _pmodel->setNamespace(prefix + _nestedPrefix);
 
71
}
 
72
 
 
73
YpR::YpR(const YpR& ypr) :
 
74
  AbstractParameterAliasable(ypr),
 
75
  AbstractSubstitutionModel(ypr),
 
76
  _pmodel(ypr._pmodel->clone()),
 
77
  _nestedPrefix(ypr.getNestedPrefix())
 
78
{}
 
79
 
 
80
void YpR::updateMatrices()
 
81
{
 
82
  updateMatrices(0,0,0,0,0,0,0,0);
 
83
}
 
84
 
 
85
 
 
86
void YpR::updateMatrices(double CgT, double cGA,
 
87
                         double TgC, double tGA,
 
88
                         double CaT, double cAG,
 
89
                         double TaC, double tAC)
 
90
{
 
91
  //  check_model(_pmodel);
 
92
 
 
93
  // Generator:
 
94
  const Alphabet* alph = _pmodel->getAlphabet();
 
95
  std::vector<int> l(4);
 
96
 
 
97
  l[0] = alph->charToInt("A");
 
98
  l[1] = alph->charToInt("G");
 
99
  l[2] = alph->charToInt("C");
 
100
  l[3] = alph->charToInt("T");
 
101
 
 
102
  unsigned int i,j,i1,i2,i3,j1,j2,j3;
 
103
 
 
104
  std::vector<double> a(4);  // a[A], a[G], a[C], a[T]
 
105
  std::vector<double> b(4);  // b[A], b[G], b[C], b[T]
 
106
 
 
107
  for (i = 0; i < 2; i++)
 
108
  {
 
109
    a[i] = _pmodel->Qij(l[1 - i],l[i]);
 
110
    b[i] = _pmodel->Qij(l[3 - i],l[i]);
 
111
    a[i + 2] = _pmodel->Qij(l[3 - i],l[i + 2]);
 
112
    b[i + 2] = _pmodel->Qij(l[1 - i],l[i + 2]);
 
113
  }
 
114
 
 
115
  // M_1
 
116
  RowMatrix<double> M1(3,3);
 
117
 
 
118
  M1(0,0) = 0;
 
119
  M1(0,1) = b[2];
 
120
  M1(0,2) = b[3];
 
121
  M1(1,0) = b[0] + b[1];
 
122
  M1(1,1) = 0;
 
123
  M1(1,2) = a[3];
 
124
  M1(2,0) = b[0] + b[1];
 
125
  M1(2,1) = a[2];
 
126
  M1(2,2) = 0;
 
127
 
 
128
  // M_2
 
129
  RowMatrix<double> M2(4,4);
 
130
 
 
131
  M2(0,0) = 0;
 
132
  M2(0,1) = a[1];
 
133
  M2(0,2) = b[2];
 
134
  M2(0,3) = b[3];
 
135
  M2(1,0) = a[0];
 
136
  M2(1,1) = 0;
 
137
  M2(1,2) = b[2];
 
138
  M2(1,3) = b[3];
 
139
  M2(2,0) = b[0];
 
140
  M2(2,1) = b[1];
 
141
  M2(2,2) = 0;
 
142
  M2(2,3) = a[3];
 
143
  M2(3,0) = b[0];
 
144
  M2(3,1) = b[1];
 
145
  M2(3,2) = a[2];
 
146
  M2(3,3) = 0;
 
147
 
 
148
  // M_3
 
149
  RowMatrix<double> M3(3,3);
 
150
 
 
151
  M3(0,0) = 0;
 
152
  M3(0,1) = a[1];
 
153
  M3(0,2) = b[2] + b[3];
 
154
  M3(1,0) = a[0];
 
155
  M3(1,1) = 0;
 
156
  M3(1,2) = b[2] + b[3];
 
157
  M3(2,0) = b[0];
 
158
  M3(2,1) = b[1];
 
159
  M3(2,2) = 0;
 
160
 
 
161
 
 
162
  for (i1 = 0; i1 < 3; i1++)
 
163
  {
 
164
    for (i2 = 0; i2 < 4; i2++)
 
165
    {
 
166
      for (i3 = 0; i3 < 3; i3++)
 
167
      {
 
168
        i = 12 * i1 + 3 * i2 + i3;
 
169
        for (j1 = 0; j1 < 3; j1++)
 
170
        {
 
171
          for (j2 = 0; j2 < 4; j2++)
 
172
          {
 
173
            for (j3 = 0; j3 < 3; j3++)
 
174
            {
 
175
              j = 12 * j1 + 3 * j2 + j3;
 
176
              if ((i1 == j1) && (i2 == j2))
 
177
                generator_(i,j) = M3(i3,j3);
 
178
              else if ((i1 == j1) && (i3 == j3))
 
179
                generator_(i,j) = M2(i2,j2);
 
180
              else if ((i2 == j2) && (i3 == j3))
 
181
                generator_(i,j) = M1(i1,j1);
 
182
              else
 
183
                generator_(i,j) = 0;
 
184
            }
 
185
          }
 
186
        }
 
187
      }
 
188
    }
 
189
  }
 
190
 
 
191
  // Introduction des dependances
 
192
 
 
193
  for (i3 = 0; i3 < 3; i3++)
 
194
  {
 
195
    generator_(15 + i3,12 + i3) += cGA * a[0]; // CG -> CA
 
196
    generator_(12 * i3 + 7,12 * i3 + 6) += cGA * a[0];
 
197
 
 
198
    generator_(15 + i3,27 + i3) += CgT * a[3]; // CG -> TG
 
199
    generator_(12 * i3 + 7,12 * i3 + 10) += CgT * a[3];
 
200
 
 
201
    generator_(27 + i3,24 + i3) += tGA * a[0]; // TG -> TA
 
202
    generator_(12 * i3 + 10,12 * i3 + 9) += tGA * a[0];
 
203
 
 
204
    generator_(27 + i3,15 + i3) += TgC * a[2]; // TG -> CG
 
205
    generator_(12 * i3 + 10,12 * i3 + 7) += TgC * a[2];
 
206
 
 
207
    generator_(12 + i3,24 + i3) += CaT * a[3]; // CA -> TA
 
208
    generator_(12 * i3 + 6,12 * i3 + 9) += CaT * a[3];
 
209
 
 
210
    generator_(12 + i3,15 + i3) += cAG * a[1]; // CA -> CG
 
211
    generator_(12 * i3 + 6,12 * i3 + 7) += cAG * a[1];
 
212
 
 
213
    generator_(24 + i3,27 + i3) += tAC * a[1]; // TA -> TG
 
214
    generator_(12 * i3 + 9,12 * i3 + 10) += tAC * a[1];
 
215
 
 
216
    generator_(24 + i3,12 + i3) += TaC * a[2]; // TA -> CA
 
217
    generator_(12 * i3 + 9,12 * i3 + 6) += TaC * a[2];
 
218
  }
 
219
 
 
220
  double x;
 
221
 
 
222
  for (i = 0; i < 36; i++)
 
223
  {
 
224
    x = 0;
 
225
    for (j = 0; j < 36; j++)
 
226
    {
 
227
      if (j != i)
 
228
        x += generator_(i,j);
 
229
    }
 
230
    generator_(i,i) = -x;
 
231
  }
 
232
 
 
233
  // calcul spectral
 
234
 
 
235
  EigenValue<double> ev(generator_);
 
236
  eigenValues_ = ev.getRealEigenValues();
 
237
  
 
238
  rightEigenVectors_ = ev.getV();
 
239
  MatrixTools::inv(rightEigenVectors_,leftEigenVectors_);
 
240
 
 
241
  iEigenValues_ = ev.getImagEigenValues();
 
242
 
 
243
  // frequence stationnaire
 
244
 
 
245
  x = 0;
 
246
  j = 0;
 
247
  while (j < 36){
 
248
    if (abs(eigenValues_[j]) < 0.000001 && abs(iEigenValues_[j]) < 0.000001) {
 
249
      eigenValues_[j]=0; //to avoid approximation problems in the future
 
250
      for (i = 0; i < 36; i++)
 
251
        {
 
252
          freq_[i] = leftEigenVectors_(j,i);
 
253
          x += freq_[i];
 
254
        }
 
255
      break;
 
256
    }
 
257
    j++;
 
258
  }
 
259
 
 
260
  for (i = 0; i < 36; i++)
 
261
    freq_[i] /= x;
 
262
 
 
263
  // mise a l'echelle
 
264
 
 
265
  x = 0;
 
266
  for (i1 = 0; i1 < 3; i1++)
 
267
  {
 
268
    for (i2 = 0; i2 < 4; i2++)
 
269
    {
 
270
      for (i3 = 0; i3 < 3; i3++)
 
271
      {
 
272
        i = 12 * i1 + 3 * i2 + i3;
 
273
        for (j2 = 0; j2 < 4; j2++)
 
274
        {
 
275
          if (j2 != i2)
 
276
          {
 
277
            j = 12 * i1 + 3 * j2 + i3;
 
278
            x += freq_[i] * generator_(i,j);
 
279
          }
 
280
        }
 
281
      }
 
282
    }
 
283
  }
 
284
 
 
285
  MatrixTools::scale(generator_,1 / x);
 
286
 
 
287
  for (i = 0; i < 36; i++)
 
288
    eigenValues_[i] /= x;
 
289
 
 
290
  isDiagonalizable_=true;
 
291
  for (i=0; i<size_ && isDiagonalizable_; i++)
 
292
    if (abs(iEigenValues_[i])> NumConstants::SMALL)
 
293
      isDiagonalizable_=false;  
 
294
}
 
295
 
 
296
 
 
297
void YpR::check_model(SubstitutionModel* const pm) const
 
298
throw (Exception)
 
299
{
 
300
  if (!pm)
 
301
    throw Exception("No Model ");
 
302
 
 
303
  const Alphabet* alph = pm->getAlphabet();
 
304
  if (alph->getAlphabetType() != "DNA alphabet")
 
305
    throw Exception("Need a DNA model");
 
306
 
 
307
  std::vector<int> l(4);
 
308
 
 
309
  l[0] = alph->charToInt("A");
 
310
  l[1] = alph->charToInt("G");
 
311
  l[2] = alph->charToInt("C");
 
312
  l[3] = alph->charToInt("T");
 
313
 
 
314
  // Check that the model is good for YpR
 
315
 
 
316
  int i;
 
317
 
 
318
  for (i = 0; i < 2; i++)
 
319
  {
 
320
    if (pm->Qij(l[2],l[i]) != pm->Qij(l[3],l[i]))
 
321
      throw Exception("Not R/Y Model " + pm->getName());
 
322
  }
 
323
  for (i = 2; i < 4; i++)
 
324
  {
 
325
    if (pm->Qij(l[0],l[i]) != pm->Qij(l[1],l[i]))
 
326
      throw Exception("Not R/Y Model " + pm->getName());
 
327
  }
 
328
}
 
329
 
 
330
void YpR::setNamespace(const std::string& prefix)
 
331
{
 
332
   AbstractSubstitutionModel::setNamespace(prefix);
 
333
  // We also need to update the namespace of the nested model:
 
334
  _pmodel->setNamespace(prefix + _nestedPrefix);
 
335
}
 
336
 
 
337
// ///////////////////////////////////////////////
 
338
// ///////////////////////////////////////////////
 
339
 
 
340
 
 
341
/******************************************************************************/
 
342
 
 
343
YpR_Sym::YpR_Sym(const RNY* alph,
 
344
                 SubstitutionModel* const pm,
 
345
                 double CgT, double TgC,
 
346
                 double CaT, double TaC) : AbstractParameterAliasable("YpR_Sym."), YpR(alph, pm,"YpR_Sym.")
 
347
{
 
348
  addParameter_(Parameter("YpR_Sym.rCgT", CgT, &Parameter::R_PLUS));
 
349
  addParameter_(Parameter("YpR_Sym.rTgC", TgC, &Parameter::R_PLUS));
 
350
  addParameter_(Parameter("YpR_Sym.rCaT", CaT, &Parameter::R_PLUS));
 
351
  addParameter_(Parameter("YpR_Sym.rTaC", TaC, &Parameter::R_PLUS));
 
352
 
 
353
  updateMatrices();
 
354
}
 
355
 
 
356
void YpR_Sym::updateMatrices()
 
357
{
 
358
   double rCgT = getParameterValue("rCgT");
 
359
   double rTgC = getParameterValue("rTgC");
 
360
   double rCaT = getParameterValue("rCaT");
 
361
   double rTaC = getParameterValue("rTaC");
 
362
 
 
363
   YpR::updateMatrices(rCgT, rCgT, rTgC, rTgC, rCaT, rCaT, rTaC, rTaC);
 
364
}
 
365
 
 
366
YpR_Sym::YpR_Sym(const YpR_Sym& ypr) : AbstractParameterAliasable(ypr), YpR(ypr,"YpR_Sym.")
 
367
{}
 
368
 
 
369
/******************************************************************************/
 
370
 
 
371
std::string YpR_Sym::getName() const
 
372
{
 
373
  return "YpR_Sym(model=" + _pmodel->getName()+")";
 
374
}
 
375
 
 
376
 
 
377
// ///////////////////////////////////////////////
 
378
// ///////////////////////////////////////////////
 
379
 
 
380
/******************************************************************************/
 
381
 
 
382
YpR_Gen::YpR_Gen(const RNY* alph,
 
383
                 SubstitutionModel* const pm,
 
384
                 double CgT, double cGA,
 
385
                 double TgC, double tGA,
 
386
                 double CaT, double cAG,
 
387
                 double TaC, double tAG) : AbstractParameterAliasable("YpR_Gen."), YpR(alph, pm,"YpR_Gen.")
 
388
{
 
389
  addParameter_(Parameter("YpR_Gen.rCgT", CgT, &Parameter::R_PLUS));
 
390
  addParameter_(Parameter("YpR_Gen.rcGA", cGA, &Parameter::R_PLUS));
 
391
  addParameter_(Parameter("YpR_Gen.rTgC", TgC, &Parameter::R_PLUS));
 
392
  addParameter_(Parameter("YpR_Gen.rtGA", tGA, &Parameter::R_PLUS));
 
393
  addParameter_(Parameter("YpR_Gen.rCaT", CaT, &Parameter::R_PLUS));
 
394
  addParameter_(Parameter("YpR_Gen.rcAG", cAG, &Parameter::R_PLUS));
 
395
  addParameter_(Parameter("YpR_Gen.rTaC", TaC, &Parameter::R_PLUS));
 
396
  addParameter_(Parameter("YpR_Gen.rtAG", tAG, &Parameter::R_PLUS));
 
397
 
 
398
  updateMatrices();
 
399
}
 
400
 
 
401
void YpR_Gen::updateMatrices()
 
402
{
 
403
   double rCgT = getParameterValue("rCgT");
 
404
   double rcGA = getParameterValue("rcGA");
 
405
   double rTgC = getParameterValue("rTgC");
 
406
   double rtGA = getParameterValue("rtGA");
 
407
   double rCaT = getParameterValue("rCaT");
 
408
   double rcAG = getParameterValue("rcAG");
 
409
   double rTaC = getParameterValue("rTaC");
 
410
   double rtAG = getParameterValue("rtAG");
 
411
 
 
412
   YpR::updateMatrices(rCgT, rcGA, rTgC, rtGA, rCaT, rcAG, rTaC, rtAG);
 
413
}
 
414
 
 
415
YpR_Gen::YpR_Gen(const YpR_Gen& ypr) : AbstractParameterAliasable(ypr), YpR(ypr,"YpR_Gen.")
 
416
{
 
417
  updateMatrices();
 
418
}
 
419
 
 
420
/******************************************************************************/
 
421
 
 
422
std::string YpR_Gen::getName() const
 
423
{
 
424
  return "YpR_Gen(model=" + _pmodel->getName()+")";
 
425
}
 
426
 
 
427