MutationProcess.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_SIMULATION_MUTATIONPROCESS_H
#define BPP_PHYL_SIMULATION_MUTATIONPROCESS_H
 
#include <Bpp/Numeric/VectorTools.h>
 
#include "../Mapping/SubstitutionRegister.h"
#include "../Model/SubstitutionModel.h"
 
namespace bpp
{
class MutationPath
{
private:
  std::shared_ptr<const Alphabet> alphabet_;
 
  std::vector<size_t> states_;
 
  std::vector<double> times_;
 
  size_t initialState_;
 
  double totalTime_;
 
public:
  MutationPath(
      std::shared_ptr<const Alphabet> alphabet,
      size_t initialState,
      double time) :
    alphabet_(alphabet), states_(), times_(), initialState_(initialState), totalTime_(time) {}
 
  MutationPath(const MutationPath& path) :
    alphabet_(path.alphabet_), states_(path.states_), times_(path.times_), initialState_(path.initialState_), totalTime_(path.totalTime_) {}
 
  MutationPath& operator=(const MutationPath& path)
  {
    alphabet_     = path.alphabet_;
    states_       = path.states_;
    times_        = path.times_;
    initialState_ = path.initialState_;
    totalTime_    = path.totalTime_;
    return *this;
  }
 
  virtual ~MutationPath() {}
 
public:
  std::shared_ptr<const Alphabet> getAlphabet() const { return alphabet_; }
 
  const Alphabet& alphabet() const { return *alphabet_; }
 
  void addEvent(size_t state, double time)
  {
    states_.push_back(state);
    times_.push_back(time);
  }
 
  /*
   * @brief Remove all mutations
   *
   */
  void clear()
  {
    states_.clear();
    times_.clear();
  }
 
  size_t getInitialState() const { return initialState_; }
 
  double getTotalTime() const { return totalTime_; }
 
  size_t getNumberOfEvents() const { return states_.size(); }
 
  template<class Scalar>
  void getEventCounts(Matrix<Scalar>& counts) const
  {
    if (counts.getNumberOfRows()    != alphabet_->getSize()
        || counts.getNumberOfColumns() != alphabet_->getSize())
      throw Exception("MutationPath::getEventCounts. Incorrect input matrix, does not match alphabet size.");
    size_t currentState = initialState_;
    for (size_t i = 0; i < states_.size(); ++i)
    {
      size_t newState = states_[i];
      counts(currentState, newState)++;
      currentState = newState;
    }
  }
 
  template<class Scalar>
  void getEventCounts(std::vector<Scalar>& counts, const SubstitutionRegisterInterface& reg) const
  {
    if (counts.size() != reg.getNumberOfSubstitutionTypes())
      throw Exception("MutationPath::getEventCounts. Incorrect input vector, does not match alphabet size.");
    size_t currentState = initialState_;
    for (size_t i = 0; i < states_.size(); ++i)
    {
      size_t newState = states_[i];
      size_t type = reg.getType(currentState, newState);
      if (type > 0) counts[type - 1]++;
      currentState = newState;
    }
  }
 
  size_t getFinalState() const
  {
    if (states_.size() == 0) return initialState_;
    else return states_[states_.size() - 1];
  }
};
 
class MutationProcess
{
public:
  MutationProcess() {}
  virtual ~MutationProcess() {}
 
public:
  virtual size_t mutate(size_t state) const = 0;
 
  virtual size_t mutate(size_t state, unsigned int n) const = 0;
 
  virtual double getTimeBeforeNextMutationEvent(size_t state) const = 0;
 
  virtual size_t evolve(size_t initialState, double time) const = 0;
 
  virtual MutationPath detailedEvolve(size_t initialState, double time) const = 0;
 
  virtual MutationPath detailedEvolve(size_t initialState, size_t finalState, double time) const = 0;
 
  virtual std::shared_ptr<const SubstitutionModelInterface> getSubstitutionModel() const = 0;
};
 
class AbstractMutationProcess :
  public virtual MutationProcess
{
protected:
  std::shared_ptr<const SubstitutionModelInterface> model_;
 
  size_t size_;
 
  VVdouble repartition_;
 
public:
  AbstractMutationProcess(std::shared_ptr<const SubstitutionModelInterface> model) :
    model_(model), size_(), repartition_()
  {}
 
public:
  size_t mutate(size_t state) const;
  size_t mutate(size_t state, unsigned int n) const;
  double getTimeBeforeNextMutationEvent(size_t state) const;
  size_t evolve(size_t initialState, double time) const;
  MutationPath detailedEvolve(size_t initialState, double time) const;
  MutationPath detailedEvolve(size_t initialState, size_t finalState, double time) const;
  std::shared_ptr<const SubstitutionModelInterface> getSubstitutionModel() const { return model_; }
};
 
class SimpleMutationProcess :
  public AbstractMutationProcess
{
public:
  // Constructor and destructor:
 
  SimpleMutationProcess(std::shared_ptr<const SubstitutionModelInterface> model);
 
  virtual ~SimpleMutationProcess();
 
  size_t evolve(size_t initialState, double time) const;
};
 
class SelfMutationProcess : public AbstractMutationProcess
{
public:
  SelfMutationProcess(size_t alphabetSize);
 
  virtual ~SelfMutationProcess();
};
} // end of namespace bpp.
#endif // BPP_PHYL_SIMULATION_MUTATIONPROCESS_H