← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Likelihood/AbstractNonHomogeneousTreeLikelihood.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/Likelihood/AbstractNonHomogeneousTreeLikelihood.cpp
6
6
//
7
7
 
8
8
/*
9
 
Copyright or © or Copr. CNRS, (November 16, 2004)
 
9
Copyright or © or Copr. Bio++ Development Team, (November 16, 2004)
10
10
 
11
11
This software is a computer program whose purpose is to provide classes
12
12
for phylogenetic data analysis.
59
59
  const Tree& tree,
60
60
  SubstitutionModelSet* modelSet,
61
61
  DiscreteDistribution* rDist,
62
 
  bool verbose)
 
62
  bool verbose,
 
63
  bool reparametrizeRoot)
63
64
  throw (Exception) :
64
65
  AbstractDiscreteRatesAcrossSitesTreeLikelihood(rDist, verbose),
65
66
  modelSet_(0),
77
78
  nbNodes_(),
78
79
  verbose_(),
79
80
  minimumBrLen_(),
 
81
  maximumBrLen_(),
80
82
  brLenConstraint_(0),
 
83
  reparametrizeRoot_(reparametrizeRoot),
81
84
  root1_(),
82
85
  root2_()
83
86
{
104
107
        nbNodes_(lik.nbNodes_),
105
108
  verbose_(lik.verbose_),
106
109
  minimumBrLen_(lik.minimumBrLen_),
 
110
  maximumBrLen_(lik.maximumBrLen_),
107
111
  brLenConstraint_(dynamic_cast<Constraint*>(lik.brLenConstraint_->clone())),
 
112
  reparametrizeRoot_(lik.reparametrizeRoot_),
108
113
  root1_(lik.root1_),
109
114
  root2_(lik.root2_)
110
115
{ 
124
129
    const AbstractNonHomogeneousTreeLikelihood& lik)
125
130
{
126
131
  AbstractDiscreteRatesAcrossSitesTreeLikelihood::operator=(lik);
127
 
  modelSet_        = lik.modelSet_;
128
 
  brLenParameters_ = lik.brLenParameters_;
129
 
  pxy_             = lik.pxy_;
130
 
  dpxy_            = lik.dpxy_;
131
 
  d2pxy_           = lik.d2pxy_;
132
 
  rootFreqs_       = lik.rootFreqs_;
133
 
  nodes_           = tree_->getNodes();
 
132
  modelSet_          = lik.modelSet_;
 
133
  brLenParameters_   = lik.brLenParameters_;
 
134
  pxy_               = lik.pxy_;
 
135
  dpxy_              = lik.dpxy_;
 
136
  d2pxy_             = lik.d2pxy_;
 
137
  rootFreqs_         = lik.rootFreqs_;
 
138
  nodes_             = tree_->getNodes();
134
139
  nodes_.pop_back(); //Remove the root node (the last added!).  
135
 
        nbSites_         = lik.nbSites_;
136
 
  nbDistinctSites_ = lik.nbDistinctSites_;
137
 
        nbClasses_       = lik.nbClasses_;
138
 
        nbStates_        = lik.nbStates_;
139
 
        nbNodes_         = lik.nbNodes_;
140
 
  verbose_         = lik.verbose_;
141
 
  minimumBrLen_    = lik.minimumBrLen_;
142
 
  brLenConstraint_.reset(dynamic_cast<Constraint*>(lik.brLenConstraint_->clone()));
143
 
  root1_           = lik.root1_;
144
 
  root2_           = lik.root2_;
 
140
        nbSites_           = lik.nbSites_;
 
141
  nbDistinctSites_   = lik.nbDistinctSites_;
 
142
        nbClasses_         = lik.nbClasses_;
 
143
        nbStates_          = lik.nbStates_;
 
144
        nbNodes_           = lik.nbNodes_;
 
145
  verbose_           = lik.verbose_;
 
146
  minimumBrLen_      = lik.minimumBrLen_;
 
147
  maximumBrLen_      = lik.maximumBrLen_;
 
148
  if (brLenConstraint_.get()) brLenConstraint_.release();
 
149
  brLenConstraint_.reset(lik.brLenConstraint_->clone());
 
150
  reparametrizeRoot_ = lik.reparametrizeRoot_;
 
151
  root1_             = lik.root1_;
 
152
  root2_             = lik.root2_;
145
153
  //Rebuild nodes index:
146
154
  for( unsigned int i = 0; i < nodes_.size(); i++)
147
155
  {
177
185
  verbose_ = verbose;
178
186
 
179
187
  minimumBrLen_ = 0.000001;
180
 
  brLenConstraint_.reset(new IncludingPositiveReal(minimumBrLen_));
 
188
  maximumBrLen_ = 10000;
 
189
  brLenConstraint_.reset(new IncludingInterval(minimumBrLen_, maximumBrLen_));
181
190
  setSubstitutionModelSet(modelSet);
182
191
}
183
192
 
270
279
 
271
280
ParameterList AbstractNonHomogeneousTreeLikelihood::getBranchLengthsParameters() const
272
281
{
273
 
  if(!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
 
282
  if (!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::getBranchLengthsParameters(). Object is not initialized.");
274
283
  return brLenParameters_.getCommonParametersWith(getParameters());
275
284
}
276
285
 
304
313
 
305
314
void AbstractNonHomogeneousTreeLikelihood::applyParameters() throw (Exception)
306
315
{
307
 
  if(!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
 
316
  if (!initialized_) throw Exception("AbstractBranchNonHomogeneousTreeLikelihood::applyParameters(). Object not initialized.");
308
317
  //Apply branch lengths:
309
 
  for(unsigned int i = 0; i < nbNodes_; i++)
 
318
  for (unsigned int i = 0; i < nbNodes_; i++)
310
319
  {
311
320
    int id = nodes_[i]->getId();
312
 
    if(id == root1_)
 
321
    if (reparametrizeRoot_ && id == root1_)
313
322
    {
314
 
      const Parameter * rootBrLen = &getParameter("BrLenRoot");
315
 
      const Parameter * rootPos = &getParameter("RootPosition");
 
323
      const Parameter* rootBrLen = &getParameter("BrLenRoot");
 
324
      const Parameter* rootPos = &getParameter("RootPosition");
316
325
      nodes_[i]->setDistanceToFather(rootBrLen->getValue() * rootPos->getValue());
317
326
    }
318
 
    else if(id == root2_)
 
327
    else if (reparametrizeRoot_ && id == root2_)
319
328
    {
320
 
      const Parameter * rootBrLen = &getParameter("BrLenRoot");
321
 
      const Parameter * rootPos = &getParameter("RootPosition");
 
329
      const Parameter* rootBrLen = &getParameter("BrLenRoot");
 
330
      const Parameter* rootPos = &getParameter("RootPosition");
322
331
      nodes_[i]->setDistanceToFather(rootBrLen->getValue() * (1. - rootPos->getValue()));
323
332
    }
324
333
    else
325
334
    {
326
 
      const Parameter * brLen = &getParameter(string("BrLen") + TextTools::toString(i));
327
 
      if(brLen != NULL) nodes_[i]->setDistanceToFather(brLen->getValue());
 
335
      const Parameter* brLen = &getParameter(string("BrLen") + TextTools::toString(i));
 
336
      if (brLen) nodes_[i]->setDistanceToFather(brLen->getValue());
328
337
    }
329
338
  }
330
339
  //Apply substitution model parameters:
356
365
        nodes_[i]->setDistanceToFather(minimumBrLen_);
357
366
        d = minimumBrLen_;
358
367
      }
 
368
      if (d > maximumBrLen_)
 
369
      {
 
370
        ApplicationTools::displayWarning("Branch length " + TextTools::toString(i) + " is too big: " + TextTools::toString(d) + ". Value is set to " + TextTools::toString(maximumBrLen_));
 
371
        nodes_[i]->setDistanceToFather(maximumBrLen_);
 
372
        d = maximumBrLen_;
 
373
      }
359
374
    }
360
 
    if (nodes_[i]->getId() == root1_)
 
375
    if (reparametrizeRoot_ && nodes_[i]->getId() == root1_)
361
376
      l1 = d;
362
 
    else if(nodes_[i]->getId() == root2_)
 
377
    else if (reparametrizeRoot_ && nodes_[i]->getId() == root2_)
363
378
      l2 = d;
364
379
    else
365
380
    {
366
381
      brLenParameters_.addParameter(Parameter("BrLen" + TextTools::toString(i), d, brLenConstraint_->clone(), true)); //Attach constraint to avoid clonage problems!
367
382
    }
368
383
  }
369
 
  brLenParameters_.addParameter(Parameter("BrLenRoot", l1 + l2, brLenConstraint_->clone(), true));
370
 
  brLenParameters_.addParameter(Parameter("RootPosition", l1 / (l1 + l2), &Parameter::PROP_CONSTRAINT_EX));
 
384
  if (reparametrizeRoot_) {
 
385
    brLenParameters_.addParameter(Parameter("BrLenRoot", l1 + l2, brLenConstraint_->clone(), true));
 
386
    brLenParameters_.addParameter(Parameter("RootPosition", l1 / (l1 + l2), &Parameter::PROP_CONSTRAINT_EX));
 
387
  }
371
388
}
372
389
 
373
390
/*******************************************************************************/