← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Distance/PGMA.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/Distance/PGMA.cpp
 
1
//
 
2
// File: PGMA.cpp
 
3
// Created by: Julien Dutheil
 
4
// Created on: Mon jul 11 11:41 2005
 
5
//
 
6
 
 
7
/*
 
8
Copyright or © 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 "PGMA.h"
 
41
#include "../NodeTemplate.h"
 
42
#include "../Tree.h"
 
43
#include "../TreeTemplate.h"
 
44
#include "../TreeTemplateTools.h"
 
45
 
 
46
using namespace bpp;
 
47
 
 
48
// From the STL:
 
49
#include <cmath>
 
50
#include <iostream>
 
51
 
 
52
using namespace std;
 
53
 
 
54
TreeTemplate<Node>* PGMA::getTree() const
 
55
{
 
56
        Node* root = TreeTemplateTools::cloneSubtree<Node>(*dynamic_cast<TreeTemplate<NodeTemplate<PGMAInfos> >*>(tree_)->getRootNode());
 
57
        return new TreeTemplate<Node>(root);
 
58
}
 
59
        
 
60
vector<unsigned int> PGMA::getBestPair() throw (Exception)
 
61
{
 
62
        vector<unsigned int> bestPair(2);
 
63
        double distMin = -std::log(0.);
 
64
        for(map<unsigned int, Node *>::iterator i = currentNodes_.begin(); i != currentNodes_.end(); i++) {
 
65
                unsigned int id = i -> first;
 
66
                map<unsigned int, Node *>::iterator j = i;
 
67
                j++;
 
68
                for(; j != currentNodes_.end(); j++) {
 
69
                        unsigned int jd = j -> first;
 
70
                        double dist = matrix_(id, jd);
 
71
                        if(dist < distMin) {
 
72
                                distMin = dist;
 
73
                                bestPair[0] = id;
 
74
                                bestPair[1] = jd;
 
75
                        }
 
76
                }
 
77
        }
 
78
 
 
79
        if(distMin == -std::log(0.)) {
 
80
                throw Exception("Unexpected error: no minimum found in the distance matrix.");
 
81
        }
 
82
 
 
83
        return bestPair;        
 
84
}
 
85
 
 
86
vector<double> PGMA::computeBranchLengthsForPair(const vector<unsigned int> & pair)
 
87
{
 
88
        vector<double> d(2);
 
89
        double dist = matrix_(pair[0], pair[1]) / 2.;
 
90
        d[0] = dist - dynamic_cast<NodeTemplate<PGMAInfos> *>(currentNodes_[pair[0]]) -> getInfos().time; 
 
91
        d[1] = dist - dynamic_cast<NodeTemplate<PGMAInfos> *>(currentNodes_[pair[1]]) -> getInfos().time; 
 
92
        return d;
 
93
}
 
94
 
 
95
double PGMA::computeDistancesFromPair(const vector<unsigned int> & pair, const vector<double> & branchLengths, unsigned int pos)
 
96
{
 
97
        double w1, w2;
 
98
        if(weighted_) {
 
99
                w1 = 1;
 
100
                w2 = 1;
 
101
        } else {
 
102
                w1 = dynamic_cast<NodeTemplate<PGMAInfos> *>(currentNodes_[pair[0]]) -> getInfos().numberOfLeaves;
 
103
                w2 = dynamic_cast<NodeTemplate<PGMAInfos> *>(currentNodes_[pair[1]]) -> getInfos().numberOfLeaves;
 
104
        }
 
105
        return (w1 * matrix_(pair[0], pos) + w2 * matrix_(pair[1], pos)) / (w1 + w2); 
 
106
}
 
107
 
 
108
void PGMA::finalStep(int idRoot)
 
109
{
 
110
        NodeTemplate<PGMAInfos>* root = new NodeTemplate<PGMAInfos>(idRoot);
 
111
        map<unsigned int, Node*>::iterator it = currentNodes_.begin();
 
112
        unsigned int i1 = it->first;
 
113
        Node* n1        = it->second;
 
114
        it++;
 
115
        unsigned int i2 = it->first;
 
116
        Node* n2        = it->second;
 
117
        double d = matrix_(i1, i2) / 2;
 
118
        root->addSon(n1);
 
119
        root->addSon(n2);
 
120
        n1->setDistanceToFather(d - dynamic_cast<NodeTemplate<PGMAInfos>*>(n1)->getInfos().time); 
 
121
        n2->setDistanceToFather(d - dynamic_cast<NodeTemplate<PGMAInfos>*>(n2)->getInfos().time); 
 
122
        tree_ = new TreeTemplate<NodeTemplate<PGMAInfos> >(root);
 
123
}
 
124
 
 
125
Node* PGMA::getLeafNode(int id, const string& name)
 
126
{
 
127
        PGMAInfos infos;
 
128
        infos.numberOfLeaves = 1;
 
129
        infos.time = 0.;
 
130
        NodeTemplate<PGMAInfos> * leaf = new NodeTemplate<PGMAInfos>(id, name);
 
131
        leaf -> setInfos(infos);
 
132
        return leaf;
 
133
}
 
134
 
 
135
Node* PGMA::getParentNode(int id, Node* son1, Node* son2)
 
136
{
 
137
        PGMAInfos infos;
 
138
        infos.numberOfLeaves = 
 
139
                dynamic_cast<NodeTemplate<PGMAInfos> *>(son1) -> getInfos().numberOfLeaves
 
140
        + dynamic_cast<NodeTemplate<PGMAInfos> *>(son2) -> getInfos().numberOfLeaves;
 
141
        infos.time = dynamic_cast<NodeTemplate<PGMAInfos> *>(son1) -> getInfos().time + son1 -> getDistanceToFather();
 
142
        Node* parent = new NodeTemplate<PGMAInfos>(id);
 
143
        dynamic_cast<NodeTemplate<PGMAInfos> *>(parent) -> setInfos(infos);
 
144
        parent->addSon(son1);
 
145
        parent->addSon(son2);
 
146
        return parent;
 
147
}
 
148