bpp-phyl  2.1.0
 All Classes Namespaces Files Functions Variables Friends Pages
MarkovModulatedSubstitutionModel.cpp
Go to the documentation of this file.
1 //
2 // File: MarkovModulatedSubstitutionModel.cpp
3 // Created by: Julien Dutheil
4 // Created on: Sat Aug 05 08:21 2006
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (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 "MarkovModulatedSubstitutionModel.h"
41 
42 #include <Bpp/Numeric/VectorTools.h>
43 #include <Bpp/Numeric/Matrix/MatrixTools.h>
44 #include <Bpp/Numeric/Matrix/EigenValue.h>
45 
46 using namespace bpp;
47 using namespace std;
48 
49 /******************************************************************************/
50 
51 MarkovModulatedSubstitutionModel::MarkovModulatedSubstitutionModel(
52  const MarkovModulatedSubstitutionModel& model) :
53  AbstractParameterAliasable(model),
54  model_ (dynamic_cast<ReversibleSubstitutionModel*>(model.model_->clone())),
55  nbStates_ (model.nbStates_),
56  nbRates_ (model.nbRates_),
57  rates_ (model.rates_),
58  ratesExchangeability_(model.ratesExchangeability_),
59  ratesFreq_ (model.ratesFreq_),
60  ratesGenerator_ (model.ratesGenerator_),
61  chars_ (model.chars_),
62  generator_ (model.generator_),
63  exchangeability_ (model.exchangeability_),
64  leftEigenVectors_ (model.leftEigenVectors_),
65  rightEigenVectors_ (model.rightEigenVectors_),
66  eigenValues_ (model.eigenValues_),
67  iEigenValues_ (model.iEigenValues_),
68  eigenDecompose_ (model.eigenDecompose_),
69  pijt_ (model.pijt_),
70  dpijt_ (model.dpijt_),
71  d2pijt_ (model.d2pijt_),
72  freq_ (model.freq_),
73  normalizeRateChanges_(model.normalizeRateChanges_),
74  nestedPrefix_ (model.nestedPrefix_)
75 {}
76 
77 MarkovModulatedSubstitutionModel& MarkovModulatedSubstitutionModel::operator=(
78  const MarkovModulatedSubstitutionModel& model)
79 {
80  AbstractParametrizable::operator=(model);
81  model_ = dynamic_cast<ReversibleSubstitutionModel*>(model.model_->clone());
82  nbStates_ = model.nbStates_;
83  nbRates_ = model.nbRates_;
84  rates_ = model.rates_;
85  ratesExchangeability_ = model.ratesExchangeability_;
86  ratesFreq_ = model.ratesFreq_;
87  ratesGenerator_ = model.ratesGenerator_;
88  chars_ = model.chars_;
89  generator_ = model.generator_;
90  exchangeability_ = model.exchangeability_;
91  leftEigenVectors_ = model.leftEigenVectors_;
92  rightEigenVectors_ = model.rightEigenVectors_;
93  eigenValues_ = model.eigenValues_;
94  iEigenValues_ = model.iEigenValues_;
95  eigenDecompose_ = model.eigenDecompose_;
96  pijt_ = model.pijt_;
97  dpijt_ = model.dpijt_;
98  d2pijt_ = model.d2pijt_;
99  freq_ = model.freq_;
100  normalizeRateChanges_ = model.normalizeRateChanges_;
101  nestedPrefix_ = model.nestedPrefix_;
102  return *this;
103 }
104 
105 /******************************************************************************/
106 
107 void MarkovModulatedSubstitutionModel::updateMatrices()
108 {
109  // ratesGenerator_ and rates_ must be initialized!
110  nbStates_ = model_->getNumberOfStates();
111  nbRates_ = rates_.getNumberOfColumns();
112  chars_ = VectorTools::rep(model_->getAlphabetChars(), nbRates_);
113  RowMatrix<double> Tmp1, Tmp2;
114  MatrixTools::diag(ratesFreq_, Tmp1);
115  MatrixTools::mult(ratesExchangeability_, Tmp1, ratesGenerator_);
116  MatrixTools::kroneckerMult(rates_, model_->getGenerator(), generator_);
117 
118  MatrixTools::MatrixTools::getId< RowMatrix<double> >(nbStates_, Tmp1);
119  MatrixTools::kroneckerMult(ratesGenerator_, Tmp1, Tmp2);
120  MatrixTools::add(generator_, Tmp2);
121 
122  MatrixTools::diag(1. / ratesFreq_, Tmp1);
123  MatrixTools::mult(rates_, Tmp1, Tmp2);
124  MatrixTools::kroneckerMult(Tmp2, model_->getExchangeabilityMatrix(), exchangeability_);
125 
126  MatrixTools::diag(1 / model_->getFrequencies(), Tmp1);
127  MatrixTools::kroneckerMult(ratesExchangeability_, Tmp1, Tmp2);
128  MatrixTools::add(exchangeability_, Tmp2);
129  freq_ = VectorTools::kroneckerMult(ratesFreq_, model_->getFrequencies());
130  if (normalizeRateChanges_)
131  {
132  // Normalization:
133  Vdouble Tmp;
134  MatrixTools::diag(generator_, Tmp);
135  double scale = -VectorTools::scalar<double, double>(Tmp, freq_);
136  MatrixTools::scale(generator_, 1. / scale);
137 
138  // Normalize exchangeability matrix too:
139  MatrixTools::scale(exchangeability_, 1. / scale);
140  }
141 
142  // Compute eigen values and vectors:
143  eigenValues_.resize(nbRates_ * nbStates_);
144  iEigenValues_.resize(nbRates_ * nbStates_);
145  rightEigenVectors_.resize(nbStates_ * nbRates_, nbStates_ * nbRates_);
146  pijt_.resize(nbStates_ * nbRates_, nbStates_ * nbRates_);
147  dpijt_.resize(nbStates_ * nbRates_, nbStates_ * nbRates_);
148  d2pijt_.resize(nbStates_ * nbRates_, nbStates_ * nbRates_);
149 
150  vector<double> modelEigenValues = model_->getEigenValues();
151  RowMatrix<double> modelRightEigenVectors = model_->getColumnRightEigenVectors();
152  for (unsigned int i = 0; i < nbStates_; i++)
153  {
154  RowMatrix<double> tmp = rates_;
155  MatrixTools::scale(tmp, modelEigenValues[i]);
156  MatrixTools::add(tmp, ratesGenerator_);
157  EigenValue<double> ev(tmp);
158  vector<double> values = ev.getRealEigenValues();
159  RowMatrix<double> vectors = ev.getV();
160  for (size_t j = 0; j < nbRates_; j++)
161  {
162  size_t c = i * nbRates_ + j; // Current eigen value index.
163  eigenValues_[c] = values[j];
164  // Compute the Kronecker product of the jth vector and the ith modelRightEigenVector.
165  for (unsigned int ii = 0; ii < nbRates_; ii++)
166  {
167  double vii = vectors(ii, j);
168  for (unsigned int jj = 0; jj < nbStates_; jj++)
169  {
170  rightEigenVectors_(ii * nbStates_ + jj, c) = vii * modelRightEigenVectors(jj, i);
171  }
172  }
173  }
174  }
175  // Now compute left eigen vectors by inversion:
176  MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
177 }
178 
179 /******************************************************************************/
180 
181 const Matrix<double>& MarkovModulatedSubstitutionModel::getPij_t(double t) const
182 {
183  if (t == 0)
184  MatrixTools::getId< RowMatrix<double> >(nbStates_ * nbRates_, pijt_);
185  else
186  MatrixTools::mult(rightEigenVectors_, VectorTools::exp(eigenValues_ * t), leftEigenVectors_, pijt_);
187  return pijt_;
188 }
189 
190 const Matrix<double>& MarkovModulatedSubstitutionModel::getdPij_dt(double t) const
191 {
192  MatrixTools::mult(rightEigenVectors_, eigenValues_ * VectorTools::exp(eigenValues_ * t), leftEigenVectors_, dpijt_);
193  return dpijt_;
194 }
195 
196 const Matrix<double>& MarkovModulatedSubstitutionModel::getd2Pij_dt2(double t) const
197 {
198  MatrixTools::mult(rightEigenVectors_, NumTools::sqr(eigenValues_) * VectorTools::exp(eigenValues_ * t), leftEigenVectors_, d2pijt_);
199  return d2pijt_;
200 }
201 
202 /******************************************************************************/
203 
204 double MarkovModulatedSubstitutionModel::getInitValue(size_t i, int state) const throw (IndexOutOfBoundsException, BadIntException)
205 {
206  if (i >= (nbStates_ * nbRates_))
207  throw IndexOutOfBoundsException("MarkovModulatedSubstitutionModel::getInitValue", i, 0, nbStates_ * nbRates_ - 1);
208  if (state < 0 || !model_->getAlphabet()->isIntInAlphabet(state))
209  throw BadIntException(state, "MarkovModulatedSubstitutionModel::getInitValue. Character " + model_->getAlphabet()->intToChar(state) + " is not allowed in model.");
210  vector<int> states = model_->getAlphabet()->getAlias(state);
211  for (size_t j = 0; j < states.size(); j++)
212  {
213  if (getAlphabetChar(i) == states[j])
214  return 1.;
215  }
216  return 0.;
217 }
218 
219 /******************************************************************************/
220 
221 void MarkovModulatedSubstitutionModel::setNamespace(const string& prefix)
222 {
223  AbstractParameterAliasable::setNamespace(prefix);
224  // We also need to update the namespace of the nested model:
225  model_->setNamespace(prefix + nestedPrefix_);
226 }
227 
228 /******************************************************************************/
229