A toolkit for working with phylogenetic data.
v0.18.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
utils/math/matrix/operators.hpp
Go to the documentation of this file.
1 #ifndef GENESIS_UTILS_MATH_MATRIX_OPERATORS_H_
2 #define GENESIS_UTILS_MATH_MATRIX_OPERATORS_H_
3 
4 /*
5  Genesis - A toolkit for working with phylogenetic data.
6  Copyright (C) 2014-2017 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 <lucas.czech@h-its.org>
23  Exelixis Lab, Heidelberg Institute for Theoretical Studies
24  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
25 */
26 
34 #include <ostream>
35 #include <sstream>
36 #include <stdexcept>
37 #include <string>
38 #include <vector>
39 
43 
44 namespace genesis {
45 namespace utils {
46 
47 // =================================================================================================
48 // Helpful Functions
49 // =================================================================================================
50 
81 std::pair<size_t, size_t> triangular_indices( size_t k, size_t n );
82 
94 size_t triangular_index( size_t i, size_t j, size_t n );
95 
115 size_t triangular_size( size_t n );
116 
117 // =================================================================================================
118 // General Matrix Operators
119 // =================================================================================================
120 
124 template <typename T>
126 {
127  auto res = Matrix<T>( data.cols(), data.rows() );
128  for( size_t r = 0; r < data.rows(); ++r ) {
129  for( size_t c = 0; c < data.cols(); ++c ) {
130  res( c, r ) = data( r, c );
131  }
132  }
133  return res;
134 }
135 
140 template <typename T>
141 bool matrix_is_symmetric( Matrix<T> const& data )
142 {
143  if( data.rows() != data.cols() ) {
144  return false;
145  }
146 
147  // We only need to check the upper triangle, and compare it to the lower triangle.
148  // Also, we use early stopping. Of course we do. Who wouldn't.
149  for( size_t i = 0; i < data.rows(); ++i ) {
150  for( size_t j = i + 1; j < data.cols(); ++j ) {
151  if( data( i, j ) != data( j, i ) ) {
152  return false;
153  }
154  }
155  }
156 
157  return true;
158 }
159 
160 template <typename T>
161 std::ostream& operator<< (std::ostream& os, const Matrix<T>& matrix)
162 {
163  for (size_t i = 0; i < matrix.rows(); ++i) {
164  for (size_t j = 0; j < matrix.cols(); ++j) {
165  os << matrix(i, j);
166  if (j < matrix.cols() - 1) {
167  os << " ";
168  }
169  }
170  os << "\n";
171  }
172  return os;
173 }
174 
179 template <typename T>
180 void print( std::ostream& out, Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
181 {
182  // If the user does not want limits, or uses wrong ones, just use everything!
183  if( rows == 0 || rows >= matrix.rows() ) {
184  rows = matrix.rows();
185  }
186  if( cols == 0 || cols >= matrix.cols() ) {
187  cols = matrix.cols();
188  }
189 
190  // Print as many rows and cols as wanted.
191  for (size_t i = 0; i < rows; ++i) {
192  for (size_t j = 0; j < cols; ++j) {
193  out << matrix(i, j);
194  if (j < matrix.cols() - 1) {
195  out << " ";
196  }
197  }
198  if( cols < matrix.cols() ) {
199  out << " ...";
200  }
201  out << "\n";
202  }
203  if( rows < matrix.rows() ) {
204  out << "...\n";
205  }
206 }
207 
216 template <typename T>
217 std::string print( Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
218 {
219  std::ostringstream out;
220  print( out, matrix, rows, cols );
221  return out.str();
222 }
223 
224 // =================================================================================================
225 // Swapping and Sorting
226 // =================================================================================================
227 
231 template <typename T>
232 void matrix_swap_rows( Matrix<T>& data, size_t row_a, size_t row_b )
233 {
234  if( row_a >= data.rows() || row_b >= data.rows() ) {
235  throw std::invalid_argument( "Invalid row index for swap_rows()." );
236  }
237 
238  using std::swap;
239  for( size_t c = 0; c < data.cols(); ++c ) {
240  swap( data( row_a, c ), data( row_b, c ) );
241  }
242 }
243 
247 template <typename T>
248 void matrix_swap_cols( Matrix<T>& data, size_t col_a, size_t col_b )
249 {
250  if( col_a >= data.rows() || col_b >= data.rows() ) {
251  throw std::invalid_argument( "Invalid column index for swap_cols()." );
252  }
253 
254  using std::swap;
255  for( size_t r = 0; r < data.rows(); ++r ) {
256  swap( data( r, col_a ), data( r, col_b ) );
257  }
258 }
259 
271 template< typename T >
273 {
274  if( data.rows() != data.cols() ) {
275  throw std::runtime_error( "Symmetric sort only works on square matrices." );
276  }
277 
278  auto result = Matrix<T>( data.rows(), data.cols() );
279  auto row_sums = matrix_row_sums( data );
280  auto si = sort_indices( row_sums.begin(), row_sums.end() );
281 
282  for( size_t i = 0; i < data.rows(); ++i ) {
283  for( size_t j = 0; j < data.cols(); ++j ) {
284  result( i, j ) = data( si[i], si[j] );
285  }
286  }
287 
288  return result;
289 }
290 
302 template< typename T >
304 {
305  if( data.rows() != data.cols() ) {
306  throw std::runtime_error( "Symmetric sort only works on square matrices." );
307  }
308 
309  auto result = Matrix<T>( data.rows(), data.cols() );
310  auto col_sums = matrix_col_sums( data );
311  auto si = sort_indices( col_sums.begin(), col_sums.end() );
312 
313  for( size_t i = 0; i < data.rows(); ++i ) {
314  for( size_t j = 0; j < data.cols(); ++j ) {
315  result( i, j ) = data( si[i], si[j] );
316  }
317  }
318 
319  return result;
320 }
321 
322 // =================================================================================================
323 // Matrix Addition
324 // =================================================================================================
325 
331 template< typename T = double, typename A = double, typename B = double >
333 {
334  if( a.rows() != b.rows() || a.cols() != b.cols() ) {
335  throw std::runtime_error( "Cannot add matrices with different dimensions." );
336  }
337 
338  auto result = Matrix<T>( a.rows(), a.cols() );
339  for( size_t r = 0; r < a.rows(); ++r ) {
340  for( size_t c = 0; c < a.cols(); ++c ) {
341  result( r, c ) = a( r, c ) + b( r, c );
342  }
343  }
344  return result;
345 }
346 
350 template< typename T = double, typename A = double, typename B = double >
351 Matrix<T> matrix_addition( Matrix<A> const& matrix, B const& scalar )
352 {
353  auto result = Matrix<T>( matrix.rows(), matrix.cols() );
354  for( size_t r = 0; r < matrix.rows(); ++r ) {
355  for( size_t c = 0; c < matrix.cols(); ++c ) {
356  result( r, c ) = matrix( r, c ) + scalar;
357  }
358  }
359  return result;
360 }
361 
367 template< typename T = double, typename A = double, typename B = double >
369 {
370  if( a.rows() != b.rows() || a.cols() != b.cols() ) {
371  throw std::runtime_error( "Cannot add matrices with different dimensions." );
372  }
373 
374  auto result = Matrix<T>( a.rows(), a.cols() );
375  for( size_t r = 0; r < a.rows(); ++r ) {
376  for( size_t c = 0; c < a.cols(); ++c ) {
377  result( r, c ) = a( r, c ) - b( r, c );
378  }
379  }
380  return result;
381 }
382 
383 // =================================================================================================
384 // Matrix Multiplication
385 // =================================================================================================
386 
393 template< typename T = double, typename A = double, typename B = double >
395 {
396  if( a.cols() != b.rows() ) {
397  throw std::runtime_error( "Cannot multiply matrices if a.cols() != b.rows()." );
398  }
399 
400  // Simple and naive. Fast enough for the few occasions were we need this.
401  // If Genesis at some point starts to need more elaborate matrix operations, it might be
402  // worth including some proper library for this.
403  auto result = Matrix<T>( a.rows(), b.cols() );
404  for( size_t r = 0; r < a.rows(); ++r ) {
405  for( size_t c = 0; c < b.cols(); ++c ) {
406  result( r, c ) = 0.0;
407  for( size_t j = 0; j < a.cols(); ++j ) {
408  result( r, c ) += a( r, j ) * b( j, c );
409  }
410  }
411  }
412 
413  return result;
414 }
415 
423 template< typename T = double, typename A = double, typename B = double >
424 std::vector<T> matrix_multiplication( std::vector<A> const& a, Matrix<B> const& b )
425 {
426  if( a.size() != b.rows() ) {
427  throw std::runtime_error( "Cannot multiply vector with matrix if a.size() != b.rows()." );
428  }
429 
430  auto result = std::vector<T>( b.cols(), 0.0 );
431  for( size_t c = 0; c < b.cols(); ++c ) {
432  for( size_t j = 0; j < a.size(); ++j ) {
433  result[ c ] += a[ j ] * b( j, c );
434  }
435  }
436 
437  return result;
438 }
439 
447 template< typename T = double, typename A = double, typename B = double >
448 std::vector<T> matrix_multiplication( Matrix<A> const& a, std::vector<B> const& b )
449 {
450  if( a.cols() != b.size() ) {
451  throw std::runtime_error( "Cannot multiply matrix with vector if a.cols() != b.size()." );
452  }
453 
454  auto result = std::vector<T>( a.rows(), 0.0 );
455  for( size_t r = 0; r < a.rows(); ++r ) {
456  for( size_t j = 0; j < a.cols(); ++j ) {
457  result[ r ] += a( r, j ) * b[ j ];
458  }
459  }
460 
461  return result;
462 }
463 
467 template< typename T = double, typename A = double, typename B = double >
468 Matrix<T> matrix_multiplication( Matrix<A> const& matrix, B const& scalar )
469 {
470  auto result = Matrix<T>( matrix.rows(), matrix.cols() );
471  for( size_t r = 0; r < matrix.rows(); ++r ) {
472  for( size_t c = 0; c < matrix.cols(); ++c ) {
473  result( r, c ) = matrix( r, c ) * scalar;
474  }
475  }
476  return result;
477 }
478 
479 } // namespace utils
480 } // namespace genesis
481 
482 #endif // include guard
Matrix< T > matrix_multiplication(Matrix< A > const &a, Matrix< B > const &b)
Calculate the product of two Matrices.
size_t cols() const
Definition: matrix.hpp:156
Provides some valuable algorithms that are not part of the C++ 11 STL.
void matrix_swap_cols(Matrix< T > &data, size_t col_a, size_t col_b)
Swap (interchange) two columns of a Matrix, given their indices.
void swap(SequenceSet &lhs, SequenceSet &rhs)
Matrix< T > matrix_sort_by_row_sum_symmetric(Matrix< T > const &data)
Sort rows and columns of a Matrix by the sum or the rows.
void matrix_swap_rows(Matrix< T > &data, size_t row_a, size_t row_b)
Swap (interchange) two rows of a Matrix, given their indices.
Matrix< T > matrix_subtraction(Matrix< A > const &a, Matrix< B > const &b)
Calculate the element-wise difference of two Matrices.
std::vector< size_t > sort_indices(RandomAccessIterator first, RandomAccessIterator last, Comparator comparator)
Get the indices to the sorted order of the given range.
Definition: algorithm.hpp:92
void swap(NexusTaxa &lhs, NexusTaxa &rhs)
Definition: taxa.hpp:207
std::vector< T > matrix_col_sums(Matrix< T > const &data)
Calculate the sum of each column and return the result as a vector.
Definition: statistics.hpp:170
Matrix statstic functions.
size_t triangular_size(size_t n)
Calculate the number of linear indices needed for a triangular Matrix of size n.
std::vector< T > matrix_row_sums(Matrix< T > const &data)
Calculate the sum of each row and return the result as a vector.
Definition: statistics.hpp:152
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.
size_t rows() const
Definition: matrix.hpp:151
Matrix< T > matrix_sort_by_col_sum_symmetric(Matrix< T > const &data)
Sort rows and columns of a Matrix by the sum or the columns.
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...
bool matrix_is_symmetric(Matrix< T > const &data)
Return whether a Matrix is symmetric, i.e., whether it is square and m[ i, j ] == m[ j...
Matrix< T > matrix_addition(Matrix< A > const &a, Matrix< B > const &b)
Calculate the element-wise sum of two Matrices.
Matrix< T > matrix_transpose(Matrix< T > const &data)
Transpose a Matrix.