NonHomogeneousSequenceSimulator.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef _NONHOMOGENEOUSSEQUENCESIMULATOR_H_
#define _NONHOMOGENEOUSSEQUENCESIMULATOR_H_
 
#include "../../Simulation/DetailedSiteSimulator.h"
#include "../../Simulation/SequenceSimulator.h"
#include "../../Tree/TreeTemplate.h"
#include "../../Tree/NodeTemplate.h"
#include "../../Model/SubstitutionModel.h"
#include "../../Likelihood/ParametrizablePhyloTree.h"
 
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Numeric/Prob/DiscreteDistribution.h>
 
// From bpp-seq:
#include <Bpp/Seq/Alphabet/Alphabet.h>
#include <Bpp/Seq/Site.h>
#include <Bpp/Seq/Container/SiteContainer.h>
 
// From the STL:
#include <map>
#include <vector>
 
#include "../Model/SubstitutionModelSet.h"
 
namespace bpp
{
class SimData
{
public:
  size_t state;
  std::vector<size_t> states;
  VVVdouble cumpxy;
  std::shared_ptr<const TransitionModelInterface> model;
 
public:
  SimData() : state(), states(), cumpxy(), model(0) {}
  SimData(const SimData& sd) : state(sd.state), states(sd.states), cumpxy(), model(sd.model) {}
  SimData& operator=(const SimData& sd)
  {
    state  = sd.state;
    states = sd.states;
    cumpxy = sd.cumpxy;
    model  = sd.model;
    return *this;
  }
};
 
typedef NodeTemplate<SimData> SNode;
 
class NonHomogeneousSequenceSimulator :
  public DetailedSiteSimulatorInterface,
  public virtual SequenceSimulatorInterface
{
private:
  std::shared_ptr<const SubstitutionModelSet> modelSet_;
  mutable std::shared_ptr<const Alphabet> alphabet_;
  std::vector<int>            supportedStates_;
  std::shared_ptr<const DiscreteDistributionInterface> rate_;
  std::shared_ptr<const Tree> templateTree_;
  mutable TreeTemplate<SNode> tree_;
  mutable std::shared_ptr<const ParametrizablePhyloTree> phyloTree_;
 
  std::vector<SNode*> leaves_;
 
  std::vector<std::string> seqNames_;
 
  size_t nbNodes_;
  size_t nbClasses_;
  size_t nbStates_;
 
  bool continuousRates_;
 
  // Should we output internal sequences as well?
  bool outputInternalSequences_;
 
public:
  NonHomogeneousSequenceSimulator(
      std::shared_ptr<const SubstitutionModelSet> modelSet,
      std::shared_ptr<const DiscreteDistributionInterface> rate,
      std::shared_ptr<const Tree> tree
      );
 
  NonHomogeneousSequenceSimulator(
      std::shared_ptr<const TransitionModelInterface> model,
      std::shared_ptr<const DiscreteDistributionInterface> rate,
      std::shared_ptr<const Tree> tree
      );
 
  virtual ~NonHomogeneousSequenceSimulator() {}
 
  NonHomogeneousSequenceSimulator(const NonHomogeneousSequenceSimulator& nhss) = default;
 
  NonHomogeneousSequenceSimulator& operator=(const NonHomogeneousSequenceSimulator& nhss) = default;
 
  NonHomogeneousSequenceSimulator* clone() const override
  {
    return new NonHomogeneousSequenceSimulator(*this);
  }
 
private:
  void init();
 
public:
  std::unique_ptr<Site> simulateSite() const override;
 
  std::unique_ptr<Site> simulateSite(size_t ancestralStateIndex) const override;
 
  std::unique_ptr<Site> simulateSite(size_t ancestralStateIndex, double rate) const override;
 
  std::unique_ptr<Site> simulateSite(double rate) const override;
 
  std::vector<std::string> getSequenceNames() const override { return seqNames_; }
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite() const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(size_t ancestralStateIndex) const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(size_t ancestralStateIndex, double rate) const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(double rate) const override;
  std::unique_ptr<SiteContainerInterface> simulate(size_t numberOfSites) const override;
 
  const SiteSimulatorInterface& siteSimulator(size_t pos) const override
  {
    throw Exception("NonHomogeneousSequenceSimulator::siteSimulator: not implemented.");
  }
  std::shared_ptr<const Alphabet> getAlphabet() const override { return alphabet_; }
 
  const Alphabet& alphabet() const override { return *alphabet_; }
  virtual std::unique_ptr<Site> simulateSite(size_t ancestralStateIndex, size_t rateClass) const;
  virtual std::unique_ptr<SiteSimulationResult> dSimulateSite(size_t ancestralStateIndex, size_t rateClass) const;
  std::shared_ptr<const SubstitutionModelSet> getSubstitutionModelSet() const { return modelSet_; }
 
  const SubstitutionModelSet& substitutionModelSet() const { return *modelSet_; }
 
  std::shared_ptr<const DiscreteDistributionInterface> getRateDistribution() const { return rate_; }
 
  const DiscreteDistributionInterface& rateDistribution() const { return *rate_; }
 
  std::shared_ptr<const Tree> getTree() const { return templateTree_; }
 
  const Tree& tree() const { return *templateTree_; }
 
  void enableContinuousRates(bool yn) { continuousRates_ = yn; }
 
  void outputInternalSequences(bool yn) override;
  void outputInternalSites(bool yn) override
  {
    throw Exception("NonHomogeneousSequenceSimulator::outputInternalSites: not implemented.");
  }
 
protected:
  size_t evolve(const SNode* node, size_t initialStateIndex, size_t rateClass) const;
 
  size_t evolve(const SNode* node, size_t initialStateIndex, double rate) const;
 
  void multipleEvolve(
      const SNode* node,
      const std::vector<size_t>& initialStateIndices,
      const std::vector<size_t>& rateClasses,
      std::vector<size_t>& finalStates) const;
 
  std::unique_ptr<SiteContainerInterface> multipleEvolve(
      const std::vector<size_t>& initialStates,
      const std::vector<size_t>& rateClasses) const;
 
  void dEvolve(size_t initialState, double rate, RASiteSimulationResult& rassr) const;
 
  void evolveInternal(SNode* node, size_t rateClass) const;
  void evolveInternal(SNode* node, double rate) const;
  void multipleEvolveInternal(SNode* node, const std::vector<size_t>& rateClasses) const;
 
  void dEvolveInternal(SNode* node, double rate, RASiteSimulationResult& rassr) const;
};
} // end of namespace bpp.
 
#endif // _NONHOMOGENEOUSSEQUENCESIMULATOR_H_