42 #include "../PatternTools.h" 45 #include <Bpp/Seq/SiteTools.h> 53 if (sites.getNumberOfSequences() == 1)
54 throw Exception(
"Error, only 1 sequence!");
55 if (sites.getNumberOfSequences() == 0)
56 throw Exception(
"Error, no sequence!");
57 if (sites.getAlphabet()->getAlphabetType()
58 != model.getAlphabet()->getAlphabetType())
59 throw AlphabetMismatchException(
"DRASDRTreeLikelihoodData::initLikelihoods. Data and model must have the same alphabet type.",
62 alphabet_ = sites.getAlphabet();
63 nbStates_ = model.getNumberOfStates();
64 nbSites_ = sites.getNumberOfSites();
72 nbDistinctSites_ = shrunkData_->getNumberOfSites();
76 const SiteContainer* sequences =
new AlignedSequenceContainer(*shrunkData_);
77 initLikelihoods(tree_->getRootNode(), *sequences, model);
81 rootLikelihoods_.resize(nbDistinctSites_);
82 rootLikelihoodsS_.resize(nbDistinctSites_);
83 rootLikelihoodsSR_.resize(nbDistinctSites_);
84 for (
size_t i = 0; i < nbDistinctSites_; i++)
86 VVdouble* rootLikelihoods_i_ = &rootLikelihoods_[i];
87 Vdouble* rootLikelihoodsS_i_ = &rootLikelihoodsS_[i];
88 rootLikelihoods_i_->resize(nbClasses_);
89 rootLikelihoodsS_i_->resize(nbClasses_);
90 for (
size_t c = 0; c < nbClasses_; c++)
92 Vdouble* rootLikelihoods_i_c_ = &(*rootLikelihoods_i_)[c];
93 rootLikelihoods_i_c_->resize(nbStates_);
94 for (
size_t x = 0; x < nbStates_; x++)
96 (*rootLikelihoods_i_c_)[x] = 1.;
112 seq = &sites.getSequence(node->getName());
114 catch (SequenceNotFoundException& snfe)
116 throw SequenceNotFoundException(
"DRASDRTreeLikelihoodData::initlikelihoods. Leaf name in tree not found in site container: ", (node->getName()));
121 leavesLikelihoods_leaf->resize(nbDistinctSites_);
122 for (
size_t i = 0; i < nbDistinctSites_; i++)
124 Vdouble* leavesLikelihoods_leaf_i = &(*leavesLikelihoods_leaf)[i];
125 leavesLikelihoods_leaf_i->resize(nbStates_);
126 int state = seq->getValue(i);
128 for (
size_t s = 0; s < nbStates_; s++)
132 ( *leavesLikelihoods_leaf_i)[s] = model.getInitValue(s, state);
133 test += ( *leavesLikelihoods_leaf_i)[s];
136 std::cerr <<
"WARNING!!! Likelihood will be 0 for this site." << std::endl;
141 size_t nbSonNodes = node->getNumberOfSons();
142 for (
size_t l = 0; l < nbSonNodes; l++)
145 initLikelihoods(node->getSon(l), sites, model);
153 int nbSons =
static_cast<int>(node->getNumberOfSons());
155 for (
int n = (node->hasFather() ? -1 : 0); n < nbSons; n++)
157 const Node* neighbor = (*node)[n];
158 VVVdouble* likelihoods_node_neighbor_ = &(*likelihoods_node_)[neighbor->
getId()];
160 likelihoods_node_neighbor_->resize(nbDistinctSites_);
164 VVdouble* leavesLikelihoods_leaf_ = &leafData_[neighbor->
getId()].getLikelihoodArray();
165 for (
size_t i = 0; i < nbDistinctSites_; i++)
167 Vdouble* leavesLikelihoods_leaf_i_ = &(*leavesLikelihoods_leaf_)[i];
168 VVdouble* likelihoods_node_neighbor_i_ = &(*likelihoods_node_neighbor_)[i];
169 likelihoods_node_neighbor_i_->resize(nbClasses_);
170 for (
size_t c = 0; c < nbClasses_; c++)
172 Vdouble* likelihoods_node_neighbor_i_c_ = &(*likelihoods_node_neighbor_i_)[c];
173 likelihoods_node_neighbor_i_c_->resize(nbStates_);
174 for (
size_t s = 0; s < nbStates_; s++)
176 (*likelihoods_node_neighbor_i_c_)[s] = (*leavesLikelihoods_leaf_i_)[s];
183 for (
size_t i = 0; i < nbDistinctSites_; i++)
185 VVdouble* likelihoods_node_neighbor_i_ = &(*likelihoods_node_neighbor_)[i];
186 likelihoods_node_neighbor_i_->resize(nbClasses_);
187 for (
size_t c = 0; c < nbClasses_; c++)
189 Vdouble* likelihoods_node_neighbor_i_c_ = &(*likelihoods_node_neighbor_i_)[c];
190 likelihoods_node_neighbor_i_c_->resize(nbStates_);
191 for (
size_t s = 0; s < nbStates_; s++)
193 (*likelihoods_node_neighbor_i_c_)[s] = 1.;
203 dLikelihoods_node_->resize(nbDistinctSites_);
204 d2Likelihoods_node_->resize(nbDistinctSites_);
226 int nbSons =
static_cast<int>(node->getNumberOfSons());
228 for (
int n = (node->hasFather() ? -1 : 0); n < nbSons; n++)
230 const Node* neighbor = (*node)[n];
233 array->resize(nbDistinctSites_);
234 for (
size_t i = 0; i < nbDistinctSites_; i++)
236 VVdouble* array_i = &(*array)[i];
237 array_i->resize(nbClasses_);
238 for (
size_t c = 0; c < nbClasses_; c++)
240 Vdouble* array_i_c = &(*array_i)[c];
241 array_i_c->resize(nbStates_);
242 for (
size_t s = 0; s < nbStates_; s++)
244 (*array_i_c)[s] = 1.;
251 size_t nbSonNodes = node->getNumberOfSons();
252 for (
size_t l = 0; l < nbSonNodes; l++)
255 reInit(node->getSon(l));
Interface for all substitution models.
const std::vector< size_t > & getIndices() const
void setNode(const Node *node)
Set the node associated to this data.
VVdouble & getLikelihoodArray()
Likelihood data structure for a node.
const std::vector< unsigned int > & getWeights() const
VVVdouble & getLikelihoodArrayForNeighbor(int neighborId)
const TreeTemplate< Node > * tree_
Vdouble & getDLikelihoodArray()
virtual bool isLeaf() const
void reInit()
Rebuild likelihood arrays at inner nodes.
virtual int getId() const
Get the node's id.
void setNode(const Node *node)
Set the node associated to this data.
SiteContainer * getSites() const
The phylogenetic node class.
const std::map< int, VVVdouble > & getLikelihoodArrays() const
Data structure for site patterns.
void initLikelihoods(const SiteContainer &sites, const SubstitutionModel &model)
Resize and initialize all likelihood arrays according to the given data set and substitution model...
Likelihood data structure for a leaf.
void eraseNeighborArrays()
Vdouble & getD2LikelihoodArray()