← Back to branch summary

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

« back to all changes in this revision

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

  • 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/Model/AbstractWordReversibleSubstitutionModel.cpp
 
1
//
 
2
// File: AbstractWordReversibleSubstitutionModel.cpp
 
3
// Created by:  Laurent Gueguen
 
4
// Created on: Jan 2009
 
5
//
 
6
 
 
7
/*
 
8
   Copyright or © or Copr. CNRS, (November 16, 2004)
 
9
   This software is a computer program whose purpose is to provide classes
 
10
   for phylogenetic data analysis.
 
11
 
 
12
   This software is governed by the CeCILL  license under French law and
 
13
   abiding by the rules of distribution of free software.  You can  use,
 
14
   modify and/ or redistribute the software under the terms of the CeCILL
 
15
   license as circulated by CEA, CNRS and INRIA at the following URL
 
16
   "http://www.cecill.info".
 
17
 
 
18
   As a counterpart to the access to the source code and  rights to copy,
 
19
   modify and redistribute granted by the license, users are provided only
 
20
   with a limited warranty  and the software's author,  the holder of the
 
21
   economic rights,  and the successive licensors  have only  limited
 
22
   liability.
 
23
 
 
24
   In this respect, the user's attention is drawn to the risks associated
 
25
   with loading,  using,  modifying and/or developing or reproducing the
 
26
   software by the user in light of its specific status of free software,
 
27
   that may mean  that it is complicated to manipulate,  and  that  also
 
28
   therefore means  that it is reserved for developers  and  experienced
 
29
   professionals having in-depth computer knowledge. Users are therefore
 
30
   encouraged to load and test the software's suitability as regards their
 
31
   requirements in conditions enabling the security of their systems and/or
 
32
   data to be ensured and,  more generally, to use and operate it in the
 
33
   same conditions as regards security.
 
34
 
 
35
   The fact that you are presently reading this means that you have had
 
36
   knowledge of the CeCILL license and that you accept its terms.
 
37
 */
 
38
 
 
39
#include "AbstractWordReversibleSubstitutionModel.h"
 
40
 
 
41
#include <Bpp/Numeric/Matrix/MatrixTools.h>
 
42
#include <Bpp/Numeric/Matrix/EigenValue.h>
 
43
#include <Bpp/Numeric/VectorTools.h>
 
44
 
 
45
// From SeqLib:
 
46
#include <Bpp/Seq/Alphabet/WordAlphabet.h>
 
47
#include <Bpp/Seq/Alphabet/AlphabetTools.h>
 
48
#include <Bpp/Seq/Container/SequenceContainerTools.h>
 
49
 
 
50
using namespace bpp;
 
51
 
 
52
// From the STL:
 
53
#include <cmath>
 
54
 
 
55
using namespace std;
 
56
 
 
57
/******************************************************************************/
 
58
 
 
59
AbstractWordReversibleSubstitutionModel::AbstractWordReversibleSubstitutionModel(
 
60
  const std::vector<SubstitutionModel*>& modelVector,
 
61
  const std::string& st) :
 
62
  AbstractReversibleSubstitutionModel(AbstractWordReversibleSubstitutionModel::extractAlph(modelVector), st),
 
63
  new_alphabet_ (true),
 
64
  VSubMod_      (),
 
65
  VnestedPrefix_(),
 
66
  Vrate_         (modelVector.size())
 
67
{
 
68
  enableEigenDecomposition(false);
 
69
  unsigned int i, j;
 
70
  unsigned int n = modelVector.size();
 
71
 
 
72
  // test whether two models are identical
 
73
 
 
74
  bool flag = 0;
 
75
  i = 0;
 
76
  j = 1;
 
77
  while (!flag && i < (n - 1))
 
78
  {
 
79
    if (modelVector[i] == modelVector[j])
 
80
      flag = 1;
 
81
    else
 
82
    {
 
83
      j++;
 
84
      if (j == n)
 
85
      {
 
86
        i++;
 
87
        j = i + 1;
 
88
      }
 
89
    }
 
90
  }
 
91
 
 
92
  if (!flag)
 
93
  {
 
94
    for (i = 0; i < n; i++)
 
95
    {
 
96
      VSubMod_.push_back(modelVector[i]);
 
97
      VnestedPrefix_.push_back(modelVector[i]->getNamespace());
 
98
      VSubMod_[i]->setNamespace(st + TextTools::toString(i+1) + "_" + VnestedPrefix_[i]);
 
99
      addParameters_(VSubMod_[i]->getParameters());
 
100
    }
 
101
  }
 
102
  else
 
103
  {
 
104
   string t = "";
 
105
    for (i = 0; i < n; i++)
 
106
    {
 
107
      VSubMod_.push_back(modelVector[0]);
 
108
      VnestedPrefix_.push_back(modelVector[0]->getNamespace());
 
109
      t += TextTools::toString(i+1);
 
110
    }
 
111
    VSubMod_[0]->setNamespace(st + t + "_" + VnestedPrefix_[0]);
 
112
    addParameters_(VSubMod_[0]->getParameters());
 
113
  }
 
114
 
 
115
  for (i = 0; i < n; i++)
 
116
  {
 
117
    Vrate_[i] = 1.0 / n;
 
118
  }
 
119
}
 
120
 
 
121
AbstractWordReversibleSubstitutionModel::AbstractWordReversibleSubstitutionModel(
 
122
  const Alphabet* alph,
 
123
  const std::string& st) :
 
124
  AbstractReversibleSubstitutionModel(alph, st),
 
125
  new_alphabet_ (false),
 
126
  VSubMod_      (),
 
127
  VnestedPrefix_(),
 
128
  Vrate_         (0)
 
129
{
 
130
  enableEigenDecomposition(false);
 
131
}
 
132
 
 
133
AbstractWordReversibleSubstitutionModel::AbstractWordReversibleSubstitutionModel(
 
134
  SubstitutionModel* pmodel,
 
135
  unsigned int num,
 
136
  const std::string& st) :
 
137
  AbstractReversibleSubstitutionModel(new WordAlphabet(pmodel->getAlphabet(), num),st),
 
138
  new_alphabet_ (true),
 
139
  VSubMod_      (),
 
140
  VnestedPrefix_(),
 
141
  Vrate_         (num)
 
142
{
 
143
  enableEigenDecomposition(false);
 
144
  unsigned int i;
 
145
 
 
146
  string t = "";
 
147
  for (i = 0; i < num; i++)
 
148
  {
 
149
   VSubMod_.push_back(pmodel);
 
150
   VnestedPrefix_.push_back(pmodel->getNamespace());
 
151
    Vrate_[i] = 1.0 / num;
 
152
    t += TextTools::toString(i+1);
 
153
  }
 
154
 
 
155
  pmodel->setNamespace(st + t + "_" + VnestedPrefix_[0]);
 
156
  addParameters_(pmodel->getParameters());
 
157
}
 
158
 
 
159
AbstractWordReversibleSubstitutionModel::AbstractWordReversibleSubstitutionModel(
 
160
  const AbstractWordReversibleSubstitutionModel& wrsm) :
 
161
  AbstractReversibleSubstitutionModel(wrsm),
 
162
  new_alphabet_ (wrsm.new_alphabet_),
 
163
  VSubMod_      (),
 
164
  VnestedPrefix_(wrsm.VnestedPrefix_),
 
165
  Vrate_         (wrsm.Vrate_)
 
166
{
 
167
   unsigned int i;
 
168
   unsigned int num = wrsm.VSubMod_.size();
 
169
 
 
170
  if (wrsm.new_alphabet_)
 
171
    alphabet_ = new WordAlphabet(*(dynamic_cast<const WordAlphabet*>(wrsm.getAlphabet())));
 
172
 
 
173
  SubstitutionModel* pSM = 0;
 
174
  if ((num > 1) & (wrsm.VSubMod_[0] == wrsm.VSubMod_[1]))
 
175
    pSM = wrsm.VSubMod_[0]->clone();
 
176
 
 
177
  for (i = 0; i < num; i++)
 
178
  {
 
179
   VSubMod_.push_back(pSM ? pSM : wrsm.VSubMod_[i]->clone());
 
180
  }
 
181
}
 
182
 
 
183
AbstractWordReversibleSubstitutionModel& AbstractWordReversibleSubstitutionModel::operator=(
 
184
  const AbstractWordReversibleSubstitutionModel& wrsm)
 
185
{
 
186
  AbstractReversibleSubstitutionModel::operator=(wrsm);
 
187
  new_alphabet_  = wrsm.new_alphabet_;
 
188
  VnestedPrefix_ = wrsm.VnestedPrefix_;
 
189
  Vrate_          = wrsm.Vrate_;
 
190
 
 
191
  unsigned int i;
 
192
  unsigned int num = wrsm.VSubMod_.size();
 
193
 
 
194
  if (wrsm.new_alphabet_)
 
195
    alphabet_ = new WordAlphabet(*(dynamic_cast<const WordAlphabet*>(wrsm.getAlphabet())));
 
196
 
 
197
  SubstitutionModel* pSM = 0;
 
198
  if ((num > 1) & (wrsm.VSubMod_[0] == wrsm.VSubMod_[1]))
 
199
    pSM = wrsm.VSubMod_[0]->clone();
 
200
 
 
201
  for (i = 0; i < num; i++)
 
202
  {
 
203
    VSubMod_[i] =  (pSM ? pSM : wrsm.VSubMod_[i]->clone());
 
204
  }
 
205
 
 
206
  return *this;
 
207
}
 
208
 
 
209
AbstractWordReversibleSubstitutionModel::~AbstractWordReversibleSubstitutionModel()
 
210
{
 
211
  if ((VSubMod_.size() > 1) && (VSubMod_[0] == VSubMod_[1]))
 
212
  {
 
213
    if (VSubMod_[0])
 
214
      delete VSubMod_[0];
 
215
  }
 
216
  else
 
217
    for (size_t i = 0; i < VSubMod_.size(); i++)
 
218
    {
 
219
      if (VSubMod_[i])
 
220
        delete VSubMod_[i];
 
221
    }
 
222
  if (new_alphabet_)
 
223
    delete alphabet_;
 
224
}
 
225
 
 
226
unsigned int AbstractWordReversibleSubstitutionModel::getNumberOfStates() const
 
227
{
 
228
  return getAlphabet()->getSize();
 
229
}
 
230
 
 
231
Alphabet* AbstractWordReversibleSubstitutionModel::extractAlph(const vector<SubstitutionModel*>& modelVector)
 
232
{
 
233
   unsigned int i;
 
234
 
 
235
   vector<const Alphabet*> vAlph;
 
236
 
 
237
  for (i = 0; i < modelVector.size(); i++)
 
238
  {
 
239
   vAlph.push_back(modelVector[i]->getAlphabet());
 
240
  }
 
241
 
 
242
  return new WordAlphabet(vAlph);
 
243
}
 
244
 
 
245
void AbstractWordReversibleSubstitutionModel::setNamespace(const std::string& prefix)
 
246
{
 
247
   AbstractReversibleSubstitutionModel::setNamespace(prefix);
 
248
 
 
249
  if (VSubMod_.size() < 2 || VSubMod_[0] == VSubMod_[1])
 
250
  {
 
251
   string t = "";
 
252
    for (unsigned int i = 0; i < VSubMod_.size(); i++)
 
253
    {
 
254
      t += TextTools::toString(i+1);
 
255
    }
 
256
    VSubMod_[0]->setNamespace(prefix + t + "_" + VnestedPrefix_[0]);
 
257
  }
 
258
  else
 
259
  {
 
260
    for (unsigned int i = 0; i < VSubMod_.size(); i++)
 
261
    {
 
262
      VSubMod_[i]->setNamespace(prefix + TextTools::toString(i+1) + "_" + VnestedPrefix_[i]);
 
263
    }
 
264
  }
 
265
}
 
266
 
 
267
/******************************************************************************/
 
268
 
 
269
void AbstractWordReversibleSubstitutionModel::updateMatrices()
 
270
{
 
271
  //First we update position specific models. This need to be done here and not
 
272
  //in fireParameterChanged, has some parameter aliases might have been defined
 
273
  //and need to be resolved first.
 
274
  if (VSubMod_.size() < 2 || VSubMod_[0] == VSubMod_[1])
 
275
    VSubMod_[0]->matchParametersValues(getParameters());
 
276
  else
 
277
    for (unsigned int i = 0; i < VSubMod_.size(); i++)
 
278
    {
 
279
      VSubMod_[i]->matchParametersValues(getParameters());
 
280
    }
 
281
 
 
282
  unsigned int nbmod = VSubMod_.size();
 
283
  unsigned int salph = getNumberOfStates();
 
284
 
 
285
  // Generator
 
286
 
 
287
  if (enableEigenDecomposition())
 
288
  {
 
289
    unsigned int i, j, n, l, k, m;
 
290
 
 
291
    vector<unsigned int> vsize;
 
292
 
 
293
    for (k = 0; k < nbmod; k++)
 
294
    {
 
295
      vsize.push_back(VSubMod_[k]->getNumberOfStates());
 
296
    }
 
297
 
 
298
    RowMatrix<double> gk, exch;
 
299
 
 
300
    m = 1;
 
301
 
 
302
    for (k = nbmod; k > 0; k--)
 
303
    {
 
304
      gk = VSubMod_[k - 1]->getGenerator();
 
305
      exch = (dynamic_cast<AbstractReversibleSubstitutionModel*>(VSubMod_[k - 1]))->getExchangeabilityMatrix();
 
306
      for (i = 0; i < vsize[k - 1]; i++)
 
307
      {
 
308
        for (j = 0; j < vsize[k - 1]; j++)
 
309
        {
 
310
          if (i != j)
 
311
          {
 
312
            n = 0;
 
313
            while (n < salph)
 
314
            { // loop on prefix
 
315
              for (l = 0; l < m; l++)
 
316
              { // loop on suffix
 
317
                generator_(n + i * m + l, n + j * m + l) = gk(i,j) * Vrate_[k - 1];
 
318
                exchangeability_(n + i * m + l, n + j * m + l) = exch(i,j) * Vrate_[k - 1];
 
319
              }
 
320
              n += m * vsize[k - 1];
 
321
            }
 
322
          }
 
323
        }
 
324
      }
 
325
      m *= vsize[k - 1];
 
326
    }
 
327
  }
 
328
 
 
329
  // modification of generator_ and freq_
 
330
 
 
331
  completeMatrices();
 
332
 
 
333
  // at that point generator_ and freq_ are done for models without
 
334
  // enableEigenDecomposition
 
335
 
 
336
  // Eigen values:
 
337
 
 
338
  if (enableEigenDecomposition())
 
339
  {
 
340
    unsigned int i, j;
 
341
    double x;
 
342
 
 
343
    unsigned int nbStop;
 
344
    Vdouble vi;
 
345
 
 
346
    for (i = 0; i < salph; i++)
 
347
    {
 
348
      x = 0;
 
349
      for (j = 0; j < salph; j++)
 
350
      {
 
351
        if (j != i)
 
352
          x += generator_(i, j);
 
353
      }
 
354
      generator_(i, i) = -x;
 
355
    }
 
356
 
 
357
    //02/03/10 Julien: this should be avoided, we have to find a way to avoid particular cases like this...
 
358
    if (AlphabetTools::isCodonAlphabet(getAlphabet()))
 
359
    {
 
360
      int gi = 0, gj = 0;
 
361
 
 
362
      const CodonAlphabet* pca = dynamic_cast<const CodonAlphabet*>(getAlphabet());
 
363
 
 
364
      RowMatrix<double> gk;
 
365
 
 
366
      nbStop = pca->numberOfStopCodons();
 
367
      gk.resize(salph - nbStop, salph - nbStop);
 
368
      for (i = 0; i < salph; i++)
 
369
      {
 
370
        if (!pca->isStop(i))
 
371
        {
 
372
          gj = 0;
 
373
          for (j = 0; j < salph; j++)
 
374
          {
 
375
            if (!pca->isStop(j))
 
376
            {
 
377
              gk(i - gi, j - gj) = generator_(i, j);
 
378
            }
 
379
            else
 
380
              gj++;
 
381
          }
 
382
        }
 
383
        else
 
384
          gi++;
 
385
      }
 
386
 
 
387
      EigenValue<double> ev(gk);
 
388
      eigenValues_ = ev.getRealEigenValues();
 
389
      vi = ev.getImagEigenValues();
 
390
 
 
391
      for (i = 0; i < nbStop; i++)
 
392
      {
 
393
        eigenValues_.push_back(0);
 
394
      }
 
395
 
 
396
      RowMatrix<double> rev = ev.getV();
 
397
      rightEigenVectors_.resize(salph, salph);
 
398
      gi = 0;
 
399
      for (i = 0; i < salph; i++)
 
400
      {
 
401
        if (pca->isStop(i))
 
402
        {
 
403
          gi++;
 
404
          for (j = 0; j < salph; j++)
 
405
          {
 
406
            rightEigenVectors_(i,j) = 0;
 
407
          }
 
408
          rightEigenVectors_(i,salph - nbStop + gi - 1) = 1;
 
409
        }
 
410
        else
 
411
        {
 
412
          for (j = 0; j < salph - nbStop; j++)
 
413
          {
 
414
            rightEigenVectors_(i, j) = rev(i - gi,j);
 
415
          }
 
416
          for (j = salph - nbStop; j < salph; j++)
 
417
          {
 
418
            rightEigenVectors_(i,j) = 0;
 
419
          }
 
420
        }
 
421
      }
 
422
    }
 
423
    else
 
424
    {
 
425
      EigenValue<double> ev(generator_);
 
426
      eigenValues_ = ev.getRealEigenValues();
 
427
      vi = ev.getImagEigenValues();
 
428
      rightEigenVectors_ = ev.getV();
 
429
      nbStop = 0;
 
430
    }
 
431
 
 
432
    MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
 
433
 
 
434
    // looking for the 0 eigenvector for which the eigen vector
 
435
    // elements are of the same sign.
 
436
    // There is a threshold value in case of approximations
 
437
 
 
438
    double seuil=0.00000001;
 
439
    bool flag=true;
 
440
    unsigned int nulleigen;
 
441
    while(flag){
 
442
      seuil*=10;
 
443
      nulleigen=0;
 
444
      int signe=0;
 
445
 
 
446
      while (flag && (nulleigen < salph - nbStop)) {
 
447
          if (abs(eigenValues_[nulleigen]) < 0.0000001 && abs(vi[nulleigen]) < 0.0000001){
 
448
            signe=0;
 
449
            i=0;
 
450
            while (signe==0 && i< salph){
 
451
              x=leftEigenVectors_(nulleigen, i);
 
452
              signe=x>seuil?1:x< -seuil?-1:0;
 
453
              i++;
 
454
            }
 
455
            if (signe==0)
 
456
              nulleigen++;
 
457
            else {
 
458
              while (i<salph){
 
459
                x=leftEigenVectors_(nulleigen, i);
 
460
                if ((signe==-1 && x>seuil) || (signe==1 && x<-seuil))
 
461
                  break;
 
462
                
 
463
                i++;
 
464
              }
 
465
              if (i<salph-nbStop)
 
466
                nulleigen++;
 
467
              else
 
468
                flag=false;
 
469
            }
 
470
          }
 
471
          else
 
472
            nulleigen++;
 
473
      }
 
474
    }
 
475
    
 
476
    for (i = 0; i < salph; i++)
 
477
    {
 
478
      freq_[i] = leftEigenVectors_(nulleigen, i);
 
479
    }
 
480
 
 
481
    x = 0;
 
482
    for (i = 0; i < salph; i++)
 
483
    {
 
484
      x += freq_[i];
 
485
    }
 
486
 
 
487
    for (i = 0; i < salph; i++)
 
488
    {
 
489
      freq_[i] /= x;
 
490
    }
 
491
 
 
492
    // normalization
 
493
 
 
494
    x = 0;
 
495
    for (i = 0; i < salph; i++)
 
496
    {
 
497
      x += freq_[i] * generator_(i, i);
 
498
    }
 
499
 
 
500
    MatrixTools::scale(generator_, -1. / x);
 
501
 
 
502
    for (i = 0; i < salph; i++)
 
503
    {
 
504
      eigenValues_[i] /= -x;
 
505
    }
 
506
  }
 
507
}
 
508
 
 
509
 
 
510
void AbstractWordReversibleSubstitutionModel::setFreq(std::map<int, double>& freqs)
 
511
{
 
512
  map<int, double> tmpFreq;
 
513
  unsigned int nbmod = VSubMod_.size();
 
514
 
 
515
  unsigned int i, j, s, k, d, size;
 
516
  d = size = getNumberOfStates();
 
517
 
 
518
  if (VSubMod_.size() < 2 || VSubMod_[0] == VSubMod_[1]){
 
519
    s = VSubMod_[0]->getAlphabet()->getSize();
 
520
    for (j = 0; j < s; j++)
 
521
        tmpFreq[j] = 0;
 
522
 
 
523
    for (i = 0; i < nbmod; i++){
 
524
      d /= s;
 
525
      for (k = 0; k < size; k++)
 
526
        tmpFreq[(k / d) % s] += freqs[k];
 
527
    }
 
528
    
 
529
    for (k=0;k<s;k++)
 
530
      tmpFreq[k]/=nbmod;
 
531
    
 
532
    VSubMod_[0]->setFreq(tmpFreq);
 
533
    matchParametersValues(VSubMod_[0]->getParameters());
 
534
  }
 
535
  else
 
536
    for (i = 0; i < nbmod; i++){
 
537
        tmpFreq.clear();
 
538
        s = VSubMod_[i]->getAlphabet()->getSize();
 
539
        d /= s;
 
540
        for (j = 0; j < s; j++)
 
541
          {
 
542
            tmpFreq[j] = 0;
 
543
          }
 
544
        for (k = 0; k < size; k++)
 
545
          {
 
546
            tmpFreq[(k / d) % s] += freqs[k];
 
547
          }
 
548
        VSubMod_[i]->setFreq(tmpFreq);
 
549
        matchParametersValues(VSubMod_[i]->getParameters());
 
550
      }
 
551
}