A toolkit for working with phylogenetic data.
v0.24.0
distance.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2018 Lucas Czech and HITS gGmbH
4 
5  This program is free software: you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 
18  Contact:
19  Lucas Czech <lucas.czech@h-its.org>
20  Exelixis Lab, Heidelberg Institute for Theoretical Studies
21  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
22 */
23 
32 
34 
35 #include <cassert>
36 
37 #ifdef GENESIS_OPENMP
38 # include <omp.h>
39 #endif
40 
41 namespace genesis {
42 namespace utils {
43 
44 // =============================================================================
45 // Distance Matrices
46 // =============================================================================
47 
49 {
50  // Init result matrix.
51  auto result = utils::Matrix<double>( data.rows(), data.rows(), 0.0 );
52 
53  // Parallel specialized code.
54  #ifdef GENESIS_OPENMP
55 
56  // We only need to calculate the upper triangle. Get the number of indices needed
57  // to describe this triangle.
58  size_t const max_k = utils::triangular_size( data.rows() );
59 
60  // Calculate.
61  #pragma omp parallel for
62  for( size_t k = 0; k < max_k; ++k ) {
63 
64  // For the given linear index, get the actual position in the Matrix.
65  auto const ij = utils::triangular_indices( k, data.rows() );
66  auto const i = ij.first;
67  auto const j = ij.second;
68 
69  // Calculate EMD and fill symmetric Matrix.
70  auto const dist = p_norm_distance(
71  data.row(i).begin(), data.row(i).end(),
72  data.row(j).begin(), data.row(j).end(),
73  p
74  );
75 
76  assert( result( i, j ) == 0.0 );
77  assert( result( j, i ) == 0.0 );
78  result( i, j ) = dist;
79  result( j, i ) = dist;
80  }
81 
82  // If no threads are available at all, use serial version.
83  #else
84 
85  for( size_t i = 0; i < data.rows() - 1; ++i ) {
86 
87  // The result is symmetric - we only calculate the upper triangle.
88  for( size_t j = i + 1; j < data.cols(); ++j ) {
89 
90  auto const dist = p_norm_distance(
91  data.row(i).begin(), data.row(i).end(),
92  data.row(j).begin(), data.row(j).end(),
93  p
94  );
95  result( i, j ) = dist;
96  result( j, i ) = dist;
97  }
98  }
99 
100  #endif
101 
102  return result;
103 }
104 
106 {
107  return p_norm_distance_matrix( data, 1.0 );
108 }
109 
111 {
112  return p_norm_distance_matrix( data, 2.0 );
113 }
114 
116 {
117  return p_norm_distance_matrix( data, std::numeric_limits<double>::infinity() );
118 }
119 
120 } // namespace utils
121 } // namespace genesis
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
Matrix operators.
double p_norm_distance(ForwardIteratorA first_a, ForwardIteratorA last_a, ForwardIteratorB first_b, ForwardIteratorB last_b, double p=2.0)
Calculate the p-norm distance between two (mathematical) vectors.
Definition: distance.hpp:287
size_t triangular_size(size_t n)
Calculate the number of linear indices needed for a triangular Matrix of size n.
MatrixRow< self_type, value_type > row(size_t row)
Matrix< double > maximum_distance_matrix(Matrix< double > const &data)
Calculate the pairwise maximum distance matrix between the rows of a given matrix.
Definition: distance.cpp:115
Matrix< double > euclidean_distance_matrix(Matrix< double > const &data)
Calculate the pairwise euclidean distance matrix between the rows of a given matrix.
Definition: distance.cpp:110
Matrix< double > p_norm_distance_matrix(Matrix< double > const &data, double p)
Calculate the pairwise distance matrix between the rows of a given matrix.
Definition: distance.cpp:48
Matrix< double > manhattan_distance_matrix(Matrix< double > const &data)
Calculate the pairwise manhatten distance matrix between the rows of a given matrix.
Definition: distance.cpp:105
std::pair< size_t, size_t > triangular_indices(size_t k, size_t n)
Given a linear index in a upper triangular Matrix, find the corresponding Matrix indices.