A toolkit for working with phylogenetic data.
v0.24.0
utils/containers/matrix/operators.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2019 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 
33 #include <cmath>
34 #include <iomanip>
35 #include <stdexcept>
36 
37 namespace genesis {
38 namespace utils {
39 
40 // ================================================================================================
41 // Helpful Functions
42 // ================================================================================================
43 
44 std::pair<size_t, size_t> triangular_indices( size_t k, size_t n )
45 {
46  // Using equations from http://stackoverflow.com/a/27088560/4184258
47  // See also https://en.wikipedia.org/wiki/Triangular_number
48 
49  size_t const i = n - 2 - static_cast<size_t>(
50  std::floor( std::sqrt( 4 * n * ( n - 1 ) - 7 - ( 8 * k )) / 2.0 - 0.5 )
51  );
52  size_t const j = k + i + 1 - n * ( n - 1 ) / 2 + ( n - i ) * (( n - i ) - 1 ) / 2;
53  return { i, j };
54 }
55 
56 size_t triangular_index( size_t i, size_t j, size_t n )
57 {
58  size_t const k = ( n * ( n - 1 ) / 2 ) - ( n - i ) * (( n - i ) - 1 ) / 2 + j - i - 1;
59  return k;
60 }
61 
62 size_t triangular_size( size_t n )
63 {
64  return ( n * n - n ) / 2;
65 }
66 
67 // ================================================================================================
68 // Output Operators
69 // ================================================================================================
70 
74 template <typename T>
75 void print_byte_matrix_( std::ostream& os, const Matrix<T>& matrix, size_t width )
76 {
77  for (size_t i = 0; i < matrix.rows(); ++i) {
78  for (size_t j = 0; j < matrix.cols(); ++j) {
79  os << std::setw(width) << std::setfill(' ') << static_cast<int>( matrix(i, j) );
80  if (j < matrix.cols() - 1) {
81  os << " ";
82  }
83  }
84  os << "\n";
85  }
86 }
87 
88 template<>
89 std::ostream& operator<< (std::ostream& os, const Matrix<signed char>& matrix)
90 {
91  print_byte_matrix_( os, matrix, 4 );
92  return os;
93 }
94 
95 template<>
96 std::ostream& operator<< (std::ostream& os, const Matrix<unsigned char>& matrix)
97 {
98  print_byte_matrix_( os, matrix, 3 );
99  return os;
100 }
101 
102 } // namespace utils
103 } // namespace genesis
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
Matrix operators.
void print_byte_matrix_(std::ostream &os, const Matrix< T > &matrix, size_t width)
Local helper function to avoid code duplication.
size_t triangular_size(size_t n)
Calculate the number of linear indices needed for a triangular Matrix of size n.
size_t triangular_index(size_t i, size_t j, size_t n)
Given indices i and j in a quadratic Matrix, find the corresponding linear index. ...
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.