SubstitutionProcessCollection.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LIKELIHOOD_SUBSTITUTIONPROCESSCOLLECTION_H
#define BPP_PHYL_LIKELIHOOD_SUBSTITUTIONPROCESSCOLLECTION_H
 
// From bpp-core:
#include <Bpp/Exceptions.h>
#include <Bpp/Numeric/ParametrizableCollection.h>
#include <Bpp/Numeric/Prob/DiscreteDistribution.h>
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Numeric/VectorTools.h>
 
// From bpp-seq:
#include <Bpp/Seq/Alphabet/Alphabet.h>
#include <Bpp/Seq/Alphabet/NucleicAlphabet.h>
 
#include "../Model/FrequencySet/FrequencySet.h"
#include "../Model/SubstitutionModel.h"
#include "ParametrizablePhyloTree.h"
#include "SubstitutionProcess.h"
#include "SubstitutionProcessCollectionMember.h"
 
 
namespace bpp
{
class SubstitutionProcessCollection :
  public virtual Clonable,
  public AbstractParameterAliasable
{
private:
  ParametrizableCollection<BranchModelInterface> modelColl_;
 
  std::map<size_t, std::vector<size_t>> mModelToSubPro_;
 
  ParametrizableCollection<FrequencySetInterface> freqColl_;
 
  std::map<size_t, std::vector<size_t>> mFreqToSubPro_;
 
  ParametrizableCollection<DiscreteDistributionInterface> distColl_;
 
  std::map<size_t, std::vector<size_t>> mVConstDist_;
 
  std::map<size_t, std::vector<size_t>> mDistToSubPro_;
 
  ParametrizableCollection<ParametrizablePhyloTree> treeColl_;
 
  std::map<size_t, std::vector<size_t>> mTreeToSubPro_;
 
  std::map<size_t, std::shared_ptr<ModelScenario>> mModelScenario_;
 
  std::map<size_t, std::shared_ptr<SubstitutionProcessCollectionMember>> mSubProcess_; // Need a specific deleter because the destructor is private.
 
public:
  SubstitutionProcessCollection() :
    AbstractParameterAliasable(""),
    modelColl_(),
    mModelToSubPro_(),
    freqColl_(),
    mFreqToSubPro_(),
    distColl_(),
    mVConstDist_(),
    mDistToSubPro_(),
    treeColl_(),
    mTreeToSubPro_(),
    mSubProcess_()
  {}
 
 
  SubstitutionProcessCollection(const SubstitutionProcessCollection& set);
 
  SubstitutionProcessCollection& operator=(const SubstitutionProcessCollection& set);
 
  virtual ~SubstitutionProcessCollection()
  {
    clear();
  }
 
  SubstitutionProcessCollection* clone() const { return new SubstitutionProcessCollection(*this); }
 
 
  void clear();
 
  void addParametrizable(std::shared_ptr<Parametrizable> parametrizable, size_t parametrizableIndex, bool withParameters = true);
 
  void replaceParametrizable(std::shared_ptr<Parametrizable> parametrizable, size_t parametrizableIndex, bool withParameters = true);
 
  void addModel(std::shared_ptr<BranchModelInterface> model, size_t modelIndex)
  {
    addParametrizable(model, modelIndex);
  }
 
  void replaceModel(std::shared_ptr<BranchModelInterface> model, size_t modelIndex)
  {
    replaceParametrizable(model, modelIndex);
  }
 
  void addFrequencies(std::shared_ptr<FrequencySetInterface> frequencies, size_t frequenciesIndex)
  {
    addParametrizable(frequencies, frequenciesIndex);
  }
 
  void replaceFrequencies(std::shared_ptr<FrequencySetInterface> frequencies, size_t frequenciesIndex)
  {
    addParametrizable(frequencies, frequenciesIndex);
  }
 
  void addDistribution(std::shared_ptr<DiscreteDistributionInterface> distribution, size_t distributionIndex)
  {
    addParametrizable(distribution, distributionIndex, (distributionIndex < 10000));
 
    if (distributionIndex >= 10000)
      mVConstDist_[distributionIndex / 10000 - 1].push_back(distributionIndex % 10000);
  }
 
 
  void addTree(std::shared_ptr<ParametrizablePhyloTree> tree, size_t treeIndex)
  {
    addParametrizable(tree, treeIndex);
  }
 
  void replaceTree(std::shared_ptr<ParametrizablePhyloTree> tree, size_t treeIndex)
  {
    replaceParametrizable(tree, treeIndex);
  }
 
  void addScenario(std::shared_ptr<ModelScenario> scen, size_t scenIndex)
  {
    mModelScenario_[scenIndex] = scen;
  }
 
  BranchModelInterface& model(size_t modelIndex)
  {
    return dynamic_cast<BranchModelInterface&>(*modelColl_[modelIndex]);
  }
 
  const BranchModelInterface& model(size_t modelIndex) const
  {
    return dynamic_cast<const BranchModelInterface&>(*modelColl_[modelIndex]);
  }
 
  std::shared_ptr<BranchModelInterface> getModel(size_t modelIndex)
  {
    return std::dynamic_pointer_cast<BranchModelInterface>(modelColl_[modelIndex]);
  }
 
  std::shared_ptr<const BranchModelInterface> getModel(size_t modelIndex) const
  {
    return std::dynamic_pointer_cast<const BranchModelInterface>(modelColl_[modelIndex]);
  }
 
  size_t getModelIndex(std::shared_ptr<BranchModelInterface> model) const
  {
    if (modelColl_.hasObject(model))
      return modelColl_.getFirstKey(model);
    else
      return 0;
  }
 
  bool hasModel(std::shared_ptr<BranchModelInterface> model) const
  {
    return modelColl_.hasObject(model);
  }
 
  FrequencySetInterface& frequencySet(size_t frequenciesIndex)
  {
    return *freqColl_[frequenciesIndex];
  }
 
  const FrequencySetInterface& frequencySet(size_t frequenciesIndex) const
  {
    return *(std::dynamic_pointer_cast<const FrequencySetInterface>(freqColl_[frequenciesIndex]));
  }
 
  std::shared_ptr<FrequencySetInterface> getFrequencySet(size_t frequenciesIndex)
  {
    return freqColl_[frequenciesIndex];
  }
 
  std::shared_ptr<const FrequencySetInterface> getFrequencySet(size_t frequenciesIndex) const
  {
    return freqColl_[frequenciesIndex];
  }
 
  std::shared_ptr<DiscreteDistributionInterface> getRateDistribution(size_t distributionIndex)
  {
    return distColl_[distributionIndex];
  }
 
  std::shared_ptr<const DiscreteDistributionInterface> getRateDistribution(size_t distributionIndex) const
  {
    return distColl_[distributionIndex];
  }
 
  DiscreteDistributionInterface& rateDistribution(size_t distributionIndex)
  {
    return *distColl_[distributionIndex];
  }
 
  const DiscreteDistributionInterface& rateDistribution(size_t distributionIndex) const
  {
    return *distColl_[distributionIndex];
  }
 
  ParametrizablePhyloTree& tree(size_t treeIndex)
  {
    return dynamic_cast<ParametrizablePhyloTree&>(*treeColl_[treeIndex]);
  }
 
  const ParametrizablePhyloTree& tree(size_t treeIndex) const
  {
    return dynamic_cast<const ParametrizablePhyloTree&>(*treeColl_[treeIndex]);
  }
 
  std::shared_ptr<ParametrizablePhyloTree> getTree(size_t treeIndex)
  {
    return std::dynamic_pointer_cast<ParametrizablePhyloTree>(treeColl_[treeIndex]);
  }
 
  std::shared_ptr<const ParametrizablePhyloTree> getTree(size_t treeIndex) const
  {
    return std::dynamic_pointer_cast<const ParametrizablePhyloTree>(treeColl_[treeIndex]);
  }
 
  bool hasModelScenario(size_t numPath) const
  {
    return mModelScenario_.find(numPath) != mModelScenario_.end();
  }
 
  std::shared_ptr<const ModelScenario> getModelScenario(size_t numPath) const
  {
    if (!hasModelScenario(numPath))
      return 0;
    else
      return mModelScenario_.at(numPath);
  }
 
  std::shared_ptr<ModelScenario> getModelScenario(size_t numPath)
  {
    if (!hasModelScenario(numPath))
      return 0;
    else
      return mModelScenario_.at(numPath);
  }
 
  std::vector<size_t> getScenarioNumbers() const
  {
    std::vector<size_t> vkeys;
 
    for (const auto& it:mModelScenario_)
    {
      vkeys.push_back(it.first);
    }
 
    return vkeys;
  }
 
 
  std::vector<size_t> getModelNumbers() const
  {
    return modelColl_.keys();
  }
 
  std::vector<size_t> getFrequenciesNumbers() const
  {
    return freqColl_.keys();
  }
 
  std::vector<size_t> getRateDistributionNumbers() const
  {
    return distColl_.keys();
  }
 
  std::vector<size_t> getTreeNumbers() const
  {
    return treeColl_.keys();
  }
 
  std::vector<size_t> getSubstitutionProcessNumbers() const
  {
    std::vector<size_t> vn;
 
    for (auto it = mSubProcess_.begin(); it != mSubProcess_.end(); it++)
    {
      vn.push_back(it->first);
    }
 
    return vn;
  }
 
  bool hasModelNumber(size_t n) const
  {
    return modelColl_.hasObject(n);
  }
 
  bool hasFrequenciesNumber(size_t n) const
  {
    return freqColl_.hasObject(n);
  }
 
  bool hasDistributionNumber(size_t n) const
  {
    return distColl_.hasObject(n);
  }
 
  bool hasTreeNumber(size_t n) const
  {
    return treeColl_.hasObject(n);
  }
 
  bool hasSubstitutionProcessNumber(size_t n) const
  {
    return mSubProcess_.find(n) != mSubProcess_.end();
  }
 
  void fireParameterChanged(const ParameterList& parameters);
 
  void setNamespace(const std::string& prefix);
 
  void aliasParameters(const std::string& p1, const std::string& p2);
 
  void unaliasParameters(const std::string& p1, const std::string& p2);
 
  void aliasParameters(std::map<std::string, std::string>& unparsedParams, bool verbose);
 
  /*
   * @brief Method to add a SubstitutionProcess.
   *
   * @param nProc the number of the Substitution process
   * @param mModBr a map associating numbers of models (from the collection) and numbers of branches
   * @param nTree the number of a Tree (from the collection)
   * @param nRate the number of a Distribution Rate (from the collection)
   * @param nFreq the number of a FrequencySet for the root (from the collection)
   *
   * @throw an Exception if the built SubstitutionModelSet is not complete or well built.
   *
   */
  void addSubstitutionProcess(size_t nProc, std::map<size_t, std::vector<unsigned int>> mModBr, size_t nTree, size_t nRate, size_t nFreq);
 
  /*
   * @brief Method to add a SubstitutionProcess.
   *
   * @param nProc the number of the Substitution process
   * @param mModBr a map associating numbers of models (from the collection) and numbers of branches
   * @param nTree the number of a Tree (from the collection)
   * @param nRate the number of a Distribution Rate (from the collection)
   *
   * @throw an Exception if the built SubstitutionModelSet is not complete or well built.
   *
   */
  void addSubstitutionProcess(size_t nProc, std::map<size_t, std::vector<unsigned int>> mModBr, size_t nTree, size_t nRate);
 
  /*
   * @brief Method to add a one per branch SubstitutionProcess. A new
   * model is created for each branch.
   *
   * @param nProc the number of the Substitution process
   * @param nMod  the number of Transition Model which is cloned on
   *        all branches.
   * @param nTree the number of a Tree (from the collection)
   * @param nRate the number of a Distribution Rate (from the
   *        collection)
   * @param nFreq the number of a FrequencySet for the root (from
   *        the collection)
   * @param sharedParameterNames the vector of the names of parameters
   *        of the model that are shared among all the models.
   *
   */
  void addOnePerBranchSubstitutionProcess(size_t nProc, size_t nMod, size_t nTree, size_t nRate, size_t nFreq, const std::vector<std::string>& sharedParameterNames);
 
  /*
   * @brief Method to add a one per branch SubstitutionProcess. A new
   * model is created for each branch.
   *
   * @param nProc the number of the Substitution process
   * @param nMod  the number of Transition Model which is cloned on
   *        all branches.
   * @param nTree the number of a Tree (from the collection)
   * @param nRate the number of a Distribution Rate (from the
   *        collection)
   * @param nFreq the number of a FrequencySet for the root (from
   *        the collection)
   * @param sharedParameterNames the vector of the names of parameters
   *        of the model that are shared among all the models.
   *
   */
  void addOnePerBranchSubstitutionProcess(size_t nProc, size_t nMod, size_t nTree, size_t nRate, const std::vector<std::string>& sharedParameterNames);
 
  size_t getNumberOfSubstitutionProcess() const { return mSubProcess_.size(); }
 
  SubstitutionProcessCollectionMember& substitutionProcess(size_t i)
  {
    return *mSubProcess_[i];
  }
 
  const SubstitutionProcessCollectionMember& substitutionProcess(size_t i) const
  {
    const auto it = mSubProcess_.find(i);
 
    return *(it->second);
  }
 
  std::shared_ptr<SubstitutionProcessCollectionMember> getSubstitutionProcess(size_t i)
  {
    return mSubProcess_[i];
  }
 
  std::shared_ptr<const SubstitutionProcessCollectionMember> getSubstitutionProcess(size_t i) const
  {
    const auto it = mSubProcess_.find(i);
 
    return it->second;
  }
 
  /*
   * @brief Get the branch lengths parameters.
   *
   * @return A ParameterList with all branch lengths.
   */
  ParameterList getBranchLengthParameters(bool independent) const
  {
    ParameterList pl = treeColl_.getParameters();
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getBranchLengthParameters(size_t nTree, bool independent) const
  {
    ParameterList pl = treeColl_.getParametersForObject(nTree);
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getSubstitutionProcessParameters(size_t nProc, bool independent) const
  {
    ParameterList pl;
 
    if (mSubProcess_.find(nProc) != mSubProcess_.end())
      pl = mSubProcess_.find(nProc)->second->getParameters();
    else
      return ParameterList();
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getSubstitutionProcessParameters() const;
 
  ParameterList getSubstitutionModelParameters(size_t nMod, bool independent) const
  {
    ParameterList pl = modelColl_.getParametersForObject(nMod);
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getSubstitutionModelParameters(bool independent) const
  {
    ParameterList pl = modelColl_.getParameters();
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
 
  ParameterList getNonDerivableParameters() const;
 
  ParameterList getRateDistributionParameters(size_t nRate, bool independent) const
  {
    ParameterList pl = distColl_.getParametersForObject(nRate);
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getRateDistributionParameters(bool independent) const
  {
    ParameterList pl = distColl_.getParameters();
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getRootFrequenciesParameters(size_t nFreq, bool independent) const
  {
    ParameterList pl = freqColl_.getParametersForObject(nFreq);
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
 
  ParameterList getRootFrequenciesParameters(bool independent) const
  {
    ParameterList pl = freqColl_.getParameters();
 
    if (independent)
      return pl.getCommonParametersWith(getIndependentParameters());
    else
      return pl;
  }
};
} // end of namespace bpp.
#endif // BPP_PHYL_LIKELIHOOD_SUBSTITUTIONPROCESSCOLLECTION_H