NNIHomogeneousTreeLikelihood.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LEGACY_LIKELIHOOD_NNIHOMOGENEOUSTREELIKELIHOOD_H
#define BPP_PHYL_LEGACY_LIKELIHOOD_NNIHOMOGENEOUSTREELIKELIHOOD_H
 
#include <Bpp/Numeric/Function/BrentOneDimension.h>
#include <Bpp/Numeric/Parametrizable.h>
#include <Bpp/Numeric/Prob/DiscreteDistribution.h>
#include <Bpp/Numeric/VectorTools.h>
 
#include "../../Tree/NNISearchable.h"
#include "DRHomogeneousTreeLikelihood.h"
 
namespace bpp
{
class BranchLikelihood :
  public FunctionInterface,
  public AbstractParametrizable
{
protected:
  const VVVdouble* array1_, * array2_;
  std::shared_ptr<const TransitionModelInterface> model_;
  std::shared_ptr<const DiscreteDistributionInterface> rDist_;
  size_t nbStates_, nbClasses_;
  VVVdouble pxy_;
  double lnL_;
  std::vector<unsigned int> weights_;
 
public:
  BranchLikelihood(const std::vector<unsigned int>& weights) :
    AbstractParametrizable(""),
    array1_(0),
    array2_(0),
    model_(0),
    rDist_(0),
    nbStates_(0),
    nbClasses_(0),
    pxy_(),
    lnL_(std::log(0.)),
    weights_(weights)
  {
    addParameter_(new Parameter("BrLen", 1, 0));
  }
 
  BranchLikelihood(const BranchLikelihood& bl) :
    AbstractParametrizable(bl),
    array1_(bl.array1_),
    array2_(bl.array2_),
    model_(bl.model_),
    rDist_(bl.rDist_),
    nbStates_(bl.nbStates_),
    nbClasses_(bl.nbClasses_),
    pxy_(bl.pxy_),
    lnL_(bl.lnL_),
    weights_(bl.weights_)
  {}
 
  BranchLikelihood& operator=(const BranchLikelihood& bl)
  {
    AbstractParametrizable::operator=(bl);
    array1_ = bl.array1_;
    array2_ = bl.array2_;
    model_ = bl.model_;
    rDist_ = bl.rDist_;
    nbStates_ = bl.nbStates_;
    nbClasses_ = bl.nbClasses_;
    pxy_ = bl.pxy_;
    lnL_ = bl.lnL_;
    weights_ = bl.weights_;
    return *this;
  }
 
  virtual ~BranchLikelihood() {}
 
  BranchLikelihood* clone() const { return new BranchLikelihood(*this); }
 
public:
  void initModel(
      std::shared_ptr<const TransitionModelInterface> model,
      std::shared_ptr<const DiscreteDistributionInterface> rDist);
 
  void initLikelihoods(const VVVdouble* array1, const VVVdouble* array2)
  {
    array1_ = array1;
    array2_ = array2;
  }
 
  void resetLikelihoods()
  {
    array1_ = 0;
    array2_ = 0;
  }
 
  void setParameters(const ParameterList& parameters)
  {
    setParametersValues(parameters);
  }
 
  double getValue() const { return lnL_; }
 
  void fireParameterChanged(const ParameterList& parameters)
  {
    computeAllTransitionProbabilities();
    computeLogLikelihood();
  }
 
protected:
  void computeAllTransitionProbabilities();
  void computeLogLikelihood();
};
 
 
class NNIHomogeneousTreeLikelihood :
  public DRHomogeneousTreeLikelihood,
  public virtual NNISearchable
{
protected:
  std::shared_ptr<BranchLikelihood> brLikFunction_;
 
  std::unique_ptr<BrentOneDimension> brentOptimizer_;
 
  mutable std::map<int, double> brLenNNIValues_;
 
  ParameterList brLenNNIParams_;
 
public:
  NNIHomogeneousTreeLikelihood(
      const Tree& tree,
      std::shared_ptr<TransitionModelInterface> model,
      std::shared_ptr<DiscreteDistributionInterface> rDist,
      bool checkRooted = true,
      bool verbose = true);
 
  NNIHomogeneousTreeLikelihood(
      const Tree& tree,
      const AlignmentDataInterface& data,
      std::shared_ptr<TransitionModelInterface> model,
      std::shared_ptr<DiscreteDistributionInterface> rDist,
      bool checkRooted = true,
      bool verbose = true);
 
  NNIHomogeneousTreeLikelihood(const NNIHomogeneousTreeLikelihood& lik);
 
  NNIHomogeneousTreeLikelihood& operator=(const NNIHomogeneousTreeLikelihood& lik);
 
  virtual ~NNIHomogeneousTreeLikelihood();
 
  NNIHomogeneousTreeLikelihood* clone() const override { return new NNIHomogeneousTreeLikelihood(*this); }
 
public:
  void setData(const AlignmentDataInterface& sites) override
  {
    DRHomogeneousTreeLikelihood::setData(sites);
    brLikFunction_ = std::make_unique<BranchLikelihood>(likelihoodData().getWeights());
  }
 
  const Tree& topology() const override { return tree(); }
 
  double getTopologyValue() const override { return getValue(); }
 
  double testNNI(int nodeId) const override;
 
  void doNNI(int nodeId) override;
 
  void topologyChangeTested(const TopologyChangeEvent& event) override
  {
    likelihoodData().reInit();
    // if(brLenNNIParams_.size() > 0)
    fireParameterChanged(brLenNNIParams_);
    brLenNNIParams_.reset();
  }
 
  void topologyChangeSuccessful(const TopologyChangeEvent& event) override
  {
    brLenNNIValues_.clear();
  }
};
} // end of namespace bpp.
#endif // BPP_PHYL_LEGACY_LIKELIHOOD_NNIHOMOGENEOUSTREELIKELIHOOD_H