42 #include "../PatternTools.h" 45 #include <Bpp/Seq/SiteTools.h> 46 #include <Bpp/Seq/Container/AlignedSequenceContainer.h> 47 #include <Bpp/Seq/Container/SequenceContainerTools.h> 48 #include <Bpp/Seq/Container/VectorSiteContainer.h> 57 if (sites.getNumberOfSequences() == 1)
58 throw Exception(
"Error, only 1 sequence!");
59 if (sites.getNumberOfSequences() == 0)
60 throw Exception(
"Error, no sequence!");
61 if (sites.getAlphabet()->getAlphabetType()
62 != model.getAlphabet()->getAlphabetType())
63 throw AlphabetMismatchException(
"DRASDRTreeLikelihoodData::initLikelihoods. Data and model must have the same alphabet type.",
66 alphabet_ = sites.getAlphabet();
67 nbStates_ = model.getNumberOfStates();
68 nbSites_ = sites.getNumberOfSites();
74 patterns = initLikelihoodsWithPatterns(tree_->getRootNode(), sites, model);
78 nbDistinctSites_ = shrunkData_->getNumberOfSites();
86 nbDistinctSites_ = shrunkData_->getNumberOfSites();
87 initLikelihoods(tree_->getRootNode(), *shrunkData_, model);
103 _likelihoods_node->resize(nbDistinctSites_);
104 _dLikelihoods_node->resize(nbDistinctSites_);
105 _d2Likelihoods_node->resize(nbDistinctSites_);
107 for (
size_t i = 0; i < nbDistinctSites_; i++)
109 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
110 VVdouble* _dLikelihoods_node_i = &(*_dLikelihoods_node)[i];
111 VVdouble* _d2Likelihoods_node_i = &(*_d2Likelihoods_node)[i];
112 _likelihoods_node_i->resize(nbClasses_);
113 _dLikelihoods_node_i->resize(nbClasses_);
114 _d2Likelihoods_node_i->resize(nbClasses_);
115 for (
size_t c = 0; c < nbClasses_; c++)
117 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
118 Vdouble* _dLikelihoods_node_i_c = &(*_dLikelihoods_node_i)[c];
119 Vdouble* _d2Likelihoods_node_i_c = &(*_d2Likelihoods_node_i)[c];
120 _likelihoods_node_i_c->resize(nbStates_);
121 _dLikelihoods_node_i_c->resize(nbStates_);
122 _d2Likelihoods_node_i_c->resize(nbStates_);
123 for (
size_t s = 0; s < nbStates_; s++)
125 (*_likelihoods_node_i_c)[s] = 1;
126 (*_dLikelihoods_node_i_c)[s] = 0;
127 (*_d2Likelihoods_node_i_c)[s] = 0;
139 seq = &sequences.getSequence(node->getName());
141 catch (SequenceNotFoundException snfe)
143 throw SequenceNotFoundException(
"DRASRTreeLikelihoodData::initTreelikelihoods. Leaf name in tree not found in site conainer: ", (node->getName()));
145 for (
size_t i = 0; i < nbDistinctSites_; i++)
147 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
148 int state = seq->getValue(i);
149 for (
size_t c = 0; c < nbClasses_; c++)
151 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
153 for (
size_t s = 0; s < nbStates_; s++)
158 (*_likelihoods_node_i_c)[s] = model.getInitValue(s, state);
159 test += (*_likelihoods_node_i_c)[s];
162 std::cerr <<
"WARNING!!! Likelihood will be 0 for this site." << std::endl;
169 std::map<int, std::vector<size_t> >* patternLinks__node = &patternLinks_[node->getId()];
170 size_t nbSonNodes = node->getNumberOfSons();
171 for (
size_t l = 0; l < nbSonNodes; l++)
174 const Node* son = (*node)[
static_cast<int>(l)];
175 initLikelihoods(son, sequences, model);
176 std::vector<size_t>* patternLinks__node_son = &(*patternLinks__node)[son->
getId()];
179 patternLinks__node_son->resize(nbDistinctSites_);
181 for (
size_t i = 0; i < nbDistinctSites_; i++)
183 (*patternLinks__node_son)[i] = i;
195 SiteContainer* subSequences = patterns->
getSites();
197 size_t nbSites = subSequences->getNumberOfSites();
205 _likelihoods_node->resize(nbSites);
206 _dLikelihoods_node->resize(nbSites);
207 _d2Likelihoods_node->resize(nbSites);
209 for (
size_t i = 0; i < nbSites; i++)
211 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
212 VVdouble* _dLikelihoods_node_i = &(*_dLikelihoods_node)[i];
213 VVdouble* _d2Likelihoods_node_i = &(*_d2Likelihoods_node)[i];
214 _likelihoods_node_i->resize(nbClasses_);
215 _dLikelihoods_node_i->resize(nbClasses_);
216 _d2Likelihoods_node_i->resize(nbClasses_);
217 for (
size_t c = 0; c < nbClasses_; c++)
219 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
220 Vdouble* _dLikelihoods_node_i_c = &(*_dLikelihoods_node_i)[c];
221 Vdouble* _d2Likelihoods_node_i_c = &(*_d2Likelihoods_node_i)[c];
222 _likelihoods_node_i_c->resize(nbStates_);
223 _dLikelihoods_node_i_c->resize(nbStates_);
224 _d2Likelihoods_node_i_c->resize(nbStates_);
225 for (
size_t s = 0; s < nbStates_; s++)
227 (*_likelihoods_node_i_c)[s] = 1;
228 (*_dLikelihoods_node_i_c)[s] = 0;
229 (*_d2Likelihoods_node_i_c)[s] = 0;
241 seq = &subSequences->getSequence(node->getName());
243 catch (SequenceNotFoundException snfe)
245 throw SequenceNotFoundException(
"HomogeneousTreeLikelihood::initTreelikelihoodsWithPatterns. Leaf name in tree not found in site conainer: ", (node->getName()));
247 for (
size_t i = 0; i < nbSites; i++)
249 VVdouble* _likelihoods_node_i = &(*_likelihoods_node)[i];
250 int state = seq->getValue(i);
251 for (
size_t c = 0; c < nbClasses_; c++)
253 Vdouble* _likelihoods_node_i_c = &(*_likelihoods_node_i)[c];
255 for (
size_t s = 0; s < nbStates_; s++)
260 (*_likelihoods_node_i_c)[s] = model.getInitValue(s, state);
261 test += (*_likelihoods_node_i_c)[s];
264 std::cerr <<
"WARNING!!! Likelihood will be 0 for this site." << std::endl;
271 std::map<int, std::vector<size_t> >* patternLinks__node = &patternLinks_[node->getId()];
274 size_t nbSonNodes = node->getNumberOfSons();
275 for (
int l = 0; l < static_cast<int>(nbSonNodes); l++)
278 const Node* son = (*node)[l];
280 std::vector<size_t>* patternLinks__node_son = &(*patternLinks__node)[son->
getId()];
283 SitePatterns* subPatterns = initLikelihoodsWithPatterns(son, *subSequences, model);
284 (*patternLinks__node_son) = subPatterns->
getIndices();
Interface for all substitution models.
const std::vector< size_t > & getIndices() const
void setNode(const Node *node)
Set the node associated to this data.
const std::vector< unsigned int > & getWeights() const
VVVdouble & getDLikelihoodArray()
virtual int getId() const
Get the node's id.
SiteContainer * getSites() const
void initLikelihoods(const SiteContainer &sites, const SubstitutionModel &model)
The phylogenetic node class.
VVVdouble & getLikelihoodArray()
VVVdouble & getD2LikelihoodArray()
Data structure for site patterns.
Likelihood data structure for a node.
virtual SitePatterns * initLikelihoodsWithPatterns(const Node *node, const SiteContainer &sequences, const SubstitutionModel &model)
This method initializes the leaves according to a sequence file.