BackwardLikelihoodTree.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_LIKELIHOOD_DATAFLOW_BACKWARDLIKELIHOODTREE_H
#define BPP_PHYL_LIKELIHOOD_DATAFLOW_BACKWARDLIKELIHOODTREE_H
 
 
#include "Bpp/Phyl/Likelihood/DataFlow/ForwardLikelihoodTree.h"
#include "Bpp/Phyl/Likelihood/DataFlow/Model.h"
#include "Bpp/Phyl/Likelihood/DataFlow/ProcessTree.h"
 
namespace bpp
{
// using RowLik = Eigen::Matrix<double, 1, Eigen::Dynamic>;
// using MatrixLik = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
 
 
using SpeciationBackward = CWiseMul<MatrixLik, ReductionOf<MatrixLik>>;
 
using MixtureBackward = CWiseAdd<MatrixLik, ReductionOf<MatrixLik>>;
 
using BackwardTransition =
    MatrixProduct<MatrixLik, Transposed<Eigen::MatrixXd>, MatrixLik>;
 
using BackwardProportion = CWiseMul<MatrixLik, std::tuple<double, MatrixLik>>;
 
// Upper Likelihood in nodes
using ConditionalLikelihood = Value<MatrixLik>;
using ConditionalLikelihoodRef = ValueRef<MatrixLik>;
 
// Upper Likelihood at top of edges
using BackwardLikelihoodAbove = Value<MatrixLik>;
using BackwardLikelihoodAboveRef = ValueRef<MatrixLik>;
 
//    using FullLikelihood = Value<Patterned<MatrixLik>>;
 
/* All the computations are set in a DataFlow context */
 
class BackwardLikelihoodTree : public AssociationDAGlobalGraphObserver<ConditionalLikelihood, BackwardLikelihoodAbove>
{
  using DAClass = AssociationDAGlobalGraphObserver<ConditionalLikelihood, BackwardLikelihoodAbove>;
 
private:
  Context& context_;
  Eigen::Index nbState_;
  Eigen::Index nbSite_;
  std::shared_ptr<ForwardLikelihoodTree> forwardTree_;
  std::shared_ptr<ProcessTree> processTree_;
  ValueRef<Eigen::RowVectorXd> rFreqs_;
  MatrixDimension likelihoodMatrixDim_;
  const StateMapInterface& statemap_;
 
public:
  BackwardLikelihoodTree(Context& c,
      std::shared_ptr<ForwardLikelihoodTree> forwardTree,
      std::shared_ptr<ProcessTree> tree,
      ValueRef<Eigen::RowVectorXd> rFreqs,
      const StateMapInterface& statemap,
      Eigen::Index nbSite) :
    DAClass(forwardTree->getGraph()),
    context_(c), nbState_(Eigen::Index(statemap.getNumberOfModelStates())), nbSite_(nbSite), forwardTree_(forwardTree), processTree_(tree), rFreqs_(rFreqs), likelihoodMatrixDim_(conditionalLikelihoodDimension (nbState_, nbSite_)), statemap_(statemap)
  {}
 
  ConditionalLikelihoodRef setRootFrequencies(const ValueRef<Eigen::RowVectorXd> rootFreqs)
  {
    auto rootFreqsEf = Convert<ExtendedFloatRowVectorXd, Eigen::RowVectorXd>::create(context_, {rootFreqs}, RowVectorDimension (Eigen::Index (nbState_)));
    auto r2 = CWiseFill<MatrixLik, RowLik>::create(context_, {rootFreqsEf}, likelihoodMatrixDim_);
 
    associateNode(r2, forwardTree_->getNodeGraphid(forwardTree_->getRoot()));
    setNodeIndex(r2, forwardTree_->getRootIndex());
    return r2;
  }
 
private:
  /*
   * @brief Compute joined likelihood BEFORE reading edge on the
   * backward process (ie at top of the edge).
   *
   */
 
  BackwardLikelihoodAboveRef makeBackwardLikelihoodAtEdge (PhyloTree::EdgeIndex index);
 
  /*
   * @brief Compute joined likelihood at node on the backward
   * process (ie with the bottom of the edges getting to the node
   * from upward).
   *
   */
 
  ConditionalLikelihoodRef makeBackwardLikelihoodAtNode (PhyloTree::NodeIndex index);
 
  /*
   * @brief the LikehoodArrays
   *
   * Beware: nodeIds are in the DAG, not the ids of the PhyloTree.
   *
   * Set in private to avoid bad usage, access through
   * LikelihoodCalculationSingleProcess.
   */
  ConditionalLikelihoodRef getBackwardLikelihoodArray(PhyloTree::NodeIndex nodeId)
  {
    if (!hasNode(nodeId))
      makeBackwardLikelihoodAtNode(nodeId);
 
    return getNode(nodeId);
  }
 
  friend class LikelihoodCalculationSingleProcess;
};
} // end of namespace bpp.
#endif // BPP_PHYL_LIKELIHOOD_DATAFLOW_BACKWARDLIKELIHOODTREE_H