AbstractSubstitutionModel.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_MODEL_ABSTRACTSUBSTITUTIONMODEL_H
#define BPP_PHYL_MODEL_ABSTRACTSUBSTITUTIONMODEL_H
 
#include <Bpp/Numeric/AbstractParameterAliasable.h>
#include <Bpp/Numeric/VectorTools.h>
#include <memory>
 
#include "SubstitutionModel.h"
 
namespace bpp
{
class AbstractLkTransitionModel :
  public virtual TransitionModelInterface
{
private:
  mutable Eigen::VectorXd lik_;
 
public:
  AbstractLkTransitionModel() : lik_() {}
  virtual ~AbstractLkTransitionModel() {}
 
public:
  const Eigen::VectorXd& Lik_t(const Eigen::VectorXd& values, double t) const override
  {
    lik_ = Eigen::VectorXd::Zero(values.size());
 
    const auto& pij = getPij_t(t);
 
    for (auto i = 0; i < values.size(); ++i)
    {
      for (auto j = 0; j < values.size(); ++j)
      {
        lik_(i) += pij(size_t(i), size_t(j)) * values(j);
      }
    }
 
    return lik_;
  }
 
  const Eigen::VectorXd& dLik_dt(const Eigen::VectorXd& values, double t) const override
  {
    lik_ = Eigen::VectorXd::Zero(values.size());
 
    const auto& pij = getdPij_dt(t);
 
    for (auto i = 0; i < values.size(); ++i)
    {
      for (auto j = 0; j < values.size(); ++j)
      {
        lik_(i) += pij(size_t(i), size_t(j)) * values(j);
      }
    }
 
    return lik_;
  }
 
  const Eigen::VectorXd& d2Lik_dt2(const Eigen::VectorXd& values, double t) const override
  {
    lik_ = Eigen::VectorXd::Zero(values.size());
 
    const auto& pij = getd2Pij_dt2(t);
 
    for (auto i = 0; i < values.size(); ++i)
    {
      for (auto j = 0; j < values.size(); ++j)
      {
        lik_(i) += pij(size_t(i), size_t(j)) * values(j);
      }
    }
 
    return lik_;
  }
 
  friend class AbstractTotallyWrappedTransitionModel;
  friend class AbstractFromSubstitutionModelTransitionModel;
  friend class InMixedSubstitutionModel;
};
 
class AbstractTransitionModel :
  public virtual AbstractLkTransitionModel,
  public virtual AbstractParameterAliasable
{
protected:
  std::shared_ptr<const Alphabet> alphabet_;
 
  std::shared_ptr<const StateMapInterface> stateMap_;
 
  size_t size_;
 
  double rate_;
 
  Vdouble freq_;
 
  mutable RowMatrix<double> pijt_;
  mutable RowMatrix<double> dpijt_;
  mutable RowMatrix<double> d2pijt_;
 
  short verboseLevel_;
 
public:
  AbstractTransitionModel(
      std::shared_ptr<const Alphabet> alpha,
      std::shared_ptr<const StateMapInterface> stateMap,
      const std::string& prefix);
 
  AbstractTransitionModel(const AbstractTransitionModel& model) = default;
 
  AbstractTransitionModel& operator=(const AbstractTransitionModel& model) = default;
 
  virtual ~AbstractTransitionModel() {}
 
public:
  const Alphabet& alphabet() const override { return *alphabet_; }
 
  std::shared_ptr<const Alphabet> getAlphabet() const override { return alphabet_; }
 
  const StateMapInterface& stateMap() const override { return *stateMap_; }
 
  std::shared_ptr<const StateMapInterface> getStateMap() const override { return stateMap_; }
 
  size_t getNumberOfStates() const override { return stateMap_->getNumberOfModelStates();}
 
  const std::vector<int>& getAlphabetStates() const override { return stateMap_->getAlphabetStates(); }
 
  std::string getAlphabetStateAsChar(size_t index) const override { return stateMap_->getAlphabetStateAsChar(index); }
 
  int getAlphabetStateAsInt(size_t index) const override { return stateMap_->getAlphabetStateAsInt(index); }
 
  std::vector<size_t> getModelStates(int code) const override { return stateMap_->getModelStates(code); }
 
  std::vector<size_t> getModelStates(const std::string& code) const override { return stateMap_->getModelStates(code); }
 
  const Vdouble& getFrequencies() const override { return freq_; }
 
  bool computeFrequencies() const override { return false; }
 
  virtual const Matrix<double>& getPij_t(double t) const override = 0;
  virtual const Matrix<double>& getdPij_dt(double t) const override = 0;
  virtual const Matrix<double>& getd2Pij_dt2(double t) const override = 0;
 
  virtual double freq(size_t i) const override { return freq_[i]; }
 
  virtual double Pij_t    (size_t i, size_t j, double t) const override { return getPij_t(t) (i, j); }
  virtual double dPij_dt  (size_t i, size_t j, double t) const override { return getdPij_dt(t) (i, j); }
  virtual double d2Pij_dt2(size_t i, size_t j, double t) const override { return getd2Pij_dt2(t) (i, j); }
 
  double getInitValue(size_t i, int state) const override;
 
  void setFreqFromData(const SequenceDataInterface& data, double pseudoCount = 0) override;
 
  virtual void setFreq(std::map<int, double>& freqs) override;
 
  const FrequencySetInterface& frequencySet() const override
  {
    throw Exception("TransitionModel::frequencySet(). No associated FrequencySet object.");
  }
 
  virtual void fireParameterChanged(const ParameterList& parameters) override
  {
    AbstractParameterAliasable::fireParameterChanged(parameters);
 
    if (parameters.hasParameter(getNamespace() + "rate"))
    {
      rate_ = parameters.getParameterValue(getNamespace() + "rate");
 
      if (parameters.size() != 1)
        updateMatrices_();
    }
    else
      updateMatrices_();
  }
 
  void addRateParameter() override;
 
  void setVerboseLevel(short level) { verboseLevel_ = level; }
 
  short verboseLevel() const { return verboseLevel_; }
 
protected:
  virtual void updateMatrices_() = 0;
 
  /*
   * @brief : To update the eq freq
   *
   */
  Vdouble& getFrequencies_() override
  {
    return freq_;
  }
 
public:
  virtual double getRate() const override;
 
  virtual void setRate(double rate) override;
};
 
 
class AbstractSubstitutionModel :
  public AbstractTransitionModel,
  public virtual SubstitutionModelInterface
{
protected:
  bool isScalable_;
 
  RowMatrix<double> generator_;
 
  bool computeFreq_;
 
  RowMatrix<double> exchangeability_;
 
  bool eigenDecompose_;
 
  Vdouble eigenValues_;
 
  Vdouble iEigenValues_;
 
  bool isDiagonalizable_;
 
  RowMatrix<double> rightEigenVectors_;
 
  bool isNonSingular_;
 
  RowMatrix<double> leftEigenVectors_;
 
  std::vector< RowMatrix<double>> vPowGen_;
 
  mutable RowMatrix<double> tmpMat_;
 
public:
  AbstractSubstitutionModel(
      std::shared_ptr<const Alphabet> alpha,
      std::shared_ptr<const StateMapInterface> stateMap,
      const std::string& prefix);
 
  AbstractSubstitutionModel(const AbstractSubstitutionModel& model) :
    AbstractTransitionModel(model),
    isScalable_(model.isScalable_),
    generator_(model.generator_),
    computeFreq_(model.computeFreq_),
    exchangeability_(model.exchangeability_),
    eigenDecompose_(model.eigenDecompose_),
    eigenValues_(model.eigenValues_),
    iEigenValues_(model.iEigenValues_),
    isDiagonalizable_(model.isDiagonalizable_),
    rightEigenVectors_(model.rightEigenVectors_),
    isNonSingular_(model.isNonSingular_),
    leftEigenVectors_(model.leftEigenVectors_),
    vPowGen_(model.vPowGen_),
    tmpMat_(model.tmpMat_)
  {}
 
  AbstractSubstitutionModel& operator=(const AbstractSubstitutionModel& model)
  {
    AbstractTransitionModel::operator=(model);
    isScalable_        = model.isScalable_;
    generator_         = model.generator_;
    computeFreq_       = model.computeFreq_;
    exchangeability_   = model.exchangeability_;
    eigenDecompose_    = model.eigenDecompose_;
    eigenValues_       = model.eigenValues_;
    iEigenValues_      = model.iEigenValues_;
    isDiagonalizable_  = model.isDiagonalizable_;
    rightEigenVectors_ = model.rightEigenVectors_;
    isNonSingular_     = model.isNonSingular_;
    leftEigenVectors_  = model.leftEigenVectors_;
    vPowGen_           = model.vPowGen_;
    tmpMat_            = model.tmpMat_;
    return *this;
  }
 
  virtual ~AbstractSubstitutionModel() {}
 
public:
  bool computeFrequencies() const { return computeFreq_; }
 
  void computeFrequencies(bool yn) { computeFreq_ = yn; }
 
  const Matrix<double>& generator() const { return generator_; }
 
  const Matrix<double>& exchangeabilityMatrix() const { return exchangeability_; }
 
  const Matrix<double>& getPij_t(double t) const;
  const Matrix<double>& getdPij_dt(double t) const;
  const Matrix<double>& getd2Pij_dt2(double t) const;
 
  double Sij(size_t i, size_t j) const { return exchangeability_(i, j); }
 
  const Vdouble& getEigenValues() const { return eigenValues_; }
 
  const Vdouble& getIEigenValues() const { return iEigenValues_; }
 
  bool isDiagonalizable() const { return isDiagonalizable_; }
 
  bool isNonSingular() const { return isNonSingular_; }
 
  const Matrix<double>& getRowLeftEigenVectors() const { return leftEigenVectors_; }
 
  const Matrix<double>& getColumnRightEigenVectors() const { return rightEigenVectors_; }
 
  virtual double Qij(size_t i, size_t j) const { return generator_(i, j); }
 
  void enableEigenDecomposition(bool yn) { eigenDecompose_ = yn; }
 
  bool enableEigenDecomposition() { return eigenDecompose_; }
 
protected:
  virtual void updateMatrices_();
 
public:
  void setScalable(bool scalable)
  {
    isScalable_ = scalable;
  }
 
  virtual bool isScalable() const
  {
    return isScalable_;
  }
 
  double getScale() const;
 
  void setScale(double scale);
 
  void normalize();
 
  void setDiagonal();
 
  friend class OneChangeRegisterTransitionModel;
};
 
 
class AbstractReversibleSubstitutionModel :
  public AbstractSubstitutionModel,
  public virtual ReversibleSubstitutionModelInterface
{
public:
  AbstractReversibleSubstitutionModel(
      std::shared_ptr<const Alphabet> alpha,
      std::shared_ptr<const StateMapInterface> stateMap,
      const std::string& prefix) :
    AbstractSubstitutionModel(alpha, stateMap, prefix)
  {
    isDiagonalizable_ = true;
    isNonSingular_    = true;
    computeFreq_      = false;
 
    // to ensure non null freq_ at construction
    for (auto& fr : freq_)
    {
      fr = 1.0 / static_cast<double>(size_);
    }
 
  }
 
  virtual ~AbstractReversibleSubstitutionModel() {}
 
  virtual AbstractReversibleSubstitutionModel* clone() const override = 0;
 
protected:
  virtual void updateMatrices_() override;
};
} // end of namespace bpp.
#endif // BPP_PHYL_MODEL_ABSTRACTSUBSTITUTIONMODEL_H