SimpleSubstitutionProcessSiteSimulator.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_SIMULATION_SIMPLESUBSTITUTIONPROCESSSITESIMULATOR_H
#define BPP_PHYL_SIMULATION_SIMPLESUBSTITUTIONPROCESSSITESIMULATOR_H
 
#include <Bpp/Numeric/Prob/DiscreteDistribution.h>
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Phyl/Likelihood/ProcessComputationTree.h>
 
#include "../Likelihood/ParametrizablePhyloTree.h"
#include "../Model/SubstitutionModel.h"
#include "DetailedSiteSimulator.h"
 
// From SeqLib:
#include <Bpp/Seq/Site.h>
 
// From the STL:
#include <map>
#include <vector>
 
#include "../Likelihood/SubstitutionProcess.h"
 
namespace bpp
{
class SimProcessNode :
  public ProcessComputationNode
{
private:
  // states during simulation
  size_t state_;
 
  // probabilities to choose, in case of mixture node or root, for all rates
  VVdouble cumProb_;
 
  // Sons in case of mixture node
  std::vector<std::shared_ptr<SimProcessNode>> sons_;
 
public:
  SimProcessNode(const ProcessComputationNode& pcn) :
    ProcessComputationNode(pcn), state_() {}
 
  friend class SimpleSubstitutionProcessSiteSimulator;
  friend class GivenDataSubstitutionProcessSiteSimulator;
};
 
class SimProcessEdge :
  public ProcessComputationEdge
{
private:
  // Cumulative pxy for all rates
  VVVdouble cumpxy_;
 
public:
  SimProcessEdge(const ProcessComputationEdge& pce) :
    ProcessComputationEdge(pce), cumpxy_() {}
 
  friend class SimpleSubstitutionProcessSiteSimulator;
  friend class GivenDataSubstitutionProcessSiteSimulator;
};
 
typedef AssociationTreeGlobalGraphObserver<SimProcessNode, SimProcessEdge>  SPTree;
 
 
class SimpleSubstitutionProcessSiteSimulator :
  public virtual DetailedSiteSimulatorInterface
{
protected:
  std::shared_ptr<const SubstitutionProcessInterface> process_;
  std::shared_ptr<const ParametrizablePhyloTree> phyloTree_;
 
 
  SPTree tree_;
 
  Vdouble qRates_;
 
  /*
   * @brief cumsum probas of the root frequencies, one per rate class
   * (this "per class" is useful for GivenDataSubstitutionProcessSiteSimulator)
   */
  VVdouble qRoots_;
 
  std::vector<size_t> seqIndexes_;
 
  std::vector<std::string> seqNames_;
 
  mutable std::map<size_t, std::shared_ptr<SimProcessNode>> speciesNodes_;
 
  size_t nbNodes_;
  size_t nbClasses_;
  size_t nbStates_;
 
  bool continuousRates_;
 
  // Should we output internal sequences as well?
  bool outputInternalSites_;
 
public:
  SimpleSubstitutionProcessSiteSimulator(
      std::shared_ptr<const SubstitutionProcessInterface> process);
 
  virtual ~SimpleSubstitutionProcessSiteSimulator() {}
 
  SimpleSubstitutionProcessSiteSimulator(const SimpleSubstitutionProcessSiteSimulator& nhss) :
    process_        (nhss.process_),
    phyloTree_      (nhss.phyloTree_),
    tree_           (nhss.tree_),
    qRates_         (nhss.qRates_),
    qRoots_         (nhss.qRoots_),
    seqIndexes_     (nhss.seqIndexes_),
    seqNames_       (nhss.seqNames_),
    speciesNodes_  (nhss.speciesNodes_),
    nbNodes_        (nhss.nbNodes_),
    nbClasses_      (nhss.nbClasses_),
    nbStates_       (nhss.nbStates_),
    continuousRates_(nhss.continuousRates_),
    outputInternalSites_(nhss.outputInternalSites_)
  {}
 
  SimpleSubstitutionProcessSiteSimulator& operator=(const SimpleSubstitutionProcessSiteSimulator& nhss)
  {
    process_        = nhss.process_;
    phyloTree_       = nhss.phyloTree_;
    tree_            = nhss.tree_;
    qRates_          = nhss.qRates_;
    qRoots_          = nhss.qRoots_;
    seqIndexes_      = nhss.seqIndexes_;
    seqNames_        = nhss.seqNames_;
    speciesNodes_   = nhss.speciesNodes_;
    nbNodes_         = nhss.nbNodes_;
    nbClasses_       = nhss.nbClasses_;
    nbStates_        = nhss.nbStates_;
    continuousRates_ = nhss.continuousRates_;
    outputInternalSites_ = nhss.outputInternalSites_;
 
    return *this;
  }
 
  SimpleSubstitutionProcessSiteSimulator* clone() const override
  {
    return new SimpleSubstitutionProcessSiteSimulator(*this);
  }
 
private:
  virtual void init();
 
public:
  std::unique_ptr<Site> simulateSite() const override;
 
  std::unique_ptr<Site> simulateSite(size_t rateClass) const override;
 
  std::unique_ptr<Site> simulateSite(double rate) const override;
 
  std::unique_ptr<Site> simulateSite(size_t ancestralStateIndex, double rate) const override;
 
  std::vector<std::string> getSequenceNames() const override
  {
    return seqNames_;
  }
 
  std::shared_ptr<const Alphabet> getAlphabet() const override { return process_->stateMap().getAlphabet(); }
 
  const Alphabet& alphabet() const override { return process_->stateMap().alphabet(); }
  std::unique_ptr<SiteSimulationResult> dSimulateSite() const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(size_t rateClass) const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(double rate) const override;
 
  std::unique_ptr<SiteSimulationResult> dSimulateSite(size_t ancestralStateIndex, double rate) const override;
 
  std::shared_ptr<const SubstitutionProcessInterface> getSubstitutionProcess() const
  {
    return process_;
  }
 
  std::shared_ptr<const ParametrizablePhyloTree> getTree() const
  {
    return phyloTree_;
  }
 
  void enableContinuousRates(bool yn) { continuousRates_ = yn; }
 
  void outputInternalSites(bool yn) override;
 
protected:
  void evolveInternal(
      std::shared_ptr<SimProcessNode> node,
      size_t rateClass, SiteSimulationResult* ssr = nullptr) const;
 
  void evolveInternal(
      std::shared_ptr<SimProcessNode> node,
      double rate,
      SiteSimulationResult* ssr = nullptr) const;
 
};
} // end of namespace bpp.
#endif // BPP_PHYL_SIMULATION_SIMPLESUBSTITUTIONPROCESSSITESIMULATOR_H