← Back to branch summary

~ubuntu-branches/ubuntu/raring/libbpp-phyl/raring

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/Distance/NeighborJoining.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/NeighborJoining.cpp
 
1
//
 
2
// File: NeighborJoining.cpp
 
3
// Created by: Julien Dutheil
 
4
//             Vincent Ranwez
 
5
// Created on: Thu jun 23 10:39 2005
 
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
#include "NeighborJoining.h"
 
42
#include "../Tree.h"
 
43
 
 
44
using namespace bpp;
 
45
 
 
46
#include <cmath>
 
47
#include <iostream>
 
48
 
 
49
using namespace std;
 
50
 
 
51
std::vector<unsigned int> NeighborJoining::getBestPair() throw (Exception)
 
52
{
 
53
  for (std::map<unsigned int, Node *>::iterator i = currentNodes_.begin(); i != currentNodes_.end(); i++) {
 
54
    unsigned int id = i -> first;
 
55
    sumDist_[id] = 0;
 
56
    for(map<unsigned int, Node *>::iterator j = currentNodes_.begin(); j != currentNodes_.end(); j++) {
 
57
      unsigned int jd = j -> first;
 
58
      sumDist_[id] += matrix_(id, jd);
 
59
    }
 
60
  }
 
61
 
 
62
  vector<unsigned int> bestPair(2);
 
63
  double critMax = 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 crit = sumDist_[id] + sumDist_[jd] - (currentNodes_.size() - 2) * matrix_(id, jd);
 
71
      //cout << "\t" << id << "\t" << jd << "\t" << crit << endl;
 
72
      if(crit > critMax) {
 
73
        critMax = crit;
 
74
        bestPair[0] = id;
 
75
        bestPair[1] = jd;
 
76
      }
 
77
    }
 
78
  }
 
79
 
 
80
  if(critMax == std::log(0.)) {
 
81
    throw Exception("Unexpected error: no maximum criterium found.");
 
82
  }
 
83
  return bestPair;  
 
84
}
 
85
 
 
86
std::vector<double> NeighborJoining::computeBranchLengthsForPair(const std::vector<unsigned int>& pair)
 
87
{
 
88
  double ratio = (sumDist_[pair[0]] - sumDist_[pair[1]]) / (currentNodes_.size() - 2);
 
89
  vector<double> d(2);
 
90
  if (positiveLengths_)
 
91
  {
 
92
    d[0] = std::max(.5 * (matrix_(pair[0], pair[1]) + ratio), 0.); 
 
93
    d[1] = std::max(.5 * (matrix_(pair[0], pair[1]) - ratio), 0.); 
 
94
  }
 
95
  else
 
96
  {
 
97
    d[0] = .5 * (matrix_(pair[0], pair[1]) + ratio); 
 
98
    d[1] = .5 * (matrix_(pair[0], pair[1]) - ratio); 
 
99
  }
 
100
  return d;
 
101
}
 
102
 
 
103
double NeighborJoining::computeDistancesFromPair(const std::vector<unsigned int>& pair, const std::vector<double>& branchLengths, unsigned int pos)
 
104
{
 
105
  return 
 
106
    positiveLengths_ ?
 
107
      std::max(.5 * (matrix_(pair[0], pos) - branchLengths[0] + matrix_(pair[1], pos) - branchLengths[1]), 0.)
 
108
    :          .5 * (matrix_(pair[0], pos) - branchLengths[0] + matrix_(pair[1], pos) - branchLengths[1]); 
 
109
}
 
110
 
 
111
void NeighborJoining::finalStep(int idRoot)
 
112
{
 
113
  Node * root = new Node(idRoot);
 
114
  map<unsigned int, Node* >::iterator it = currentNodes_.begin();
 
115
  unsigned int i1 = it->first;
 
116
  Node * n1       = it->second;
 
117
  it++;
 
118
  unsigned int i2 = it->first;
 
119
  Node * n2       = it->second;
 
120
  if(currentNodes_.size() == 2)
 
121
  {
 
122
    //Rooted
 
123
    double d = matrix_(i1, i2) / 2;
 
124
    root->addSon(n1);
 
125
    root->addSon(n2);
 
126
    n1->setDistanceToFather(d);
 
127
    n2->setDistanceToFather(d);
 
128
  }
 
129
  else
 
130
  {
 
131
    //Unrooted
 
132
    it++;
 
133
    unsigned int i3 = it->first;
 
134
    Node * n3       = it->second;
 
135
    double d1 = positiveLengths_ ?
 
136
        std::max(matrix_(i1, i2) + matrix_(i1, i3) - matrix_(i2, i3), 0.)
 
137
      :          matrix_(i1, i2) + matrix_(i1, i3) - matrix_(i2, i3);
 
138
    double d2 = positiveLengths_ ?
 
139
        std::max(matrix_(i2, i1) + matrix_(i2, i3) - matrix_(i1, i3), 0.)
 
140
      :          matrix_(i2, i1) + matrix_(i2, i3) - matrix_(i1, i3);
 
141
    double d3 = positiveLengths_ ?
 
142
        std::max(matrix_(i3, i1) + matrix_(i3, i2) - matrix_(i1, i2), 0.)
 
143
      :          matrix_(i3, i1) + matrix_(i3, i2) - matrix_(i1, i2);
 
144
    root->addSon(n1);
 
145
    root->addSon(n2);
 
146
    root->addSon(n3);
 
147
    n1->setDistanceToFather(d1/2.);
 
148
    n2->setDistanceToFather(d2/2.);
 
149
    n3->setDistanceToFather(d3/2.);
 
150
  }
 
151
  tree_ = new TreeTemplate<Node>(root);
 
152
}
 
153