41 #include <Bpp/Numeric/VectorTools.h> 42 #include <Bpp/Numeric/Random/RandomTools.h> 49 vector<size_t> ancestors(nbDistinctSites_);
50 probs.resize(nbDistinctSites_);
53 if (likelihood_->getTree().isLeaf(nodeId))
55 VVdouble larray = likelihood_->getLikelihoodData()->getLeafLikelihoods(nodeId);
56 for (
size_t i = 0; i < nbDistinctSites_; ++i)
58 Vdouble* probs_i = &probs[i];
59 probs_i->resize(nbStates_);
60 size_t j = VectorTools::whichMax(larray[i]);
69 likelihood_->computeLikelihoodAtNode(nodeId, larray);
70 for (
size_t i = 0; i < nbDistinctSites_; i++)
72 VVdouble* larray_i = &larray[i];
73 Vdouble* probs_i = &probs[i];
74 probs_i->resize(nbStates_);
75 for (
size_t c = 0; c < nbClasses_; c++)
77 Vdouble* larray_i_c = &(*larray_i)[c];
78 for (
size_t x = 0; x < nbStates_; x++)
80 (*probs_i)[x] += (*larray_i_c)[x] * r_[c] / l_[i];
86 r = RandomTools::giveRandomNumberBetweenZeroAndEntry(1.);
87 for (
size_t j = 0; j < nbStates_; j++)
89 cumProb += (*probs_i)[j];
98 ancestors[i] = VectorTools::whichMax(*probs_i);
106 map<int, vector<size_t> > ancestors;
108 AlignedSequenceContainer* data =
new AlignedSequenceContainer(*likelihood_->getLikelihoodData()->getShrunkData());
109 recursiveMarginalAncestralStates(tree_.getRootNode(), ancestors, *data);
116 string name = tree_.hasNodeName(nodeId) ? tree_.getNodeName(nodeId) : (
"" + TextTools::toString(nodeId));
117 const vector<size_t>* rootPatternLinks = &likelihood_->getLikelihoodData()->getRootArrayPositions();
118 const SubstitutionModel* model = likelihood_->getSubstitutionModel(tree_.getNodesId()[0], 0);
119 vector<size_t> states;
120 vector<int> allStates(nbSites_);
121 VVdouble patternedProbs;
124 states = getAncestralStatesForNode(nodeId, patternedProbs, sample);
125 probs->resize(nbSites_);
126 for (
size_t i = 0; i < nbSites_; i++)
129 (*probs)[i] = patternedProbs[(*rootPatternLinks)[i]];
134 states = getAncestralStatesForNode(nodeId, patternedProbs, sample);
135 for (
size_t i = 0; i < nbSites_; i++)
140 return new BasicSequence(name, allStates, alphabet_);
145 map<
int, vector<size_t> >& ancestors,
146 AlignedSequenceContainer& data)
const 150 vector<int> content = data.getContent(node->
getName());
151 vector<size_t>* v = &ancestors[node->
getId()];
152 v->resize(content.size());
156 const SubstitutionModel* model = likelihood_->getSubstitutionModel(tree_.getNodesId()[0], 0);
157 for (
size_t i = 0; i < content.size(); i++)
164 ancestors[node->
getId()] = getAncestralStatesForNode(node->
getId());
167 recursiveMarginalAncestralStates(node->
getSon(i), ancestors, data);
174 AlignedSequenceContainer* asc =
new AlignedSequenceContainer(alphabet_);
175 vector<int> ids = tree_.getInnerNodesId();
176 for (
size_t i = 0; i < ids.size(); i++)
178 Sequence* seq = getAncestralSequenceForNode(ids[i], NULL, sample);
179 asc->addSequence(*seq);
Interface for all substitution models.
virtual std::vector< size_t > getModelStates(int code) const =0
Get the state in the model corresponding to a particular state in the alphabet.
virtual int getAlphabetStateAsInt(size_t index) const =0
virtual std::string getName() const
Get the name associated to this node, if there is one, otherwise throw a NodeException.
std::map< int, std::vector< size_t > > getAllAncestralStates() const
Get all ancestral states for all nodes.
virtual const Node * getSon(size_t pos) const
void recursiveMarginalAncestralStates(const Node *node, std::map< int, std::vector< size_t > > &ancestors, AlignedSequenceContainer &data) const
AlignedSequenceContainer * getAncestralSequences() const
Get all the ancestral sequences for all nodes.
virtual bool isLeaf() const
virtual int getId() const
Get the node's id.
The phylogenetic node class.
virtual size_t getNumberOfSons() const
Sequence * getAncestralSequenceForNode(int nodeId, VVdouble *probs, bool sample) const
Get the ancestral sequence for a given node.
std::vector< size_t > getAncestralStatesForNode(int nodeId, VVdouble &probs, bool sample) const
Get ancestral states for a given node as a vector of int.