← 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/AbstractCodonSubstitutionModel.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/AbstractCodonSubstitutionModel.cpp
 
1
//
 
2
// File: CodonSubstitutionModel.cpp
 
3
// Created by:  Laurent Gueguen
 
4
// Created on: Feb 2009
 
5
//
 
6
 
 
7
/*
 
8
   Copyright or © or Copr. Bio++ Development Team, (November 16, 2004)
 
9
   This software is a computer program whose purpose is to provide classes
 
10
   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 only
 
20
   with a limited warranty  and the software's author,  the holder of the
 
21
   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 their
 
31
   requirements in conditions enabling the security of their systems and/or
 
32
   data to be ensured and,  more generally, to use and operate it in the
 
33
   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 "AbstractCodonSubstitutionModel.h"
 
40
 
 
41
using namespace bpp;
 
42
 
 
43
using namespace std;
 
44
 
 
45
/******************************************************************************/
 
46
 
 
47
AbstractCodonSubstitutionModel::AbstractCodonSubstitutionModel(
 
48
                                                               const CodonAlphabet* palph,
 
49
                                                               NucleotideSubstitutionModel* pmod,
 
50
                                                               const std::string& st,
 
51
                                                               bool paramRates) :
 
52
  AbstractParameterAliasable(st),
 
53
  AbstractSubstitutionModel(palph, st),
 
54
  AbstractWordSubstitutionModel(palph, st),
 
55
  hasParametrizedRates_(paramRates)
 
56
{
 
57
  enableEigenDecomposition(1);
 
58
 
 
59
  unsigned int i;
 
60
  for (i = 0; i < 3; i++)
 
61
  {
 
62
   VSubMod_.push_back(pmod);
 
63
   VnestedPrefix_.push_back(pmod->getNamespace());
 
64
  }
 
65
 
 
66
  pmod->setNamespace(st + "123_" + VnestedPrefix_[0]);
 
67
  pmod->enableEigenDecomposition(0);
 
68
  addParameters_(pmod->getParameters());
 
69
 
 
70
  Vrate_.resize(3);
 
71
  for (i = 0; i < 3; i++)
 
72
  {
 
73
    Vrate_[i] = 1.0 / 3;
 
74
  }
 
75
 
 
76
 
 
77
  if (hasParametrizedRates_){
 
78
 
 
79
  // relative rates
 
80
    for (i = 0; i < 2; i++)
 
81
      {
 
82
        addParameter_(Parameter(st+"relrate" + TextTools::toString(i + 1), 1.0 / (3 - i), &Parameter::PROP_CONSTRAINT_EX));
 
83
      }
 
84
    
 
85
  }
 
86
     
 
87
}
 
88
 
 
89
AbstractCodonSubstitutionModel::AbstractCodonSubstitutionModel(
 
90
                                                               const CodonAlphabet* palph,
 
91
                                                               NucleotideSubstitutionModel* pmod1,
 
92
                                                               NucleotideSubstitutionModel* pmod2,
 
93
                                                               NucleotideSubstitutionModel* pmod3,
 
94
                                                               const std::string& st,
 
95
                                                               bool paramRates) :
 
96
  AbstractParameterAliasable(st),
 
97
  AbstractSubstitutionModel(palph, st),
 
98
  AbstractWordSubstitutionModel(palph, st),
 
99
  hasParametrizedRates_(paramRates)
 
100
{
 
101
  unsigned int i;
 
102
  enableEigenDecomposition(1);
 
103
 
 
104
  if ((pmod1 == pmod2) || (pmod2 == pmod3) || (pmod1 == pmod3))
 
105
  {
 
106
    for (i = 0; i < 3; i++){
 
107
      VSubMod_.push_back(pmod1);
 
108
      VnestedPrefix_.push_back(pmod1->getNamespace());
 
109
    }
 
110
 
 
111
    pmod1->setNamespace(st + "123_" + VnestedPrefix_[0]);
 
112
    pmod1->enableEigenDecomposition(0);
 
113
    addParameters_(pmod1->getParameters());
 
114
  }
 
115
  else {
 
116
    VSubMod_.push_back(pmod1);
 
117
    VnestedPrefix_.push_back(pmod1->getNamespace());
 
118
    VSubMod_[0]->setNamespace(st + "1_" + VnestedPrefix_[0]);
 
119
    VSubMod_[0]->enableEigenDecomposition(0);
 
120
    addParameters_(pmod1->getParameters());
 
121
    
 
122
    VSubMod_.push_back(pmod2);
 
123
    VnestedPrefix_.push_back(pmod2->getNamespace());
 
124
    VSubMod_[1]->setNamespace(st + "2_" + VnestedPrefix_[1]);
 
125
    VSubMod_[1]->enableEigenDecomposition(0);
 
126
    addParameters_(pmod2->getParameters());
 
127
    
 
128
    VSubMod_.push_back(pmod3);
 
129
    VnestedPrefix_.push_back(pmod3->getNamespace());
 
130
    VSubMod_[2]->setNamespace(st + "3_" + VnestedPrefix_[2]);
 
131
    VSubMod_[2]->enableEigenDecomposition(0);
 
132
    addParameters_(pmod3->getParameters());
 
133
  }
 
134
 
 
135
  Vrate_.resize(3);
 
136
  for (i = 0; i < 3; i++)
 
137
  {
 
138
    Vrate_[i] = 1.0 / 3;
 
139
  }
 
140
 
 
141
  if (hasParametrizedRates_){
 
142
    // relative rates
 
143
    for (i = 0; i < 2; i++)
 
144
      {
 
145
        addParameter_(Parameter(st+"relrate" + TextTools::toString(i + 1), 1.0 / (3 - i), &Parameter::PROP_CONSTRAINT_EX));
 
146
      }
 
147
  }
 
148
}
 
149
 
 
150
void AbstractCodonSubstitutionModel::updateMatrices()
 
151
{
 
152
  if (hasParametrizedRates_){
 
153
    int i, nbmod = VSubMod_.size();
 
154
    double x;
 
155
    int k;
 
156
    for (k = nbmod - 1; k >= 0; k--)
 
157
      {
 
158
        x = 1.0;
 
159
        for (i = 0; i < k; i++)
 
160
          {
 
161
            x *= 1 - getParameterValue("relrate" + TextTools::toString(i + 1));
 
162
          }
 
163
        if (k != nbmod - 1)
 
164
          x *= getParameterValue("relrate" + TextTools::toString(k + 1));
 
165
        Vrate_[k] = x;
 
166
      }
 
167
  }
 
168
  
 
169
  AbstractWordSubstitutionModel::updateMatrices();
 
170
}
 
171
 
 
172
void AbstractCodonSubstitutionModel::completeMatrices()
 
173
{
 
174
  unsigned int i, j;
 
175
  unsigned int salph = getNumberOfStates();
 
176
  
 
177
  const CodonAlphabet* ca = dynamic_cast<const CodonAlphabet*>(alphabet_);
 
178
 
 
179
  for (i = 0; i < salph; i++)
 
180
    {
 
181
      for (j = 0; j < salph; j++)
 
182
        {
 
183
          if (ca->isStop(i) || ca->isStop(j))
 
184
            {
 
185
              generator_(i, j) = 0;
 
186
            }
 
187
          else
 
188
            generator_(i, j) *= getCodonsMulRate(i,j);
 
189
        }
 
190
    }  
 
191
}
 
192
 
 
193
  
 
194
 
 
195
 
 
196