A library for working with phylogenetic and population genetic data.
v0.32.0
utils/containers/matrix/operators.hpp
Go to the documentation of this file.
1 #ifndef GENESIS_UTILS_CONTAINERS_MATRIX_OPERATORS_H_
2 #define GENESIS_UTILS_CONTAINERS_MATRIX_OPERATORS_H_
3 
4 /*
5  Genesis - A toolkit for working with phylogenetic data.
6  Copyright (C) 2014-2023 Lucas Czech
7 
8  This program is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  This program is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 
21  Contact:
22  Lucas Czech <lczech@carnegiescience.edu>
23  Department of Plant Biology, Carnegie Institution For Science
24  260 Panama Street, Stanford, CA 94305, USA
25 */
26 
34 #include <cassert>
35 #include <ostream>
36 #include <sstream>
37 #include <stdexcept>
38 #include <string>
39 #include <type_traits>
40 #include <utility>
41 #include <vector>
42 
45 
46 namespace genesis {
47 namespace utils {
48 
49 // =================================================================================================
50 // Helpful Functions
51 // =================================================================================================
52 
83 std::pair<size_t, size_t> triangular_indices( size_t k, size_t n );
84 
96 size_t triangular_index( size_t i, size_t j, size_t n );
97 
117 size_t triangular_size( size_t n );
118 
119 // =================================================================================================
120 // General Matrix Operators
121 // =================================================================================================
122 
126 template <typename T>
128 {
129  auto res = Matrix<T>( mat.cols(), mat.rows() );
130  for( size_t r = 0; r < mat.rows(); ++r ) {
131  for( size_t c = 0; c < mat.cols(); ++c ) {
132  res( c, r ) = mat( r, c );
133  }
134  }
135  return res;
136 }
137 
145 template <typename T>
147 {
148  // We use a follow-the-cycles implementation inspired by https://stackoverflow.com/a/9320349
149  // That description seems to use a flipped notation by expecting n x m matrix.
150 
151  // Need additional storage to see which entries have been cycled already.
152  std::vector<bool> visited( mat.data().size() );
153  size_t const div = mat.data().size() - 1;
154 
155  // Cycle through the matrix until we have moved every entry to its destination.
156  size_t cycle = 0;
157  while( ++cycle < mat.data().size() ) {
158  if( visited[cycle] ){
159  continue;
160  }
161  size_t cur = cycle;
162  do {
163  cur = cur == div ? div : (mat.rows() * cur) % div;
164  std::swap( mat.data_[cur], mat.data_[cycle]);
165  visited[cur] = true;
166  } while( cur != cycle );
167  }
168 
169  // Finally we need to update the dimensions of the matrix.
170  std::swap( mat.rows_, mat.cols_ );
171 }
172 
177 template <typename T>
178 bool is_square( Matrix<T> const& data )
179 {
180  return data.rows() == data.cols();
181 }
182 
187 template <typename T>
188 bool is_symmetric( Matrix<T> const& data )
189 {
190  if( data.rows() != data.cols() ) {
191  return false;
192  }
193 
194  // We only need to check the upper triangle, and compare it to the lower triangle.
195  // Also, we use early stopping. Of course we do. Who wouldn't.
196  for( size_t i = 0; i < data.rows(); ++i ) {
197  for( size_t j = i + 1; j < data.cols(); ++j ) {
198  if( data( i, j ) != data( j, i ) ) {
199  return false;
200  }
201  }
202  }
203 
204  return true;
205 }
206 
210 template<typename T>
211 std::ostream& operator<< ( std::ostream& os, const Matrix<T>& matrix )
212 {
213  for (size_t i = 0; i < matrix.rows(); ++i) {
214  for (size_t j = 0; j < matrix.cols(); ++j) {
215  os << matrix(i, j);
216  if (j < matrix.cols() - 1) {
217  os << " ";
218  }
219  }
220  os << "\n";
221  }
222  return os;
223 }
224 
228 template<>
229 std::ostream& operator<< (std::ostream& os, const Matrix<signed char>& matrix);
230 
234 template<>
235 std::ostream& operator<< (std::ostream& os, const Matrix<unsigned char>& matrix);
236 
241 template <typename T>
242 void print( std::ostream& out, Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
243 {
244  // If the user does not want limits, or uses wrong ones, just use everything!
245  if( rows == 0 || rows >= matrix.rows() ) {
246  rows = matrix.rows();
247  }
248  if( cols == 0 || cols >= matrix.cols() ) {
249  cols = matrix.cols();
250  }
251 
252  // Print as many rows and cols as wanted.
253  for (size_t i = 0; i < rows; ++i) {
254  for (size_t j = 0; j < cols; ++j) {
255 
256  // We need some special printing for int char types...
257  if( std::is_same<T, signed char>::value || std::is_same<T, unsigned char>::value ) {
258  out << static_cast<int>( matrix( i, j ));
259  } else {
260  out << matrix(i, j);
261  }
262  if (j < matrix.cols() - 1) {
263  out << " ";
264  }
265  }
266  if( cols < matrix.cols() ) {
267  out << " ...";
268  }
269  out << "\n";
270  }
271  if( rows < matrix.rows() ) {
272  out << "...\n";
273  }
274 }
275 
284 template <typename T>
285 std::string print( Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
286 {
287  std::ostringstream out;
288  print( out, matrix, rows, cols );
289  return out.str();
290 }
291 
292 // =================================================================================================
293 // Swapping
294 // =================================================================================================
295 
299 template <typename T>
300 void swap_rows( Matrix<T>& data, size_t row_a, size_t row_b )
301 {
302  if( row_a >= data.rows() || row_b >= data.rows() ) {
303  throw std::invalid_argument( "Invalid row index for swap_rows()." );
304  }
305 
306  using std::swap;
307  for( size_t c = 0; c < data.cols(); ++c ) {
308  swap( data( row_a, c ), data( row_b, c ) );
309  }
310 }
311 
315 template <typename T>
316 void swap_cols( Matrix<T>& data, size_t col_a, size_t col_b )
317 {
318  if( col_a >= data.rows() || col_b >= data.rows() ) {
319  throw std::invalid_argument( "Invalid column index for swap_cols()." );
320  }
321 
322  using std::swap;
323  for( size_t r = 0; r < data.rows(); ++r ) {
324  swap( data( r, col_a ), data( r, col_b ) );
325  }
326 }
327 
328 } // namespace utils
329 } // namespace genesis
330 
331 #endif // include guard
genesis::placement::swap
void swap(Sample &lhs, Sample &rhs)
Definition: sample.cpp:104
algorithm.hpp
Provides some valuable algorithms that are not part of the C++ 11 STL.
genesis::utils::triangular_index
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.
Definition: utils/containers/matrix/operators.cpp:56
genesis::utils::Matrix::cols
size_t cols() const
Definition: containers/matrix.hpp:181
genesis::utils::is_symmetric
bool is_symmetric(Matrix< T > const &data)
Return whether a Matrix is symmetric, i.e., whether it is square and m[ i, j ] == m[ j,...
Definition: utils/containers/matrix/operators.hpp:188
genesis::utils::print
void print(std::ostream &out, Matrix< T > const &matrix, size_t rows=10, size_t cols=10)
Print a Matrix to an out stream. See print( Matrix<T> const&, size_t, size_t ) for details.
Definition: utils/containers/matrix/operators.hpp:242
genesis::utils::transpose_inplace
void transpose_inplace(Matrix< T > &mat)
Transpose a Matrix inplace, without allocating a new Matrix.
Definition: utils/containers/matrix/operators.hpp:146
genesis::utils::operator<<
std::ostream & operator<<(std::ostream &os, Color const &color)
Write a textual representation of the Color the a stream, in the format "(r, g, b,...
Definition: utils/color/functions.cpp:129
genesis::utils::triangular_indices
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.
Definition: utils/containers/matrix/operators.cpp:44
genesis::utils::swap
void swap(Color &lhs, Color &rhs)
Definition: color.hpp:207
genesis::utils::Matrix
Definition: placement/function/emd.hpp:53
genesis::utils::is_square
bool is_square(Matrix< T > const &data)
Return whether a Matrix is a square matrix, that is, whether its number of rows and number of columns...
Definition: utils/containers/matrix/operators.hpp:178
genesis::utils::Matrix::data
container_type const & data() const
Definition: containers/matrix.hpp:196
matrix.hpp
genesis
Container namespace for all symbols of genesis in order to keep them separate when used as a library.
Definition: placement/formats/edge_color.cpp:42
genesis::utils::swap_rows
void swap_rows(Matrix< T > &data, size_t row_a, size_t row_b)
Swap (interchange) two rows of a Matrix, given their indices.
Definition: utils/containers/matrix/operators.hpp:300
genesis::utils::Matrix::rows
size_t rows() const
Definition: containers/matrix.hpp:176
genesis::utils::triangular_size
size_t triangular_size(size_t n)
Calculate the number of linear indices needed for a triangular Matrix of size n.
Definition: utils/containers/matrix/operators.cpp:62
genesis::utils::swap_cols
void swap_cols(Matrix< T > &data, size_t col_a, size_t col_b)
Swap (interchange) two columns of a Matrix, given their indices.
Definition: utils/containers/matrix/operators.hpp:316
genesis::utils::transpose
Matrix< T > transpose(Matrix< T > const &mat)
Transpose a Matrix.
Definition: utils/containers/matrix/operators.hpp:127