← Back to branch summary

~ubuntu-branches/ubuntu/precise/libbpp-phyl/precise

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Likelihood/DRASDRTreeLikelihoodData.h

  • Committer: Bazaar Package Importer
  • Author(s): Julien Dutheil
  • Date: 2011-06-09 11:00:00 UTC
  • Revision ID: james.westby@ubuntu.com-20110609110000-yvx78svv6w7xxgph
Tags: upstream-2.0.2
Import upstream version 2.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Bpp/Phyl/Likelihood/DRASDRTreeLikelihoodData.h
 
1
//
 
2
// File: DRASDRTreeLikelihoodData.h
 
3
// Created by: Julien Dutheil
 
4
// Created on: Sat Dec 30 14:20 2006
 
5
// From file DRHomogeneousTreeLikelihood.h
 
6
//
 
7
 
 
8
/*
 
9
Copyright or © or Copr. CNRS, (November 16, 2004)
 
10
 
 
11
This software is a computer program whose purpose is to provide classes
 
12
for phylogenetic data analysis.
 
13
 
 
14
This software is governed by the CeCILL  license under French law and
 
15
abiding by the rules of distribution of free software.  You can  use, 
 
16
modify and/ or redistribute the software under the terms of the CeCILL
 
17
license as circulated by CEA, CNRS and INRIA at the following URL
 
18
"http://www.cecill.info". 
 
19
 
 
20
As a counterpart to the access to the source code and  rights to copy,
 
21
modify and redistribute granted by the license, users are provided only
 
22
with a limited warranty  and the software's author,  the holder of the
 
23
economic rights,  and the successive licensors  have only  limited
 
24
liability. 
 
25
 
 
26
In this respect, the user's attention is drawn to the risks associated
 
27
with loading,  using,  modifying and/or developing or reproducing the
 
28
software by the user in light of its specific status of free software,
 
29
that may mean  that it is complicated to manipulate,  and  that  also
 
30
therefore means  that it is reserved for developers  and  experienced
 
31
professionals having in-depth computer knowledge. Users are therefore
 
32
encouraged to load and test the software's suitability as regards their
 
33
requirements in conditions enabling the security of their systems and/or 
 
34
data to be ensured and,  more generally, to use and operate it in the 
 
35
same conditions as regards security. 
 
36
 
 
37
The fact that you are presently reading this means that you have had
 
38
knowledge of the CeCILL license and that you accept its terms.
 
39
*/
 
40
 
 
41
#ifndef _DRASDRHOMOGENEOUSTREELIKELIHOODDATA_H_
 
42
#define _DRASDRHOMOGENEOUSTREELIKELIHOODDATA_H_
 
43
 
 
44
#include "AbstractTreeLikelihoodData.h"
 
45
#include "../Model/SubstitutionModel.h"
 
46
#include "../PatternTools.h"
 
47
#include "../SitePatterns.h"
 
48
 
 
49
//From SeqLib:
 
50
#include <Bpp/Seq/Container/AlignedSequenceContainer.h>
 
51
 
 
52
// From the STL:
 
53
#include <map>
 
54
 
 
55
namespace bpp
 
56
{
 
57
 
 
58
/**
 
59
 * @brief Likelihood data structure for a leaf.
 
60
 * 
 
61
 * This class is for use with the DRASDRTreeLikelihoodData class.
 
62
 * 
 
63
 * Store the likelihoods arrays associated to a leaf.
 
64
 * 
 
65
 * @see DRASDRTreeLikelihoodData
 
66
 */
 
67
class DRASDRTreeLikelihoodLeafData :
 
68
  public virtual TreeLikelihoodNodeData
 
69
{
 
70
  private:
 
71
    mutable VVdouble leafLikelihood_;
 
72
    const Node* leaf_;
 
73
 
 
74
  public:
 
75
    DRASDRTreeLikelihoodLeafData() : leafLikelihood_(), leaf_(0) {}
 
76
 
 
77
    DRASDRTreeLikelihoodLeafData(const DRASDRTreeLikelihoodLeafData& data) :
 
78
      leafLikelihood_(data.leafLikelihood_), leaf_(data.leaf_) {}
 
79
    
 
80
    DRASDRTreeLikelihoodLeafData& operator=(const DRASDRTreeLikelihoodLeafData& data)
 
81
    {
 
82
      leafLikelihood_ = data.leafLikelihood_;
 
83
      leaf_           = data.leaf_;
 
84
      return *this;
 
85
    }
 
86
 
 
87
#ifndef NO_VIRTUAL_COV
 
88
    DRASDRTreeLikelihoodLeafData*
 
89
#else
 
90
    Clonable*
 
91
#endif
 
92
    clone() const
 
93
    { 
 
94
      return new DRASDRTreeLikelihoodLeafData(*this);
 
95
    }
 
96
 
 
97
  public:
 
98
    const Node* getNode() const { return leaf_; }
 
99
    void setNode(const Node* node) { leaf_ = node; }
 
100
 
 
101
    VVdouble& getLikelihoodArray()  { return leafLikelihood_;  }
 
102
};
 
103
 
 
104
/**
 
105
 * @brief Likelihood data structure for a node.
 
106
 * 
 
107
 * This class is for use with the DRASDRTreeLikelihoodData class.
 
108
 * 
 
109
 * Store for each neighbor node an array with conditionnal likelihoods.
 
110
 *
 
111
 * @see DRASDRTreeLikelihoodData
 
112
 */
 
113
class DRASDRTreeLikelihoodNodeData :
 
114
  public virtual TreeLikelihoodNodeData
 
115
{
 
116
  private:
 
117
    /**
 
118
     * @brief This contains all likelihood values used for computation.
 
119
     *
 
120
     * <pre>
 
121
     * x[b][i][c][s]
 
122
     *   |------------> Neighbor node of n (id)
 
123
      *     |---------> Site i
 
124
     *         |------> Rate class c
 
125
     *            |---> Ancestral state s
 
126
     * </pre>
 
127
     * We call this the <i>likelihood array</i> for each node.
 
128
     */
 
129
 
 
130
    mutable std::map<int, VVVdouble> nodeLikelihoods_;
 
131
    /**
 
132
     * @brief This contains all likelihood first order derivatives values used for computation.
 
133
     *
 
134
     * <pre>
 
135
     * x[i]
 
136
     *   |---------> Site i
 
137
     * </pre> 
 
138
     * We call this the <i>dLikelihood array</i> for each node.
 
139
     */
 
140
    mutable Vdouble nodeDLikelihoods_;
 
141
  
 
142
    /**
 
143
     * @brief This contains all likelihood second order derivatives values used for computation.
 
144
     *
 
145
     * <pre>
 
146
     * x[i]
 
147
         |---------> Site i
 
148
     * </pre> 
 
149
     * We call this the <i>d2Likelihood array</i> for each node.
 
150
     */
 
151
    mutable Vdouble nodeD2Likelihoods_;
 
152
    
 
153
    const Node* node_;
 
154
 
 
155
  public:
 
156
    DRASDRTreeLikelihoodNodeData() : nodeLikelihoods_(), nodeDLikelihoods_(), nodeD2Likelihoods_(), node_(0) {}
 
157
    
 
158
    DRASDRTreeLikelihoodNodeData(const DRASDRTreeLikelihoodNodeData& data) :
 
159
      nodeLikelihoods_(data.nodeLikelihoods_),
 
160
      nodeDLikelihoods_(data.nodeDLikelihoods_),
 
161
      nodeD2Likelihoods_(data.nodeD2Likelihoods_),
 
162
      node_(data.node_)
 
163
    {}
 
164
    
 
165
    DRASDRTreeLikelihoodNodeData& operator=(const DRASDRTreeLikelihoodNodeData& data)
 
166
    {
 
167
      nodeLikelihoods_   = data.nodeLikelihoods_;
 
168
      nodeDLikelihoods_  = data.nodeDLikelihoods_;
 
169
      nodeD2Likelihoods_ = data.nodeD2Likelihoods_;
 
170
      node_              = data.node_;
 
171
      return *this;
 
172
    }
 
173
 
 
174
    virtual ~DRASDRTreeLikelihoodNodeData() {}
 
175
 
 
176
#ifndef NO_VIRTUAL_COV
 
177
    DRASDRTreeLikelihoodNodeData*
 
178
#else 
 
179
    Clonable*
 
180
#endif
 
181
    clone() const
 
182
    { 
 
183
      return new DRASDRTreeLikelihoodNodeData(*this);
 
184
    }
 
185
 
 
186
  public:
 
187
    const Node* getNode() const { return node_; }
 
188
    
 
189
    void setNode(const Node* node) { node_ = node; }
 
190
 
 
191
    const std::map<int, VVVdouble>& getLikelihoodArrays() const { return nodeLikelihoods_; }
 
192
    
 
193
    std::map<int, VVVdouble>& getLikelihoodArrays() { return nodeLikelihoods_; }
 
194
    
 
195
    VVVdouble& getLikelihoodArrayForNeighbor(int neighborId)
 
196
    {
 
197
      return nodeLikelihoods_[neighborId];
 
198
    }
 
199
    
 
200
    const VVVdouble& getLikelihoodArrayForNeighbor(int neighborId) const
 
201
    {
 
202
      return nodeLikelihoods_[neighborId];
 
203
    }
 
204
    
 
205
    Vdouble& getDLikelihoodArray() { return nodeDLikelihoods_;  }
 
206
    
 
207
    const Vdouble& getDLikelihoodArray() const  {  return nodeDLikelihoods_;  }
 
208
    
 
209
    Vdouble& getD2LikelihoodArray()  {  return nodeD2Likelihoods_; }
 
210
    
 
211
    const Vdouble& getD2LikelihoodArrayForNeighbor() const  { return nodeD2Likelihoods_; }
 
212
 
 
213
    bool isNeighbor(int neighborId) const
 
214
    {
 
215
      return nodeLikelihoods_.find(neighborId) != nodeLikelihoods_.end();
 
216
    }
 
217
 
 
218
    void eraseNeighborArrays()
 
219
    {
 
220
      nodeLikelihoods_.erase(nodeLikelihoods_.begin(), nodeLikelihoods_.end());
 
221
      nodeDLikelihoods_.erase(nodeDLikelihoods_.begin(), nodeDLikelihoods_.end());
 
222
      nodeD2Likelihoods_.erase(nodeD2Likelihoods_.begin(), nodeD2Likelihoods_.end());
 
223
    }
 
224
};
 
225
 
 
226
/**
 
227
 * @brief Likelihood data structure for rate across sites models, using a double-recursive algorithm.
 
228
 */
 
229
class DRASDRTreeLikelihoodData :
 
230
  public virtual AbstractTreeLikelihoodData
 
231
{
 
232
  private:
 
233
 
 
234
    mutable std::map<int, DRASDRTreeLikelihoodNodeData> nodeData_;
 
235
    mutable std::map<int, DRASDRTreeLikelihoodLeafData> leafData_;
 
236
    mutable VVVdouble rootLikelihoods_;
 
237
    mutable VVdouble  rootLikelihoodsS_;
 
238
    mutable Vdouble   rootLikelihoodsSR_;
 
239
 
 
240
    SiteContainer* shrunkData_;
 
241
    unsigned int nbSites_; 
 
242
    unsigned int nbStates_;
 
243
    unsigned int nbClasses_;
 
244
    unsigned int nbDistinctSites_; 
 
245
 
 
246
  public:
 
247
    DRASDRTreeLikelihoodData(const TreeTemplate<Node>* tree, unsigned int nbClasses) :
 
248
      AbstractTreeLikelihoodData(tree),
 
249
      nodeData_(), leafData_(), rootLikelihoods_(), rootLikelihoodsS_(), rootLikelihoodsSR_(),
 
250
      shrunkData_(0), nbSites_(0), nbStates_(0), nbClasses_(nbClasses), nbDistinctSites_(0)
 
251
    {}
 
252
 
 
253
    DRASDRTreeLikelihoodData(const DRASDRTreeLikelihoodData& data):
 
254
      AbstractTreeLikelihoodData(data),
 
255
      nodeData_(data.nodeData_), leafData_(data.leafData_),
 
256
      rootLikelihoods_(data.rootLikelihoods_),
 
257
      rootLikelihoodsS_(data.rootLikelihoodsS_),
 
258
      rootLikelihoodsSR_(data.rootLikelihoodsSR_),
 
259
      shrunkData_(0),
 
260
      nbSites_(data.nbSites_), nbStates_(data.nbStates_),
 
261
      nbClasses_(data.nbClasses_), nbDistinctSites_(data.nbDistinctSites_)
 
262
    {
 
263
      if (data.shrunkData_)
 
264
        shrunkData_ = dynamic_cast<SiteContainer*>(data.shrunkData_->clone());
 
265
    }
 
266
 
 
267
    DRASDRTreeLikelihoodData& operator=(const DRASDRTreeLikelihoodData& data)
 
268
    {
 
269
      AbstractTreeLikelihoodData::operator=(data);
 
270
      nodeData_          = data.nodeData_;
 
271
      leafData_          = data.leafData_;
 
272
      rootLikelihoods_   = data.rootLikelihoods_;
 
273
      rootLikelihoodsS_  = data.rootLikelihoodsS_;
 
274
      rootLikelihoodsSR_ = data.rootLikelihoodsSR_;
 
275
      nbSites_           = data.nbSites_;
 
276
      nbStates_          = data.nbStates_;
 
277
      nbClasses_         = data.nbClasses_;
 
278
      nbDistinctSites_   = data.nbDistinctSites_;
 
279
      if (shrunkData_) delete shrunkData_;
 
280
      if (data.shrunkData_)
 
281
        shrunkData_      = dynamic_cast<SiteContainer *>(data.shrunkData_->clone());
 
282
      else
 
283
        shrunkData_      = 0;
 
284
      return *this;
 
285
    }
 
286
 
 
287
    virtual ~DRASDRTreeLikelihoodData() { delete shrunkData_; }
 
288
 
 
289
    DRASDRTreeLikelihoodData* clone() const { return new DRASDRTreeLikelihoodData(*this); }
 
290
 
 
291
  public:
 
292
    /**
 
293
     * @brief Set the tree associated to the data.
 
294
     *
 
295
     * All node data will be actualized accordingly by calling the setNode() method on the corresponding nodes.
 
296
     * @warning: the old tree and the new tree must be two clones! And particularly, they have to share the
 
297
     * same topology and nodes id.
 
298
     *
 
299
     * @param tree The tree to be associated to this data.
 
300
     */
 
301
    void setTree(const TreeTemplate<Node>* tree)
 
302
    { 
 
303
      tree_ = tree;
 
304
      for (std::map<int, DRASDRTreeLikelihoodNodeData>::iterator it = nodeData_.begin(); it != nodeData_.end(); it++)
 
305
      {
 
306
        int id = it->second.getNode()->getId();
 
307
        it->second.setNode(tree_->getNode(id));
 
308
      }
 
309
      for (std::map<int, DRASDRTreeLikelihoodLeafData>::iterator it = leafData_.begin(); it != leafData_.end(); it++)
 
310
      {
 
311
        int id = it->second.getNode()->getId();
 
312
        it->second.setNode(tree_->getNode(id));
 
313
      }
 
314
    }
 
315
 
 
316
    DRASDRTreeLikelihoodNodeData& getNodeData(int nodeId)
 
317
    { 
 
318
      return nodeData_[nodeId];
 
319
    }
 
320
    
 
321
    const DRASDRTreeLikelihoodNodeData& getNodeData(int nodeId) const
 
322
    { 
 
323
      return nodeData_[nodeId];
 
324
    }
 
325
    
 
326
    DRASDRTreeLikelihoodLeafData& getLeafData(int nodeId)
 
327
    { 
 
328
      return leafData_[nodeId];
 
329
    }
 
330
    
 
331
    const DRASDRTreeLikelihoodLeafData& getLeafData(int nodeId) const
 
332
    { 
 
333
      return leafData_[nodeId];
 
334
    }
 
335
    
 
336
    unsigned int getArrayPosition(int parentId, int sonId, unsigned int currentPosition) const
 
337
    {
 
338
      return currentPosition;
 
339
    }
 
340
 
 
341
    const std::map<int, VVVdouble>& getLikelihoodArrays(int nodeId) const 
 
342
    {
 
343
      return nodeData_[nodeId].getLikelihoodArrays();
 
344
    }
 
345
    
 
346
    std::map<int, VVVdouble>& getLikelihoodArrays(int nodeId)
 
347
    {
 
348
      return nodeData_[nodeId].getLikelihoodArrays();
 
349
    }
 
350
 
 
351
    VVVdouble& getLikelihoodArray(int parentId, int neighborId)
 
352
    {
 
353
      return nodeData_[parentId].getLikelihoodArrayForNeighbor(neighborId);
 
354
    }
 
355
    
 
356
    const VVVdouble& getLikelihoodArray(int parentId, int neighborId) const
 
357
    {
 
358
      return nodeData_[parentId].getLikelihoodArrayForNeighbor(neighborId);
 
359
    }
 
360
    
 
361
    Vdouble& getDLikelihoodArray(int nodeId)
 
362
    {
 
363
      return nodeData_[nodeId].getDLikelihoodArray();
 
364
    }
 
365
    
 
366
    const Vdouble& getDLikelihoodArray(int nodeId) const
 
367
    {
 
368
      return nodeData_[nodeId].getDLikelihoodArray();
 
369
    }
 
370
    
 
371
    Vdouble& getD2LikelihoodArray(int nodeId)
 
372
    {
 
373
      return nodeData_[nodeId].getD2LikelihoodArray();
 
374
    }
 
375
 
 
376
    const Vdouble& getD2LikelihoodArray(int nodeId) const
 
377
    {
 
378
      return nodeData_[nodeId].getD2LikelihoodArray();
 
379
    }
 
380
 
 
381
    VVdouble& getLeafLikelihoods(int nodeId)
 
382
    {
 
383
      return leafData_[nodeId].getLikelihoodArray();
 
384
    }
 
385
    
 
386
    const VVdouble& getLeafLikelihoods(int nodeId) const
 
387
    {
 
388
      return leafData_[nodeId].getLikelihoodArray();
 
389
    }
 
390
    
 
391
    VVVdouble& getRootLikelihoodArray() { return rootLikelihoods_; }
 
392
    const VVVdouble & getRootLikelihoodArray() const { return rootLikelihoods_; }
 
393
    
 
394
    VVdouble& getRootSiteLikelihoodArray() { return rootLikelihoodsS_; }
 
395
    const VVdouble& getRootSiteLikelihoodArray() const { return rootLikelihoodsS_; }
 
396
    
 
397
    Vdouble& getRootRateSiteLikelihoodArray() { return rootLikelihoodsSR_; }
 
398
    const Vdouble& getRootRateSiteLikelihoodArray() const { return rootLikelihoodsSR_; }
 
399
 
 
400
    unsigned int getNumberOfDistinctSites() const { return nbDistinctSites_; }
 
401
    
 
402
    unsigned int getNumberOfSites() const { return nbSites_; }
 
403
    
 
404
    unsigned int getNumberOfStates() const { return nbStates_; }
 
405
    
 
406
    unsigned int getNumberOfClasses() const { return nbClasses_; }
 
407
 
 
408
    const SiteContainer* getShrunkData() const { return shrunkData_; }
 
409
    
 
410
    /**
 
411
     * @brief Resize and initialize all likelihood arrays according to the given data set and substitution model.
 
412
     *
 
413
     * @param sites The sequences to use as data.
 
414
     * @param model The substitution model to use.
 
415
     * @throw Exception if an error occures.
 
416
     */
 
417
    void initLikelihoods(const SiteContainer& sites, const SubstitutionModel& model) throw (Exception);
 
418
    
 
419
    /**
 
420
     * @brief Rebuild likelihood arrays at inner nodes.
 
421
     *
 
422
     * This method is to be called when the topology of the tree has changed.
 
423
     * Node arrays relationship are rebuilt according to the new topology of the tree.
 
424
     * The leaves likelihood remain unchanged, so as for the first and second order derivatives.
 
425
     */
 
426
    void reInit() throw (Exception);
 
427
    
 
428
    void reInit(const Node* node) throw (Exception);
 
429
 
 
430
  protected:
 
431
    /**
 
432
     * @brief This method initializes the leaves according to a sequence container.
 
433
     *
 
434
     * Here the container shrunkData_ is used.
 
435
     * Likelihood is set to 1 for the state corresponding to the sequence site,
 
436
     * otherwise it is set to 0.
 
437
     *
 
438
     * All likelihood arrays at each nodes are initialized according to alphabet
 
439
     * size and sequences length, and filled with 1.
 
440
     *
 
441
     * NB: This method is recursive.
 
442
     *
 
443
     * @param node  The node defining the subtree to analyse.
 
444
     * @param sites The sequence container to use.
 
445
     * @param model The model, used for initializing leaves' likelihoods.
 
446
     */
 
447
    void initLikelihoods(const Node* node, const SiteContainer& sites, const SubstitutionModel& model) throw (Exception);
 
448
    
 
449
};
 
450
 
 
451
} //end of namespace bpp.
 
452
 
 
453
#endif //_DRASDRHOMOGENEOUSTREELIKELIHOODDATA_H_
 
454