SubstitutionModelSet.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LEGACY_MODEL_SUBSTITUTIONMODELSET_H
#define BPP_PHYL_LEGACY_MODEL_SUBSTITUTIONMODELSET_H
 
#include <Bpp/Exceptions.h>
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Numeric/VectorTools.h>
 
#include "../../Model/AbstractSubstitutionModel.h"
#include "../../Model/FrequencySet/FrequencySet.h"
#include "../../Tree/Tree.h"
 
// From bpp-seq:
#include <Bpp/Seq/Alphabet/Alphabet.h>
#include <Bpp/Seq/Alphabet/NucleicAlphabet.h>
 
// From the STL:
#include <vector>
#include <map>
#include <algorithm>
#include <memory>
#include <typeinfo>
 
namespace bpp
{
class SubstitutionModelSet :
  public AbstractParameterAliasable
{
protected:
  std::shared_ptr<const Alphabet> alphabet_;
 
  size_t nbStates_;
 
  std::vector< std::shared_ptr<TransitionModelInterface>> modelSet_;
 
private:
  std::shared_ptr<FrequencySetInterface> rootFrequencies_;
 
  mutable std::map<int, size_t> nodeToModel_;
  mutable std::map<size_t, std::vector<int>> modelToNodes_;
 
  std::vector<ParameterList> modelParameters_;
 
  bool stationarity_;
 
public:
  SubstitutionModelSet(std::shared_ptr<const Alphabet> alpha) :
    AbstractParameterAliasable(""),
    alphabet_(alpha),
    nbStates_(0),
    modelSet_(),
    rootFrequencies_(),
    nodeToModel_(),
    modelToNodes_(),
    modelParameters_(),
    stationarity_(true)
  {}
 
  SubstitutionModelSet(
      std::shared_ptr<const Alphabet> alpha,
      std::shared_ptr<FrequencySetInterface> rootFreqs) :
    AbstractParameterAliasable(""),
    alphabet_(alpha),
    nbStates_(0),
    modelSet_(),
    rootFrequencies_(),
    nodeToModel_(),
    modelToNodes_(),
    modelParameters_(),
    stationarity_(true)
  {
    setRootFrequencies(rootFreqs);
  }
 
  void clear();
 
  bool isStationary() const
  {
    return stationarity_;
  }
 
  void setRootFrequencies(std::shared_ptr<FrequencySetInterface> rootFreqs);
 
  SubstitutionModelSet(const SubstitutionModelSet& set);
 
  SubstitutionModelSet& operator=(const SubstitutionModelSet& set);
 
  virtual ~SubstitutionModelSet() {}
 
  SubstitutionModelSet* clone() const { return new SubstitutionModelSet(*this); }
 
public:
  size_t getNumberOfStates() const
  {
    return nbStates_;
  }
 
  virtual void fireParameterChanged(const ParameterList& parameters);
 
  size_t getNumberOfModels() const { return modelSet_.size(); }
 
  bool hasMixedTransitionModel() const;
 
  std::shared_ptr<const TransitionModelInterface> getModel(size_t i) const
  {
    if (i >= modelSet_.size()) throw IndexOutOfBoundsException("SubstitutionModelSet::getModel(i).", i, 0, modelSet_.size() - 1);
    return modelSet_[i];
  }
 
  const TransitionModelInterface& model(size_t i) const
  {
    if (i >= modelSet_.size()) throw IndexOutOfBoundsException("SubstitutionModelSet::model(i).", i, 0, modelSet_.size() - 1);
    return *modelSet_[i];
  }
 
  std::shared_ptr<TransitionModelInterface> getModel(size_t i)
  {
    if (i >= modelSet_.size()) throw IndexOutOfBoundsException("SubstitutionModelSet::getModel(i).", i, 0, modelSet_.size() - 1);
    return modelSet_[i];
  }
 
  TransitionModelInterface& model(size_t i)
  {
    if (i >= modelSet_.size()) throw IndexOutOfBoundsException("SubstitutionModelSet::model(i).", i, 0, modelSet_.size() - 1);
    return *modelSet_[i];
  }
  std::shared_ptr<const SubstitutionModelInterface> getSubstitutionModel(size_t i) const
  {
    auto m = std::dynamic_pointer_cast<const SubstitutionModelInterface>(getModel(i));
    if (m) return m;
    else
      throw Exception("SubstitutionModelSet::getSubstitutionModel : " + model(i).getName() + " is  not a substitution model.");
  }
 
 
  std::shared_ptr<SubstitutionModelInterface> getSubstitutionModel(size_t i)
  {
    auto m = std::dynamic_pointer_cast<SubstitutionModelInterface>(getModel(i));
    if (m) return m;
    else
      throw Exception("SubstitutionModelSet::getSubstitutionModel : " + model(i).getName() + " is  not a substitution model.");
  }
 
  const SubstitutionModelInterface& substitutionModel(size_t i) const
  {
    try
    {
      auto& m = dynamic_cast<const SubstitutionModelInterface&>(model(i));
      return m;
    }
    catch (Exception& ex)
    {
      throw Exception("SubstitutionModelSet::substitutionModel : " + model(i).getName() + " is not a substitution model.");
    }
  }
 
  bool hasOnlySubstitutionModels() const
  {
    for (const auto& mod : modelSet_)
    {
      auto model = std::dynamic_pointer_cast<const SubstitutionModelInterface>(mod);
      if (!model) return false;
    }
    return true;
  }
 
  size_t getModelIndexForNode(int nodeId) const
  {
    auto i = nodeToModel_.find(nodeId);
    if (i == nodeToModel_.end())
      throw Exception("SubstitutionModelSet::getModelIndexForNode(). No model associated to node with id " + TextTools::toString(nodeId));
    return i->second;
  }
 
  std::shared_ptr<const TransitionModelInterface> getModelForNode(int nodeId) const
  {
    auto i = nodeToModel_.find(nodeId);
    if (i == nodeToModel_.end())
      throw Exception("SubstitutionModelSet::getModelForNode(). No model associated to node with id " + TextTools::toString(nodeId));
    return modelSet_[i->second];
  }
  std::shared_ptr<TransitionModelInterface> getModelForNode(int nodeId)
  {
    auto i = nodeToModel_.find(nodeId);
    if (i == nodeToModel_.end())
      throw Exception("SubstitutionModelSet::getModelForNode(). No model associated to node with id " + TextTools::toString(nodeId));
    return modelSet_[i->second];
  }
 
  std::shared_ptr<const SubstitutionModelInterface> getSubstitutionModelForNode(int nodeId) const
  {
    return std::dynamic_pointer_cast<const SubstitutionModelInterface>(getModelForNode(nodeId));
  }
 
  std::shared_ptr<SubstitutionModelInterface> getSubstitutionModelForNode(int nodeId)
  {
    return std::dynamic_pointer_cast<SubstitutionModelInterface>(getModelForNode(nodeId));
  }
 
  const std::vector<int>& getNodesWithModel(size_t i) const
  {
    if (i >= modelSet_.size()) throw IndexOutOfBoundsException("SubstitutionModelSet::getNodesWithModel().", i, 0, modelSet_.size());
    return modelToNodes_[i];
  }
 
  std::vector<int> getNodesWithParameter(const std::string& name) const;
 
  void addModel(std::shared_ptr<TransitionModelInterface> model, const std::vector<int>& nodesId); // , const std::vector<std::string>& newParams);
 
  void setNodeToModel(size_t modelIndex, int nodeId); // Keren: added on my own to allow alternation of nodes assignments to existing nodes
 
  void resetModelToNodeIds();  // Keren: added on my own to allow alternation of nodes assignments to existing nodes
 
  void replaceModel(size_t modelIndex, std::shared_ptr<TransitionModelInterface> model);
 
  void listModelNames(std::ostream& out = std::cout) const;
 
  const std::shared_ptr<FrequencySetInterface> getRootFrequencySet() const { return rootFrequencies_; }
 
  const FrequencySetInterface& rootFrequencySet() const
  {
    if (rootFrequencies_)
    {
      return *rootFrequencies_;
    }
    else
    {
      throw NullPointerException("SubstitutionModelSet::rootFrequencySet. No associated root frequencies.");
    }
  }
 
  std::vector<double> getRootFrequencies() const
  {
    if (stationarity_)
      return modelSet_[0]->getFrequencies();
    else
      return rootFrequencies_->getFrequencies();
  }
 
  ParameterList getRootFrequenciesParameters() const
  {
    if (stationarity_)
      return ParameterList();
    else
      return rootFrequencies_->getParameters();
  }
 
  ParameterList getNodeParameters() const
  {
    ParameterList pl;
    for (size_t i = stationarity_ ? 0 : rootFrequencies_->getNumberOfParameters();
        i < getNumberOfParameters(); i++)
    {
      pl.addParameter(getParameter_(i));
    }
    return pl;
  }
 
  ParameterList getModelParameters(size_t modelIndex) const;
 
  std::shared_ptr<const Alphabet> getAlphabet() const { return alphabet_; }
 
  const Alphabet& alphabet() const { return *alphabet_; }
 
  const std::vector<int>& getAlphabetStates() const
  {
    return model(0).getAlphabetStates();
  }
 
  const StateMapInterface& stateMap() const
  {
    return model(0).stateMap();
  }
 
  std::shared_ptr<const StateMapInterface> getStateMap() const
  {
    return model(0).getStateMap();
  }
 
  std::vector<size_t> getModelStates(int code) const
  {
    return model(0).getModelStates(code);
  }
 
  std::vector<size_t> getModelStates(const std::string& code) const
  {
    return model(0).getModelStates(code);
  }
 
  int getAlphabetStateAsInt(size_t index) const
  {
    return model(0).getAlphabetStateAsInt(index);
  }
 
  std::string getAlphabetStateAsChar(size_t index) const
  {
    return model(0).getAlphabetStateAsChar(index);
  }
 
  bool isFullySetUpFor(const Tree& tree, bool throwEx = true) const
  {
    return checkOrphanModels(throwEx)
           //           && checkOrphanParameters(throwEx)
           && checkOrphanNodes(tree, throwEx)
           && checkUnknownNodes(tree, throwEx);
  }
 
protected:
  void updateRootFrequencies()
  {
    if (!stationarity_)
      rootFrequencies_->matchParametersValues(getParameters());
  }
 
  bool checkOrphanModels(bool throwEx) const;
 
  bool checkOrphanNodes(const Tree& tree, bool throwEx) const;
 
  bool checkUnknownNodes(const Tree& tree, bool throwEx) const;
};
} // end of namespace bpp.
#endif // BPP_PHYL_LEGACY_MODEL_SUBSTITUTIONMODELSET_H