← Back to branch summary

~ubuntu-branches/ubuntu/saucy/libbpp-phyl/saucy

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Model/MixedSubstitutionModelSet.cpp

  • Committer: Package Import Robot
  • Author(s): Julien Dutheil
  • Date: 2013-02-09 16:00:00 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20130209160000-5v65ba68z8032nzj
Tags: 2.0.3-1
* Reorganized model hierarchy
* New pairwise models
* Several bugs fixed

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Bpp/Phyl/Model/MixedSubstitutionModelSet.cpp
 
1
//
 
2
// File: MixedSubstitutionModelSet.cpp
 
3
// Created by: Laurent Guéguen
 
4
// Created on: mercredi 25 mai 2011, à 22h 12
 
5
//
 
6
 
 
7
/*
 
8
   Copyright or <A9> or Copr. CNRS, (November 16, 2004)
 
9
 
 
10
   This software is a computer program whose purpose is to provide classes
 
11
   for phylogenetic data analysis.
 
12
 
 
13
   This software is governed by the CeCILL  license under French law and
 
14
   abiding by the rules of distribution of free software.  You can  use,
 
15
   modify and/ or redistribute the software under the terms of the CeCILL
 
16
   license as circulated by CEA, CNRS and INRIA at the following URL
 
17
   "http://www.cecill.info".
 
18
 
 
19
   As a counterpart to the access to the source code and  rights to copy,
 
20
   modify and redistribute granted by the license, users are provided only
 
21
   with a limited warranty  and the software's author,  the holder of the
 
22
   economic rights,  and the successive licensors  have only  limited
 
23
   liability.
 
24
 
 
25
   In this respect, the user's attention is drawn to the risks associated
 
26
   with loading,  using,  modifying and/or developing or reproducing the
 
27
   software by the user in light of its specific status of free software,
 
28
   that may mean  that it is complicated to manipulate,  and  that  also
 
29
   therefore means  that it is reserved for developers  and  experienced
 
30
   professionals having in-depth computer knowledge. Users are therefore
 
31
   encouraged to load and test the software's suitability as regards their
 
32
   requirements in conditions enabling the security of their systems and/or
 
33
   data to be ensured and,  more generally, to use and operate it in the
 
34
   same conditions as regards security.
 
35
 
 
36
   The fact that you are presently reading this means that you have had
 
37
   knowledge of the CeCILL license and that you accept its terms.
 
38
 */
 
39
 
 
40
#include "MixedSubstitutionModelSet.h"
 
41
#include "MixedSubstitutionModel.h"
 
42
#include "MixtureOfASubstitutionModel.h"
 
43
 
 
44
using namespace bpp;
 
45
using namespace std;
 
46
 
 
47
MixedSubstitutionModelSet::MixedSubstitutionModelSet(const MixedSubstitutionModelSet& set) :
 
48
  SubstitutionModelSet(set), vpHyperNodes_()
 
49
{
 
50
  for (unsigned int i=0;i<set.vpHyperNodes_.size();i++)
 
51
    vpHyperNodes_.push_back(new HyperNode(*set.vpHyperNodes_[i]));  
 
52
}
 
53
 
 
54
MixedSubstitutionModelSet::MixedSubstitutionModelSet(const SubstitutionModelSet& set) :
 
55
  SubstitutionModelSet(set), vpHyperNodes_()
 
56
{
 
57
}
 
58
 
 
59
MixedSubstitutionModelSet::~MixedSubstitutionModelSet() 
 
60
{
 
61
  for (unsigned int i=0;i<vpHyperNodes_.size();i++)
 
62
    delete vpHyperNodes_[i];
 
63
}
 
64
 
 
65
void MixedSubstitutionModelSet::clear() 
 
66
{
 
67
  SubstitutionModelSet::clear();
 
68
  for (unsigned int i=0;i<vpHyperNodes_.size();i++)
 
69
    delete vpHyperNodes_[i];
 
70
}
 
71
 
 
72
MixedSubstitutionModelSet& MixedSubstitutionModelSet::operator=(const MixedSubstitutionModelSet& set)
 
73
{
 
74
  SubstitutionModelSet::operator=(set);
 
75
  for (unsigned int i=0;i<vpHyperNodes_.size();i++)
 
76
    if (vpHyperNodes_[i]!=NULL)
 
77
      delete vpHyperNodes_[i];
 
78
  vpHyperNodes_.clear();
 
79
  
 
80
  for (unsigned int i=0;i<set.vpHyperNodes_.size();i++)
 
81
    vpHyperNodes_.push_back(new HyperNode(*set.vpHyperNodes_[i]));
 
82
  
 
83
  return *this;
 
84
}
 
85
 
 
86
void MixedSubstitutionModelSet::addEmptyHyperNode()
 
87
{
 
88
  vpHyperNodes_.push_back(new HyperNode(this));
 
89
}
 
90
 
 
91
void MixedSubstitutionModelSet::addHyperNode(const HyperNode& hn)
 
92
{
 
93
  vpHyperNodes_.push_back(new HyperNode(hn));
 
94
}
 
95
 
 
96
bool MixedSubstitutionModelSet::complete()
 
97
{
 
98
  MixedSubstitutionModelSet::HyperNode nhn(this);
 
99
  unsigned int i;
 
100
  for ( i=0;i<vpHyperNodes_.size();i++)
 
101
    nhn+=*vpHyperNodes_[i];
 
102
 
 
103
  unsigned int nbm=getNumberOfModels();
 
104
  for (i=0;i<nbm;i++){
 
105
    const MixedSubstitutionModel* pSM=dynamic_cast<const MixedSubstitutionModel*>(getModel(i));
 
106
    if (pSM!=NULL){
 
107
      if (nhn.getNode(i).size()!=pSM->getNumberOfModels())
 
108
        break;
 
109
    }
 
110
  }
 
111
 
 
112
  if (i==nbm)
 
113
    return false;
 
114
  
 
115
  addEmptyHyperNode();
 
116
  for ( i=0;i<nbm;i++){
 
117
    const MixedSubstitutionModel* pSM=dynamic_cast<const MixedSubstitutionModel*>(getModel(i));
 
118
    if (pSM!=NULL){
 
119
      const MixedSubstitutionModelSet::HyperNode::Node& nd=nhn.getNode(i);
 
120
      int snd=(int)nd.size();
 
121
      int vs=(int)pSM->getNumberOfModels();
 
122
      Vint an;
 
123
 
 
124
      int j(0), k(0);
 
125
      while (j<vs){
 
126
        while ((k<snd) && (nd[k]<j))
 
127
          k++;
 
128
        if ((k>=snd) || (nd[k]>j))
 
129
          an.push_back(j);
 
130
        j++;
 
131
      }
 
132
      addToHyperNode(i,an);
 
133
    }
 
134
  }
 
135
  
 
136
  return true;
 
137
}
 
138
 
 
139
void MixedSubstitutionModelSet::addToHyperNode(int nM, const Vint& vnS, int nH)
 
140
{
 
141
  if (nH>=(int)vpHyperNodes_.size())
 
142
    throw BadIntegerException("Bad HyperNode number",nH);
 
143
  if (nH<0)
 
144
    nH=vpHyperNodes_.size()-1;
 
145
  
 
146
  if (nM<0 || nM>=(int)getNumberOfModels())
 
147
    throw BadIntegerException("Bad Mixed Model number",nM);
 
148
 
 
149
  vpHyperNodes_[nH]->addToModel(nM,vnS);
 
150
}
 
151
 
 
152
bool MixedSubstitutionModelSet::hasExclusivePaths() const
 
153
{
 
154
  HyperNode tthn(this);
 
155
 
 
156
  unsigned int nhn=getNumberOfHyperNodes();
 
157
  for (unsigned int i=0;i<nhn;i++){
 
158
    if (tthn.intersects(getHyperNode(i)))
 
159
      return false;
 
160
    else
 
161
      tthn+=getHyperNode(i);
 
162
    }
 
163
  
 
164
  return true;
 
165
}
 
166
 
 
167
void MixedSubstitutionModelSet::fireParameterChanged(const ParameterList& parameters)
 
168
{
 
169
  SubstitutionModelSet::fireParameterChanged(parameters);
 
170
 
 
171
  // should be restricted only when probability related parameters are changed 
 
172
  computeHyperNodesProbabilities();
 
173
}
 
174
 
 
175
 
 
176
void MixedSubstitutionModelSet::computeHyperNodesProbabilities()
 
177
{
 
178
  unsigned int nbm=getNumberOfModels();
 
179
 
 
180
  // Looking for the first Mixed model
 
181
  
 
182
  unsigned int fmM=0;
 
183
  
 
184
  MixedSubstitutionModel* pfSM=0;
 
185
  for (fmM=0;fmM<nbm;fmM++){
 
186
    pfSM=dynamic_cast<MixedSubstitutionModel*>(getModel(fmM));
 
187
    if (pfSM!=NULL)
 
188
      break;
 
189
  }
 
190
  if (fmM==nbm)
 
191
    return;
 
192
 
 
193
  // Compute the probabilities of the hypernodes from the first mixed
 
194
  // model
 
195
  
 
196
  unsigned int nbh=getNumberOfHyperNodes();
 
197
 
 
198
  
 
199
  for (unsigned int nh=0;nh<nbh;nh++){
 
200
    MixedSubstitutionModelSet::HyperNode& h=getHyperNode(nh);
 
201
    const MixedSubstitutionModelSet::HyperNode::Node& fnd=h.getNode(fmM);
 
202
    
 
203
    double fprob=0;
 
204
    for (unsigned int i=0;i<fnd.size();i++)
 
205
      fprob+=pfSM->getNProbability(fnd[i]);
 
206
    h.setProbability(fprob);
 
207
  }
 
208
 
 
209
  // Sets the new probabilities & rates of the mixmodels
 
210
  
 
211
  for (unsigned int iM=fmM+1;iM<nbm;iM++){
 
212
    pfSM=dynamic_cast<MixedSubstitutionModel*>(getModel(iM));
 
213
    if (pfSM!=NULL){
 
214
      for (unsigned int nh=0;nh<nbh;nh++){
 
215
        MixedSubstitutionModelSet::HyperNode& h=getHyperNode(nh);
 
216
        const MixedSubstitutionModelSet::HyperNode::Node& fnd=h.getNode(iM);
 
217
        double prob=0;
 
218
        for (unsigned int j=0;j<fnd.size();j++)
 
219
          prob+=pfSM->getNProbability(fnd[j]);
 
220
 
 
221
        // sets the real probabilities
 
222
        for (unsigned int j=0;j<fnd.size();j++)
 
223
          pfSM->setNProbability(fnd[j],h.getProbability()*pfSM->getNProbability(fnd[j])/prob);
 
224
      }
 
225
 
 
226
      // normalizes Vrates with the real probabilities
 
227
 
 
228
      pfSM->normalizeVRates();
 
229
 
 
230
      // sets the conditional probabilities
 
231
 
 
232
      for (unsigned int nh=0;nh<nbh;nh++){
 
233
        MixedSubstitutionModelSet::HyperNode& h=getHyperNode(nh);
 
234
        const MixedSubstitutionModelSet::HyperNode::Node& fnd=h.getNode(iM);
 
235
        for (unsigned int j=0;j<fnd.size();j++)
 
236
          pfSM->setNProbability(fnd[j],pfSM->getNProbability(fnd[j])/h.getProbability());
 
237
      }
 
238
    }
 
239
  }  
 
240
}
 
241
 
 
242
double MixedSubstitutionModelSet::getHyperNodeProbability(const HyperNode& hn) const
 
243
{
 
244
  unsigned int nbm=getNumberOfModels();
 
245
  double fprob=1;
 
246
  
 
247
  for (unsigned int fmM=0;fmM<nbm;fmM++){
 
248
    
 
249
    const MixedSubstitutionModelSet::HyperNode::Node& fnd=hn.getNode(fmM);
 
250
    const MixedSubstitutionModel* pfSM=dynamic_cast<const MixedSubstitutionModel*>(getModel(fmM));
 
251
    if (pfSM!=NULL){
 
252
      double x=0;
 
253
 
 
254
      for (unsigned int i=0;i<fnd.size();i++)
 
255
        x+=pfSM->getNProbability(fnd[i]);
 
256
 
 
257
      fprob*=x;
 
258
    }
 
259
  }
 
260
 
 
261
  return fprob;
 
262
}
 
263
 
 
264
/**********************************************************/
 
265
/*************** HYPERNODE ********************************/
 
266
/***********************************************************/
 
267
 
 
268
 
 
269
MixedSubstitutionModelSet::HyperNode::HyperNode(const MixedSubstitutionModelSet *pMSMS): vNumbers_(pMSMS->getNumberOfModels()), vUnused_(), proba_(1)
 
270
{
 
271
  for (unsigned int i=0;i<pMSMS->getNumberOfModels();i++){
 
272
    const MixedSubstitutionModel* pSM=dynamic_cast<const MixedSubstitutionModel*>(pMSMS->getModel(i));
 
273
    if (pSM==NULL)
 
274
      vUnused_.push_back(i);
 
275
  }
 
276
}
 
277
 
 
278
MixedSubstitutionModelSet::HyperNode::HyperNode(const HyperNode& hn): vNumbers_(hn.vNumbers_), vUnused_(hn.vUnused_), proba_(hn.proba_)
 
279
{}
 
280
 
 
281
 
 
282
MixedSubstitutionModelSet::HyperNode& MixedSubstitutionModelSet::HyperNode::operator=(const MixedSubstitutionModelSet::HyperNode& hn)
 
283
{
 
284
  vNumbers_.clear();
 
285
  vNumbers_.resize(hn.vNumbers_.size());
 
286
  for (unsigned int i=0;i<hn.vNumbers_.size();i++)
 
287
    vNumbers_[i]=hn.vNumbers_[i];
 
288
  vUnused_.clear();
 
289
  vUnused_.resize(hn.vUnused_.size());
 
290
  for (unsigned int i=0;i<hn.vUnused_.size();i++)
 
291
    vUnused_[i]=hn.vUnused_[i];
 
292
 
 
293
  proba_=hn.proba_;
 
294
  
 
295
  return *this;
 
296
}
 
297
 
 
298
void MixedSubstitutionModelSet::HyperNode::addToModel(int nM, const Vint& vnS)
 
299
{
 
300
  if ((nM<0) || (nM>=(int)vNumbers_.size()))
 
301
    throw BadIntegerException("Bad Mixed model Number",nM);
 
302
      
 
303
  vNumbers_[nM].insertN(vnS);
 
304
}
 
305
 
 
306
void MixedSubstitutionModelSet::HyperNode::setModel(int nM, const Vint& vnS)
 
307
{
 
308
  if ((nM<0) || (nM>=(int)vNumbers_.size()))
 
309
    throw BadIntegerException("Bad Mixed model Number",nM);
 
310
      
 
311
  vNumbers_[nM]=vnS;
 
312
}
 
313
 
 
314
bool MixedSubstitutionModelSet::HyperNode::isComplete() const
 
315
{
 
316
  int k;
 
317
  int vUs=(int)vUnused_.size();
 
318
  for (int i=0;i<(int)vNumbers_.size();i++){
 
319
    for (k=0;k<vUs;k++)
 
320
      if (vUnused_[k]==i)
 
321
        break;
 
322
    if ((k==vUs) && vNumbers_[i].size()==0)
 
323
      return false;
 
324
  }
 
325
  return true;
 
326
}
 
327
 
 
328
bool MixedSubstitutionModelSet::HyperNode::operator<=(const HyperNode& hn) const
 
329
{
 
330
  for (unsigned int i=0;i<vNumbers_.size();i++){
 
331
    if (!( vNumbers_[i]<=hn.vNumbers_[i]))
 
332
      return false;
 
333
  }
 
334
 
 
335
  return true;
 
336
}
 
337
 
 
338
bool MixedSubstitutionModelSet::HyperNode::intersects(const HyperNode& hn) const
 
339
{
 
340
  for (unsigned int i=0;i<vNumbers_.size();i++)
 
341
    if (vNumbers_[i].intersects(hn.vNumbers_[i]))
 
342
      return true;
 
343
  
 
344
  return false;
 
345
}
 
346
 
 
347
bool MixedSubstitutionModelSet::HyperNode::operator>=(const HyperNode& hn) const
 
348
{
 
349
  return hn>=*this;
 
350
}
 
351
 
 
352
MixedSubstitutionModelSet::HyperNode& MixedSubstitutionModelSet::HyperNode::operator+=(const HyperNode& hn)
 
353
{
 
354
  for (unsigned int i=0;i<vNumbers_.size();i++)
 
355
    vNumbers_[i]+=hn.vNumbers_[i];
 
356
  
 
357
  return *this;
 
358
}
 
359
 
 
360
/**********************************************************/
 
361
/******************** NODE ********************************/
 
362
/***********************************************************/
 
363
 
 
364
void MixedSubstitutionModelSet::HyperNode::Node::insertN(const Vint& vn)
 
365
{
 
366
  vector<int>::iterator it;
 
367
  vector<int>::const_iterator it2;
 
368
  
 
369
  for (it2=vn.begin();it2!=vn.end();it2++){
 
370
    for (it=vNumb_.begin();it!=vNumb_.end();it++)
 
371
      if (*it>=*it2)
 
372
        break;
 
373
    if (it==vNumb_.end())
 
374
        vNumb_.push_back(*it2);
 
375
    else if (*it!=*it2)
 
376
      vNumb_.insert(it,*it2);
 
377
  }
 
378
}
 
379
 
 
380
MixedSubstitutionModelSet::HyperNode::Node& MixedSubstitutionModelSet::HyperNode::Node::operator+=(const Node& n)
 
381
{
 
382
  insertN(n.vNumb_);
 
383
 
 
384
  return *this;
 
385
}
 
386
 
 
387
bool MixedSubstitutionModelSet::HyperNode::Node::operator<=(const Node& n) const
 
388
{
 
389
  vector<int>::const_iterator it(vNumb_.begin());
 
390
  vector<int>::const_iterator it2(n.vNumb_.begin());
 
391
  
 
392
  for (;it!=vNumb_.end();it++){
 
393
    while (it2!=n.vNumb_.end()  && (*it2<*it))
 
394
      it2++;
 
395
    if (it2==n.vNumb_.end() || (*it2>*it))
 
396
      return false;
 
397
    it++;
 
398
  }
 
399
  return true;
 
400
}
 
401
 
 
402
bool MixedSubstitutionModelSet::HyperNode::Node::operator>=(const Node& n) const
 
403
{
 
404
  return n<=*this;
 
405
}
 
406
 
 
407
bool MixedSubstitutionModelSet::HyperNode::Node::intersects(const Node& n) const
 
408
{
 
409
  vector<int>::const_iterator it(vNumb_.begin());
 
410
  vector<int>::const_iterator it2(n.vNumb_.begin());
 
411
  
 
412
  for (;it!=vNumb_.end();it++){
 
413
    while (it2!=n.vNumb_.end()  && (*it2<*it))
 
414
      it2++;
 
415
 
 
416
    if (it2==n.vNumb_.end())
 
417
      return false;
 
418
    if (*it2==*it)
 
419
      return true;
 
420
    it++;
 
421
  }
 
422
  return false;
 
423
}
 
424