AbstractNonHomogeneousTreeLikelihood.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LEGACY_LIKELIHOOD_ABSTRACTNONHOMOGENEOUSTREELIKELIHOOD_H
#define BPP_PHYL_LEGACY_LIKELIHOOD_ABSTRACTNONHOMOGENEOUSTREELIKELIHOOD_H
 
 
#include "AbstractDiscreteRatesAcrossSitesTreeLikelihood.h"
#include "NonHomogeneousTreeLikelihood.h"
 
// From the STL:
#include <memory>
 
namespace bpp
{
class AbstractNonHomogeneousTreeLikelihood :
  public virtual NonHomogeneousTreeLikelihood,
  public AbstractDiscreteRatesAcrossSitesTreeLikelihood
{
public:
  class ConstNonHomogeneousSiteModelIterator :
    public ConstSiteModelIterator
  {
private:
    std::vector<ConstNoPartitionSiteModelDescription> siteModelDescriptions_;
    size_t index_;
    size_t nbModels_;
 
public:
    ConstNonHomogeneousSiteModelIterator(std::shared_ptr<const SubstitutionModelSet> modelSet) :
      siteModelDescriptions_(), index_(0), nbModels_(modelSet->getNumberOfModels())
    {
      for (size_t i = 0; i < nbModels_; ++i)
      {
        siteModelDescriptions_.push_back(ConstNoPartitionSiteModelDescription(modelSet->getModel(i), modelSet->getNodesWithModel(i)));
      }
    }
 
public:
    ConstSiteModelDescription* next()
    {
      if (!hasNext())
        throw Exception("AbstractNonHomogeneousTreeLikelihood::ConstHomogeneousSiteModelIterator::next(). No more site in the set.");
      return &siteModelDescriptions_[index_++];
    }
 
    bool hasNext() const { return index_ < nbModels_; }
  };
 
protected:
  std::shared_ptr<SubstitutionModelSet> modelSet_;
  ParameterList brLenParameters_;
 
  mutable std::map<int, VVVdouble> pxy_;
 
  mutable std::map<int, VVVdouble> dpxy_;
 
  mutable std::map<int, VVVdouble> d2pxy_;
 
  std::vector<double> rootFreqs_;
 
  std::vector<Node*> nodes_;
 
  mutable std::map<int, const Node*> idToNode_;
 
  // some values we'll need:
  size_t nbSites_,         // the number of sites in the container
      nbDistinctSites_, // the number of distinct sites
      nbClasses_,       // the number of rate classes
      nbStates_,        // the number of states in the alphabet
      nbNodes_;         // the number of nodes in the tree
 
  bool verbose_;
 
  double minimumBrLen_;
  double maximumBrLen_;
  std::shared_ptr<ConstraintInterface> brLenConstraint_;
 
  bool reparametrizeRoot_;
  int root1_, root2_;
 
public:
  AbstractNonHomogeneousTreeLikelihood(
      const Tree& tree,
      std::shared_ptr<SubstitutionModelSet> modelSet,
      std::shared_ptr<DiscreteDistributionInterface> rDist,
      bool verbose = true,
      bool reparametrizeRoot = true);
 
  AbstractNonHomogeneousTreeLikelihood(const AbstractNonHomogeneousTreeLikelihood& lik);
 
  AbstractNonHomogeneousTreeLikelihood& operator=(const AbstractNonHomogeneousTreeLikelihood& lik);
 
  virtual ~AbstractNonHomogeneousTreeLikelihood() {}
 
private:
  void init_(
      const Tree& tree,
      std::shared_ptr<SubstitutionModelSet> modelSet,
      std::shared_ptr<DiscreteDistributionInterface> rDist,
      bool verbose);
 
public:
  size_t getNumberOfStates() const override { return modelSet_->getNumberOfStates(); }
 
  const std::vector<int>& getAlphabetStates() const override { return modelSet_->getAlphabetStates(); }
 
  int getAlphabetStateAsInt(size_t i) const override { return modelSet_->getAlphabetStateAsInt(i); }
 
  std::string getAlphabetStateAsChar(size_t i) const override { return modelSet_->getAlphabetStateAsChar(i); }
 
  void initialize() override;
 
  ParameterList getBranchLengthsParameters() const override;
 
  ParameterList getSubstitutionModelParameters() const override;
 
  ParameterList getRateDistributionParameters() const override
  {
    return AbstractDiscreteRatesAcrossSitesTreeLikelihood::getRateDistributionParameters();
  }
 
  std::shared_ptr<const TransitionModelInterface> getModelForNode(int nodeId) const override
  {
    return modelSet_->getModelForNode(nodeId);
  }
 
  std::shared_ptr<TransitionModelInterface> getModelForNode(int nodeId) override
  {
    return modelSet_->getModelForNode(nodeId);
  }
 
  const std::vector<double>& getRootFrequencies(size_t siteIndex) const override { return rootFreqs_; }
 
  VVVdouble getTransitionProbabilitiesPerRateClass(int nodeId, size_t siteIndex) const override { return pxy_[nodeId]; }
 
  ConstBranchModelIterator* getNewBranchModelIterator(int nodeId) const override
  {
    return new ConstNoPartitionBranchModelIterator(modelSet_->getModelForNode(nodeId), nbDistinctSites_);
  }
 
  ConstSiteModelIterator* getNewSiteModelIterator(size_t siteIndex) const override
  {
    return new ConstNonHomogeneousSiteModelIterator(modelSet_);
  }
 
  std::shared_ptr<const SubstitutionModelSet> getSubstitutionModelSet() const override { return modelSet_; }
 
  const SubstitutionModelSet& substitutionModelSet() const override { return *modelSet_; }
 
  std::shared_ptr<SubstitutionModelSet> getSubstitutionModelSet() override { return modelSet_; }
 
  SubstitutionModelSet& substitutionModelSet() override { return *modelSet_; }
 
  void setSubstitutionModelSet(std::shared_ptr<SubstitutionModelSet> modelSet) override;
 
  ParameterList getRootFrequenciesParameters() const override
  {
    return modelSet_->getRootFrequenciesParameters();
  }
public:
  // Specific methods:
 
  virtual void initParameters();
 
  virtual void applyParameters();
 
  virtual void initBranchLengthsParameters(bool verbose = true);
 
  virtual void setMinimumBranchLength(double minimum)
  {
    if (minimum > maximumBrLen_)
      throw Exception("AbstractNonHomogeneousTreeLikelihood::setMinimumBranchLength. Minimum branch length should be lower than the maximum one: " + TextTools::toString(maximumBrLen_));
    minimumBrLen_ = minimum;
    brLenConstraint_ = std::make_shared<IntervalConstraint>(minimumBrLen_, maximumBrLen_, true, true);
    initBranchLengthsParameters();
  }
 
  virtual void setMaximumBranchLength(double maximum)
  {
    if (maximum < minimumBrLen_)
      throw Exception("AbstractNonHomogeneousTreeLikelihood::setMaximumBranchLength. Maximum branch length should be higher than the minimum one: " + TextTools::toString(minimumBrLen_));
    maximumBrLen_ = maximum;
    brLenConstraint_ = std::make_shared<IntervalConstraint>(minimumBrLen_, maximumBrLen_, true, true);
    initBranchLengthsParameters();
  }
 
  virtual double getMinimumBranchLength() const { return minimumBrLen_; }
  virtual double getMaximumBranchLength() const { return maximumBrLen_; }
 
protected:
  virtual void computeAllTransitionProbabilities();
  virtual void computeTransitionProbabilitiesForNode(const Node* node);
};
} // end of namespace bpp.
#endif // BPP_PHYL_LEGACY_LIKELIHOOD_ABSTRACTNONHOMOGENEOUSTREELIKELIHOOD_H