DRASRTreeLikelihoodData.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#include "../../PatternTools.h"
#include "DRASRTreeLikelihoodData.h"
 
// From SeqLib:
#include <Bpp/Seq/SiteTools.h>
#include <Bpp/Seq/Container/AlignedSequenceContainer.h>
#include <Bpp/Seq/Container/SequenceContainerTools.h>
#include <Bpp/Seq/Container/VectorSiteContainer.h>
 
using namespace bpp;
using namespace std;
 
/******************************************************************************/
 
void DRASRTreeLikelihoodData::initLikelihoods(
    const AlignmentDataInterface& sites,
    const TransitionModelInterface& model)
{
  if (sites.getNumberOfSequences() == 1)
    throw Exception("Error, only 1 sequence!");
  if (sites.getNumberOfSequences() == 0)
    throw Exception("Error, no sequence!");
  if (sites.getAlphabet()->getAlphabetType()
      != model.getAlphabet()->getAlphabetType())
    throw AlphabetMismatchException("DRASDRTreeLikelihoodData::initLikelihoods. Data and model must have the same alphabet type.",
          sites.getAlphabet(),
          model.getAlphabet());
  alphabet_ = sites.getAlphabet();
  nbStates_ = model.getNumberOfStates();
  nbSites_  = sites.getNumberOfSites();
  shared_ptr<SitePatterns> patterns;
 
  if (usePatterns_)
  {
    patterns = initLikelihoodsWithPatterns(tree_->getRootNode(), sites, model);
    shrunkData_       = patterns->getSites();
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_.resize((size_t)patterns->getIndices().size());
    SitePatterns::IndicesType::Map(&rootPatternLinks_[0], patterns->getIndices().size()) = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();
  }
  else
  {
    patterns = std::make_unique<SitePatterns>(sites);
    shrunkData_       = patterns->getSites();
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_.resize(size_t(patterns->getIndices().size()));
    SitePatterns::IndicesType::Map(&rootPatternLinks_[0], patterns->getIndices().size()) = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();
    initLikelihoods(tree_->getRootNode(), *shrunkData_, model);
  }
}
 
/******************************************************************************/
 
void DRASRTreeLikelihoodData::initLikelihoods(
    const Node* node,
    const AlignmentDataInterface& sequences,
    const TransitionModelInterface& model)
{
  // Initialize likelihood vector:
  DRASRTreeLikelihoodNodeData* nodeData = &nodeData_[node->getId()];
  nodeData->setNode(node);
  VVVdouble* _likelihoods_node = &nodeData->getLikelihoodArray();
  VVVdouble* _dLikelihoods_node = &nodeData->getDLikelihoodArray();
  VVVdouble* _d2Likelihoods_node = &nodeData->getD2LikelihoodArray();
 
  _likelihoods_node->resize(nbDistinctSites_);
  _dLikelihoods_node->resize(nbDistinctSites_);
  _d2Likelihoods_node->resize(nbDistinctSites_);
 
  for (size_t i = 0; i < nbDistinctSites_; i++)
  {
    VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
    VVdouble* _dLikelihoods_node_i = &(*_dLikelihoods_node)[i];
    VVdouble* _d2Likelihoods_node_i = &(*_d2Likelihoods_node)[i];
    _likelihoods_node_i->resize(nbClasses_);
    _dLikelihoods_node_i->resize(nbClasses_);
    _d2Likelihoods_node_i->resize(nbClasses_);
    for (size_t c = 0; c < nbClasses_; c++)
    {
      Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
      Vdouble* _dLikelihoods_node_i_c = &(*_dLikelihoods_node_i)[c];
      Vdouble* _d2Likelihoods_node_i_c = &(*_d2Likelihoods_node_i)[c];
      _likelihoods_node_i_c->resize(nbStates_);
      _dLikelihoods_node_i_c->resize(nbStates_);
      _d2Likelihoods_node_i_c->resize(nbStates_);
      for (size_t s = 0; s < nbStates_; s++)
      {
        (*_likelihoods_node_i_c)[s] = 1; // All likelihoods are initialized to 1.
        (*_dLikelihoods_node_i_c)[s] = 0; // All dLikelihoods are initialized to 0.
        (*_d2Likelihoods_node_i_c)[s] = 0; // All d2Likelihoods are initialized to 0.
      }
    }
  }
 
  // Now initialize likelihood values and pointers:
  if (node->isLeaf())
  {
    size_t posSeq;
    try
    {
      posSeq = sequences.getSequencePosition(node->getName());
    }
    catch (SequenceNotFoundException& snfe)
    {
      throw SequenceNotFoundException("DRASRTreeLikelihoodData::initTreelikelihoods. Leaf name in tree not found in site container: ", (node->getName()));
    }
    for (size_t i = 0; i < nbDistinctSites_; i++)
    {
      VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
      for (size_t c = 0; c < nbClasses_; c++)
      {
        Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
        double test = 0.;
 
        for (size_t s = 0; s < nbStates_; s++)
        {
          // Leaves likelihood are set to 1 if the char correspond to the site in the sequence,
          // otherwise value set to 0:
 
          (*_likelihoods_node_i_c)[s] = sequences.getStateValueAt(i, posSeq, model.getAlphabetStateAsInt(s));
          test += (*_likelihoods_node_i_c)[s];
        }
        if (test < 0.000001)
          std::cerr << "WARNING!!! Likelihood will be 0 for site " << i << std::endl;
      }
    }
  }
  else
  {
    // 'node' is an internal node.
    std::map<int, std::vector<size_t>>* patternLinks__node = &patternLinks_[node->getId()];
    size_t nbSonNodes = node->getNumberOfSons();
    for (size_t l = 0; l < nbSonNodes; l++)
    {
      // For each son node,
      const Node* son = (*node)[static_cast<int>(l)];
      initLikelihoods(son, sequences, model);
      std::vector<size_t>* patternLinks__node_son = &(*patternLinks__node)[son->getId()];
 
      // Init map:
      patternLinks__node_son->resize(nbDistinctSites_);
 
      for (size_t i = 0; i < nbDistinctSites_; i++)
      {
        (*patternLinks__node_son)[i] = i;
      }
    }
  }
}
 
/******************************************************************************/
 
std::unique_ptr<SitePatterns> DRASRTreeLikelihoodData::initLikelihoodsWithPatterns(
    const Node* node,
    const AlignmentDataInterface& sequences,
    const TransitionModelInterface& model)
{
  vector<const Node*> leaves = TreeTemplateTools::getLeaves(*node);
 
  auto tmp = PatternTools::getSequenceSubset(sequences, *node);
 
  auto patterns = std::make_unique<SitePatterns>(*tmp);
 
  auto subSequences = patterns->getSites();
  size_t nbSites = subSequences->getNumberOfSites();
 
  // Initialize likelihood vector:
  DRASRTreeLikelihoodNodeData* nodeData = &nodeData_[node->getId()];
  nodeData->setNode(node);
  VVVdouble* _likelihoods_node = &nodeData->getLikelihoodArray();
  VVVdouble* _dLikelihoods_node = &nodeData->getDLikelihoodArray();
  VVVdouble* _d2Likelihoods_node = &nodeData->getD2LikelihoodArray();
  _likelihoods_node->resize(nbSites);
  _dLikelihoods_node->resize(nbSites);
  _d2Likelihoods_node->resize(nbSites);
 
  for (size_t i = 0; i < nbSites; i++)
  {
    VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
    VVdouble* _dLikelihoods_node_i = &(*_dLikelihoods_node)[i];
    VVdouble* _d2Likelihoods_node_i = &(*_d2Likelihoods_node)[i];
    _likelihoods_node_i->resize(nbClasses_);
    _dLikelihoods_node_i->resize(nbClasses_);
    _d2Likelihoods_node_i->resize(nbClasses_);
    for (size_t c = 0; c < nbClasses_; c++)
    {
      Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
      Vdouble* _dLikelihoods_node_i_c = &(*_dLikelihoods_node_i)[c];
      Vdouble* _d2Likelihoods_node_i_c = &(*_d2Likelihoods_node_i)[c];
      _likelihoods_node_i_c->resize(nbStates_);
      _dLikelihoods_node_i_c->resize(nbStates_);
      _d2Likelihoods_node_i_c->resize(nbStates_);
      for (size_t s = 0; s < nbStates_; s++)
      {
        (*_likelihoods_node_i_c)[s] = 1; // All likelihoods are initialized to 1.
        (*_dLikelihoods_node_i_c)[s] = 0; // All dLikelihoods are initialized to 0.
        (*_d2Likelihoods_node_i_c)[s] = 0; // All d2Likelihoods are initialized to 0.
      }
    }
  }
 
  // Now initialize likelihood values and pointers:
 
  if (node->isLeaf())
  {
    size_t posSeq;
    try
    {
      posSeq = subSequences->getSequencePosition(node->getName());
    }
    catch (SequenceNotFoundException& snfe)
    {
      throw SequenceNotFoundException("HomogeneousTreeLikelihood::initTreelikelihoodsWithPatterns. Leaf name in tree not found in site container: ", (node->getName()));
    }
 
    for (size_t i = 0; i < nbSites; i++)
    {
      VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
 
      for (size_t c = 0; c < nbClasses_; c++)
      {
        Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
        double test = 0.;
        for (size_t s = 0; s < nbStates_; s++)
        {
          // Leaves likelihood are set to 1 if the char correspond to the site in the sequence,
          // otherwise value set to 0:
          // cout << "i=" << i << "\tc=" << c << "\ts=" << s << endl;
          (*_likelihoods_node_i_c)[s] = subSequences->getStateValueAt(i, posSeq, model.getAlphabetStateAsInt(s));
          test += (*_likelihoods_node_i_c)[s];
        }
        if (test < 0.000001)
          std::cerr << "WARNING!!! Likelihood will be 0 for site " << i << std::endl;
      }
    }
  }
  else
  {
    // 'node' is an internal node.
    std::map<int, std::vector<size_t>>* patternLinks__node = &patternLinks_[node->getId()];
 
    // Now initialize pattern links:
    size_t nbSonNodes = node->getNumberOfSons();
    for (int l = 0; l < static_cast<int>(nbSonNodes); l++)
    {
      // For each son node,
      const Node* son = (*node)[l];
 
      std::vector<size_t>* patternLinks__node_son = &(*patternLinks__node)[son->getId()];
 
      // Initialize subtree 'l' and retrieves corresponding
      // subSequences:
 
      shared_ptr<SitePatterns> subPatterns(initLikelihoodsWithPatterns(son, *subSequences, model));
 
      patternLinks__node_son->resize(size_t(subPatterns->getIndices().size()));
      SitePatterns::IndicesType::Map(&((*patternLinks__node_son)[0]), subPatterns->getIndices().size()) = subPatterns->getIndices();
    }
  }
  return patterns;
}
 
/******************************************************************************/