NonHomogeneousSubstitutionProcess.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LIKELIHOOD_NONHOMOGENEOUSSUBSTITUTIONPROCESS_H
#define BPP_PHYL_LIKELIHOOD_NONHOMOGENEOUSSUBSTITUTIONPROCESS_H
 
 
#include "../Model/FrequencySet/FrequencySet.h"
#include "AbstractAutonomousSubstitutionProcess.h"
 
// From bpp-core:
#include <Bpp/Exceptions.h>
#include <Bpp/Numeric/Matrix/Matrix.h>
#include <Bpp/Numeric/AbstractParameterAliasable.h>
 
// From the STL:
#include <vector>
#include <map>
// #include <algorithm>
#include <memory>
#include <numeric>
 
namespace bpp
{
class NonHomogeneousSubstitutionProcess :
  public AbstractAutonomousSubstitutionProcess
{
private:
  std::vector<std::shared_ptr<BranchModelInterface>> modelSet_;
 
  std::shared_ptr<DiscreteDistributionInterface> rDist_;
 
  mutable std::map<unsigned int, size_t> nodeToModel_;
  mutable std::map<size_t, std::vector<unsigned int>> modelToNodes_;
 
  std::vector<ParameterList> modelParameters_;
 
public:
  NonHomogeneousSubstitutionProcess(
      std::shared_ptr<DiscreteDistributionInterface> rdist,
      std::shared_ptr<const PhyloTree> tree = 0,
      std::shared_ptr<FrequencySetInterface> rootFreqs = 0) :
    AbstractParameterAliasable(""),
    AbstractAutonomousSubstitutionProcess(tree, rootFreqs),
    modelSet_(),
    rDist_(rdist),
    nodeToModel_(),
    modelToNodes_(),
    modelParameters_()
  {
    if (rDist_)
      addParameters_(rDist_->getIndependentParameters());
  }
 
  NonHomogeneousSubstitutionProcess(
      std::shared_ptr<DiscreteDistributionInterface> rdist,
      std::shared_ptr<ParametrizablePhyloTree> tree,
      std::shared_ptr<FrequencySetInterface> rootFreqs = 0) :
    AbstractParameterAliasable(""),
    AbstractAutonomousSubstitutionProcess(tree, rootFreqs),
    modelSet_(),
    rDist_(rdist),
    nodeToModel_(),
    modelToNodes_(),
    modelParameters_()
  {
    if (rDist_)
      addParameters_(rDist_->getIndependentParameters());
  }
 
  NonHomogeneousSubstitutionProcess(const NonHomogeneousSubstitutionProcess& set);
 
  NonHomogeneousSubstitutionProcess& operator=(const NonHomogeneousSubstitutionProcess& set);
 
  virtual ~NonHomogeneousSubstitutionProcess()
  {
    clear();
  }
 
  NonHomogeneousSubstitutionProcess* clone() const override { return new NonHomogeneousSubstitutionProcess(*this); }
 
  void clear();
 
  void fireParameterChanged(const ParameterList& parameters) override;
 
  const StateMapInterface& stateMap() const override
  {
    if (modelSet_.size() == 0)
      throw Exception("NonHomogeneousSubstitutionProcess::getStateMap : no model associated");
    else
      return modelSet_[0]->stateMap();
  }
 
  std::shared_ptr<const StateMapInterface> getStateMap() const override
  {
    if (modelSet_.size() == 0)
      throw Exception("NonHomogeneousSubstitutionProcess::getStateMap : no model associated");
    else
      return modelSet_[0]->getStateMap();
  }
 
  size_t getNumberOfModels() const override { return modelSet_.size(); }
 
  bool hasMixedTransitionModel() const;
 
  void setModelScenario(std::shared_ptr<ModelScenario> modelscenario) override;
 
  std::vector<size_t> getModelNumbers() const override
  {
    std::vector<size_t> v(getNumberOfModels());
    std::iota(std::begin(v), std::end(v), 1);
    return v;
  }
 
  const BranchModelInterface& model(size_t n) const override
  {
    if ((n == 0) || (n > modelSet_.size())) throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::model().", 1, modelSet_.size(), n);
    return *modelSet_[n - 1];
  }
 
  std::shared_ptr<const BranchModelInterface> getModel(size_t n) const override
  {
    if ((n == 0) || (n > modelSet_.size())) throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::getModel().", 1, modelSet_.size(), n);
    return modelSet_[n - 1];
  }
 
  std::shared_ptr<const BranchModelInterface> getModel(size_t n)
  {
    if ((n == 0) || (n > modelSet_.size())) throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::getModel().", 1, modelSet_.size(), n);
    return modelSet_[n - 1];
  }
 
  size_t getModelNumberForNode(unsigned int nodeId) const override
  {
    const auto i = nodeToModel_.find(nodeId);
    if (i == nodeToModel_.end())
      throw Exception("NonHomogeneousSubstitutionProcess::getModelNumberForNode(). No model associated to node with id " + TextTools::toString(nodeId));
    return i->second + 1;
  }
 
  std::shared_ptr<const BranchModelInterface> getModelForNode(unsigned int nodeId) const override
  {
    std::map<unsigned int, size_t>::const_iterator i = nodeToModel_.find(nodeId);
    if (i == nodeToModel_.end())
      throw Exception("NonHomogeneousSubstitutionProcess::getModelForNode(). No model associated to node with id " + TextTools::toString(nodeId));
    return modelSet_[i->second];
  }
 
  const std::vector<unsigned int> getNodesWithModel(size_t i) const override
  {
    if (i > modelSet_.size()) throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::getNodesWithModel().", i, 0, modelSet_.size());
    return modelToNodes_[i - 1];
  }
 
  void addModel(std::shared_ptr<BranchModelInterface> model, const std::vector<unsigned int>& nodesId);
 
  void setModel(std::shared_ptr<BranchModelInterface> model, size_t modelIndex);
 
  void setModelToNode(size_t modelIndex, unsigned int nodeNumber);
 
  void listModelNames(std::ostream& out = std::cout) const;
 
  ParameterList getBranchLengthParameters(bool independent) const override
  {
    if (getParametrizablePhyloTree())
      return getParametrizablePhyloTree()->getParameters();
    else
      return ParameterList();
  }
 
  ParameterList getRateDistributionParameters(bool independent) const override
  {
    return rDist_.get() ? (independent ? rDist_->getIndependentParameters() : rDist_->getParameters()) : ParameterList();
  }
 
  const DiscreteDistributionInterface& rateDistribution() const override
  {
    if (!rDist_)
      throw NullPointerException("NonHomogeneousSubstitutionProcess::rateDistribution. No associated rate distribution.");
    return *rDist_;
  }
 
  DiscreteDistributionInterface& rateDistribution() override
  {
    if (!rDist_)
      throw NullPointerException("NonHomogeneousSubstitutionProcess::rateDistribution. No associated rate distribution.");
    return *rDist_;
  }
 
  std::shared_ptr<const DiscreteDistributionInterface> getRateDistribution() const override
  {
    return rDist_ ? rDist_ : nullptr;
  }
 
  std::shared_ptr<DiscreteDistributionInterface> getRateDistribution() override
  {
    return rDist_ ? rDist_ : nullptr;
  }
 
  ParameterList getSubstitutionModelParameters(bool independent) const override;
 
  bool isFullySetUp(bool throwEx = true) const
  {
    return checkOrphanNodes(throwEx)
           && checkUnknownNodes(throwEx);
  }
 
protected:
  bool checkOrphanNodes(bool throwEx) const;
 
  bool checkUnknownNodes(bool throwEx) const;
 
public:
  const std::vector<double>& getRootFrequencies() const override
  {
    if (!hasRootFrequencySet())
    {
      if (std::dynamic_pointer_cast<const TransitionModelInterface>(modelSet_[0]))
        return std::dynamic_pointer_cast<const TransitionModelInterface>(modelSet_[0])->getFrequencies();
      else
        throw Exception("NonHomogeneousSubstitutionProcess::getRootFrequencies not callable.");
    }
    else
      return rootFrequencySet().getFrequencies();
  }
 
  const BranchModelInterface& model(unsigned int nodeId, size_t classIndex) const override
  {
    return *modelSet_[nodeToModel_[nodeId]];
  }
 
  std::shared_ptr<const BranchModelInterface> getModel(unsigned int nodeId, size_t classIndex) const override
  {
    return modelSet_[nodeToModel_[nodeId]];
  }
 
  double getProbabilityForModel(size_t classIndex) const override
  {
    if (classIndex >= (rDist_ ? rDist_->getNumberOfCategories() : 1))
      throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::getProbabilityForModel.", classIndex, 0, rDist_->getNumberOfCategories());
    return rDist_ ? rDist_->getProbability(classIndex) : 1.;
  }
 
  Vdouble getClassProbabilities() const override
  {
    Vdouble vProb;
 
    if (!rDist_)
      vProb.push_back(1.);
    else
      for (size_t i = 0; i < rDist_->getNumberOfCategories(); i++)
      {
        vProb.push_back(rDist_->getProbability(i));
      }
 
    return vProb;
  }
 
  double getRateForModel(size_t classIndex) const override
  {
    if (classIndex >= (rDist_ ? rDist_->getNumberOfCategories() : 1))
      throw IndexOutOfBoundsException("NonHomogeneousSubstitutionProcess::getRateForModel.", classIndex, 0, (rDist_ ? rDist_->getNumberOfCategories() : 1));
    return rDist_ ? rDist_->getCategory(classIndex) : 1;
  }
 
  static std::unique_ptr<AutonomousSubstitutionProcessInterface> createHomogeneousSubstitutionProcess(
      std::shared_ptr<BranchModelInterface> model,
      std::shared_ptr<DiscreteDistributionInterface> rdist,
      std::shared_ptr<PhyloTree> tree,
      std::shared_ptr<FrequencySetInterface> rootFreqs = 0,
      std::shared_ptr<ModelScenario> scenario = 0
      );
 
  static std::unique_ptr<NonHomogeneousSubstitutionProcess> createNonHomogeneousSubstitutionProcess(
      std::shared_ptr<BranchModelInterface> model,
      std::shared_ptr<DiscreteDistributionInterface> rdist,
      std::shared_ptr<PhyloTree> tree,
      std::shared_ptr<FrequencySetInterface> rootFreqs,
      const std::vector<std::string>& globalParameterNames,
      std::shared_ptr<ModelScenario> scenario = 0
      );
 
 
};
} // end of namespace bpp.
#endif // BPP_PHYL_LIKELIHOOD_NONHOMOGENEOUSSUBSTITUTIONPROCESS_H