10 #include "../Likelihood/DataFlow/ForwardLikelihoodTree.h"
11 #include "../Likelihood/DataFlow/LikelihoodCalculationSingleProcess.h"
15 using namespace numeric;
27 const vector<uint>& edgeIds,
29 short unresolvedOption,
34 throw Exception(
"RewardMappingTools::computeSubstitutionVectors(). Likelihood object is not initialized.");
38 if (edgeIds.size() == 0)
39 return make_unique<ProbabilisticRewardMapping>(
48 map<const SubstitutionModelInterface*, shared_ptr<Reward>> mModReward;
50 for (
auto speciesId : edgeIds)
52 const auto& dagIndexes = rltc.
getEdgesIds(speciesId, 0);
54 for (
auto id : dagIndexes)
56 const auto& edge = processTree->getEdge(
id);
59 auto model = edge->getModel();
61 auto nMod = edge->getNMod();
63 auto tm = dynamic_pointer_cast<const TransitionModelInterface>(model->targetValue());
65 shared_ptr<const SubstitutionModelInterface> sm =
nullptr;
69 sm = dynamic_pointer_cast<const SubstitutionModelInterface>(tm);
72 throw Exception(
"SubstitutionMappingTools::computeCounts : SubstitutionVectors possible only for SubstitutionModels, not in branch " +
TextTools::toString(speciesId) +
". Got model " + tm->getName());
76 size_t nmod = nMod->targetValue();
78 auto ttm = dynamic_pointer_cast<const MixedTransitionModelInterface>(tm);
80 throw Exception(
"SubstitutionMappingTools::computeCounts : Expecting Mixed model in branch " +
TextTools::toString(speciesId) +
". Got model " + tm->getName());
82 sm = dynamic_pointer_cast<const SubstitutionModelInterface>(ttm->getNModel(nmod));
88 if (mModReward.find(sm.get()) == mModReward.end())
90 mModReward[sm.get()] = shared_ptr<Reward>(reward.
clone());
91 mModReward[sm.get()]->setSubstitutionModel(sm);
103 size_t nbNodes = edgeIds.size();
118 unique_ptr<ProbabilisticRewardMapping::mapTree::EdgeIterator> brIt = rewards->allEdgesIterator();
120 Eigen::MatrixXd rpxy;
123 for ( ; !brIt->end(); brIt->next())
128 shared_ptr<PhyloBranchReward> br = **brIt;
131 uint speciesId = rewards->getEdgeIndex(br);
138 for (
size_t ncl = 0; ncl < nbClasses; ncl++)
146 const auto& dagIndexes = rltc.
getEdgesIds(speciesId, ncl);
149 for (
auto id : dagIndexes)
151 auto edge = processTree->getEdge(
id);
153 auto tm = dynamic_pointer_cast<const TransitionModelInterface>(edge->model().targetValue());
155 auto nMod = edge->getNMod();
157 shared_ptr<const SubstitutionModelInterface> sm =
nullptr;
160 sm = dynamic_pointer_cast<const SubstitutionModelInterface>(tm);
163 size_t nmod = nMod->targetValue();
165 auto ttm = dynamic_pointer_cast<const MixedTransitionModelInterface>(tm);
166 sm = dynamic_pointer_cast<const SubstitutionModelInterface>(ttm->getNModel(nmod));
169 auto subReward = mModReward[sm.get()];
180 const Eigen::MatrixXd& pxy = edge->getTransitionMatrix()->accessValueConst();
184 subReward->storeAllRewards(edge->getBrLen()->getValue(), rpxy);
186 rpxy.array() *= pxy.array();
190 auto rew = rpxy * likelihoodsBotEdge;
192 auto bb = (
cwise(likelihoodsTopEdge) *
cwise(rew)).colwise().sum();
195 Eigen::VectorXd ff(likelihoodsBotEdge.cols());
196 switch (unresolvedOption)
202 for (
auto i = 0; i < ff.size(); i++)
204 const auto& s = likelihoodsBotEdge.col(i).sum();
219 auto cc = bb /
cwise(likelihoodsFather);
225 rewardsForCurrentNode += rewardsForCurrentClass * pr;
229 for (
size_t i = 0; i < nbDistinctSites; ++i)
231 (*br)(i) =
convert(rewardsForCurrentNode(Eigen::Index(i)));
251 throw IOException(
"RewardMappingTools::writeToFile. Can't write to stream.");
260 unique_ptr<ProbabilisticRewardMapping::mapTree::EdgeIterator> brIt = rewards.
allEdgesIterator();
262 for ( ; !brIt->end(); brIt->next())
264 const shared_ptr<PhyloBranchReward> br = **brIt;
266 out << rewards.
getEdgeIndex(br) <<
"\t" << br->getLength();
270 out <<
"\t" << br->getSiteReward(rewards.
getSiteIndex(i));
284 vector<string> ids = data->getColumn(0);
285 data->deleteColumn(0);
286 data->deleteColumn(0);
288 size_t nbSites = data->getNumberOfColumns();
291 unique_ptr<ProbabilisticRewardMapping::mapTree::EdgeIterator> brIt = rewards.
allEdgesIterator();
293 for ( ; !brIt->end(); brIt->next())
295 const shared_ptr<PhyloBranchReward> br = **brIt;
298 for (
size_t j = 0; j < nbSites; j++)
305 for (
size_t i = 0; i < nbSites; i++)
307 string siteTxt = data->getColumnName(i);
309 if (siteTxt.substr(0, 4) ==
"Site")
310 site = TextTools::to<int>(siteTxt.substr(4));
312 site = TextTools::to<int>(siteTxt);
328 shared_ptr<PhyloBranchReward> br = smap.
getEdge((uint)branchIndex);
330 for (
size_t i = 0; i < nbSites; ++i)
342 unique_ptr<ProbabilisticRewardMapping::mapTree::EdgeIterator> brIt = smap.
allEdgesIterator();
346 for ( ; !brIt->end(); brIt->next())
348 v += (**brIt)->getSiteReward(siteIndex);
void setSitePosition(size_t index, int position)
Set the position of a given site.
size_t getNumberOfSites() const
virtual std::unique_ptr< EdgeIterator > allEdgesIterator()=0
virtual EdgeIndex getEdgeIndex(const std::shared_ptr< E > edgeObject) const=0
virtual std::shared_ptr< E > getEdge(EdgeIndex edgeIndex) const=0
virtual int getCoordinate() const=0
static std::unique_ptr< DataTable > read(std::istream &in, const std::string &sep="\t", bool header=true, int rowNames=-1)
const char * what() const noexcept override
static Self Zero(Eigen::Index rows, Eigen::Index cols)
SiteLikelihoodsRef getLikelihoodsAtNodeForClass(uint nodeId, size_t nCat)
Get shrunked conditional likelihood matrix at Node (ie just above the node), for a given rate class.
const SubstitutionProcessInterface & substitutionProcess() const
Return the ref to the SubstitutionProcess.
std::shared_ptr< ForwardLikelihoodTree > getForwardLikelihoodTree(size_t nCat)
ConditionalLikelihoodRef getBackwardLikelihoodsAtEdgeForClass(uint edgeId, size_t nCat)
ConditionalLikelihoodRef getForwardLikelihoodsAtNodeForClass(uint nodeId, size_t nCat)
const PatternType & getRootArrayPositions() const
size_t getNumberOfDistinctSites() const
const DAGindexes & getEdgesIds(uint speciesId, size_t nCat) const
Get indexes of the non-empty edges in the Likelihood DAG that have a given species index for a given ...
std::shared_ptr< ProcessTree > getTreeNode(size_t nCat)
Get process tree for a rate category.
virtual bool isInitialized() const
Data storage class for probabilistic rewards mappings.
virtual void setNumberOfSites(size_t numberOfSites) override
const size_t getSiteIndex(size_t site) const
double getProbaAtNode(PhyloTree::NodeIndex nodeId)
virtual Reward * clone() const =0
This interface describes the substitution process along the tree and sites of the alignment.
virtual const ParametrizablePhyloTree & parametrizablePhyloTree() const =0
virtual size_t getNumberOfClasses() const =0
virtual double getProbabilityForModel(size_t classIndex) const =0
virtual std::shared_ptr< const ParametrizablePhyloTree > getParametrizablePhyloTree() const =0
virtual const CoreSiteInterface & site(size_t siteIndex) const=0
double toDouble(const std::string &s, char dec='.', char scientificNotation='e')
std::string toString(T t)
R convert(const F &from, const Dimension< R > &)
Defines the basic types of data flow nodes.