1 #ifndef GENESIS_UTILS_MATH_MATRIX_H_ 2 #define GENESIS_UTILS_MATH_MATRIX_H_ 54 template<
typename T >
58 if( data.
rows() == 0 || data.
cols() == 0 ) {
60 std::numeric_limits<double>::quiet_NaN(),
61 std::numeric_limits<double>::quiet_NaN()
67 std::numeric_limits<double>::infinity(),
68 - std::numeric_limits<double>::infinity()
73 for(
auto const& e : data ) {
74 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
75 ret.
min = std::min( ret.min, e );
76 ret.max = std::max( ret.max, e );
84 std::numeric_limits<double>::quiet_NaN(),
85 std::numeric_limits<double>::quiet_NaN()
98 template<
typename T >
101 auto ret = std::vector<MinMaxPair<T>>( data.
cols(), {
102 std::numeric_limits<double>::quiet_NaN(),
103 std::numeric_limits<double>::quiet_NaN()
107 if( data.
rows() == 0 ) {
112 for(
size_t c = 0; c < data.
cols(); ++c ) {
113 ret[ c ].min = std::numeric_limits<double>::infinity();
114 ret[ c ].max = - std::numeric_limits<double>::infinity();
117 for(
size_t r = 0; r < data.
rows(); ++r ) {
120 auto const& e = data( r, c );
121 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
122 ret[ c ].min = std::min( ret[ c ].min, data( r, c ) );
123 ret[ c ].max = std::max( ret[ c ].max, data( r, c ) );
130 ret[ c ].min = std::numeric_limits<double>::quiet_NaN();
131 ret[ c ].max = std::numeric_limits<double>::quiet_NaN();
144 template<
typename T >
147 auto ret = std::vector<MinMaxPair<T>>( data.
rows(), {
148 std::numeric_limits<double>::quiet_NaN(),
149 std::numeric_limits<double>::quiet_NaN()
153 if( data.
cols() == 0 ) {
157 for(
size_t r = 0; r < data.
rows(); ++r ) {
159 ret[ r ].min = std::numeric_limits<double>::infinity();
160 ret[ r ].max = - std::numeric_limits<double>::infinity();
164 for(
size_t c = 0; c < data.
cols(); ++c ) {
165 auto const& e = data( r, c );
166 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
167 ret[ r ].min = std::min( ret[ r ].min, data( r, c ) );
168 ret[ r ].max = std::max( ret[ r ].max, data( r, c ) );
175 ret[ r ].min = std::numeric_limits<double>::quiet_NaN();
176 ret[ r ].max = std::numeric_limits<double>::quiet_NaN();
189 template<
typename T >
194 for(
auto const& e : data ) {
195 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
208 template<
typename T >
212 auto row_sums = std::vector<T>( data.
rows(), T{} );
213 for(
size_t i = 0; i < data.
rows(); ++i ) {
214 for(
size_t j = 0; j < data.
cols(); ++j ) {
215 auto const& e = data( i, j );
216 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
230 template<
typename T >
234 auto col_sums = std::vector<T>( data.
cols(), T{} );
235 for(
size_t i = 0; i < data.
rows(); ++i ) {
236 for(
size_t j = 0; j < data.
cols(); ++j ) {
237 auto const& e = data( i, j );
238 if( std::isfinite( e ) || ! ignore_non_finite_values ) {
297 bool scale_means =
true,
298 bool scale_std =
true 320 bool scale_means =
true,
321 bool scale_std =
true 340 double epsilon = 1e-5
383 template<
typename T >
387 throw std::runtime_error(
"Symmetric sort only works on square matrices." );
394 for(
size_t i = 0; i < data.
rows(); ++i ) {
395 for(
size_t j = 0; j < data.
cols(); ++j ) {
396 result( i, j ) = data( si[i], si[j] );
414 template<
typename T >
418 throw std::runtime_error(
"Symmetric sort only works on square matrices." );
425 for(
size_t i = 0; i < data.
rows(); ++i ) {
426 for(
size_t j = 0; j < data.
cols(); ++j ) {
427 result( i, j ) = data( si[i], si[j] );
442 template<
typename T >
446 throw std::runtime_error(
"Symmetric sort only works on square matrices." );
451 auto find_max = [](
Matrix<T> const& mat,
size_t start ){
452 size_t max_r = start;
453 size_t max_c = start;
454 for(
size_t r = start; r < mat.
rows(); ++r ) {
455 for(
size_t c = start; c < mat.
cols(); ++c ) {
456 if( mat(r, c) > mat( max_r, max_c ) ) {
462 return std::make_pair( max_r, max_c );
468 for(
size_t i = 0; i < mat.
rows(); ++i ) {
469 auto const max_pair = find_max( mat, i );
470 matrix_swap_rows( mat, i, max_pair.first );
471 matrix_swap_cols( mat, i, max_pair.second );
485 template<
typename T =
double,
typename A =
double,
typename B =
double >
489 throw std::runtime_error(
"Cannot add matrices with different dimensions." );
493 for(
size_t r = 0; r < a.
rows(); ++r ) {
494 for(
size_t c = 0; c < a.
cols(); ++c ) {
495 result( r, c ) = a( r, c ) + b( r, c );
504 template<
typename T =
double,
typename A =
double,
typename B =
double >
508 for(
size_t r = 0; r < matrix.
rows(); ++r ) {
509 for(
size_t c = 0; c < matrix.
cols(); ++c ) {
510 result( r, c ) = matrix( r, c ) + scalar;
521 template<
typename T =
double,
typename A =
double,
typename B =
double >
525 throw std::runtime_error(
"Cannot add matrices with different dimensions." );
529 for(
size_t r = 0; r < a.
rows(); ++r ) {
530 for(
size_t c = 0; c < a.
cols(); ++c ) {
531 result( r, c ) = a( r, c ) - b( r, c );
547 template<
typename T =
double,
typename A =
double,
typename B =
double >
551 throw std::runtime_error(
"Cannot multiply matrices if a.cols() != b.rows()." );
558 for(
size_t r = 0; r < a.
rows(); ++r ) {
559 for(
size_t c = 0; c < b.
cols(); ++c ) {
560 result( r, c ) = 0.0;
561 for(
size_t j = 0; j < a.
cols(); ++j ) {
562 result( r, c ) += a( r, j ) * b( j, c );
577 template<
typename T =
double,
typename A =
double,
typename B =
double >
580 if( a.size() != b.
rows() ) {
581 throw std::runtime_error(
"Cannot multiply vector with matrix if a.size() != b.rows()." );
584 auto result = std::vector<T>( b.
cols(), 0.0 );
585 for(
size_t c = 0; c < b.
cols(); ++c ) {
586 for(
size_t j = 0; j < a.size(); ++j ) {
587 result[ c ] += a[ j ] * b( j, c );
601 template<
typename T =
double,
typename A =
double,
typename B =
double >
604 if( a.
cols() != b.size() ) {
605 throw std::runtime_error(
"Cannot multiply matrix with vector if a.cols() != b.size()." );
608 auto result = std::vector<T>( a.
rows(), 0.0 );
609 for(
size_t r = 0; r < a.
rows(); ++r ) {
610 for(
size_t j = 0; j < a.
cols(); ++j ) {
611 result[ r ] += a( r, j ) * b[ j ];
621 template<
typename T =
double,
typename A =
double,
typename B =
double >
625 for(
size_t r = 0; r < matrix.
rows(); ++r ) {
626 for(
size_t c = 0; c < matrix.
cols(); ++c ) {
627 result( r, c ) = matrix( r, c ) * scalar;
636 #endif // include guard Matrix< T > matrix_multiplication(Matrix< A > const &a, Matrix< B > const &b)
Calculate the product of two Matrices.
std::vector< MinMaxPair< T > > matrix_row_minmax(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the row-wise min and max values of a Matrix.
Matrix< double > sums_of_squares_and_cross_products_matrix(Matrix< double > const &data)
Calculate the Sums of Squares and Cross Products Matrix (SSCP Matrix).
Store a pair of min and max values.
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.
std::vector< MeanStddevPair > standardize_rows(Matrix< double > &data, bool scale_means, bool scale_std)
Standardize the rows of a Matrix by subtracting the mean and scaling to unit variance.
Matrix< T > matrix_subtraction(Matrix< A > const &a, Matrix< B > const &b)
Calculate the element-wise difference of two Matrices.
std::vector< T > matrix_row_sums(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the sum of each row and return the result as a vector.
double sum(const Histogram &h)
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
std::vector< size_t > sort_indices(RandomAccessIterator first, RandomAccessIterator last, Comparator comparator)
Get the indices to the sorted order of the given range.
std::vector< T > matrix_col_sums(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the sum of each column and return the result as a vector.
T matrix_sum(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the sum of all elements in a Matrix.
MinMaxPair< T > matrix_minmax(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the min and max values of a Matrix.
Matrix< double > covariance_matrix(Matrix< double > const &data)
Calculate the covariance Matrix of a given data Matrix.
std::vector< MinMaxPair< double > > normalize_cols(Matrix< double > &data)
Normalize the columns of a Matrix so that all values are in the range [ 0.0, 1.0 ].
std::vector< MinMaxPair< double > > normalize_rows(Matrix< double > &data)
Normalize the rows of a Matrix so that all values are in the range [ 0.0, 1.0 ].
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.
std::vector< MinMaxPair< T > > matrix_col_minmax(Matrix< T > const &data, bool ignore_non_finite_values=true)
Calculate the column-wise min and max values of a Matrix.
std::vector< size_t > filter_constant_columns(utils::Matrix< double > &data, double epsilon)
Filter out columns that have nearly constant values, measured using an epsilon.
Matrix< T > matrix_addition(Matrix< A > const &a, Matrix< B > const &b)
Calculate the element-wise sum of two Matrices.
Matrix< T > matrix_sort_diagonal_symmetric(Matrix< T > const &data)
Sort a Matrix so that the highest entries are on the diagonal.
std::vector< MeanStddevPair > standardize_cols(Matrix< double > &data, bool scale_means, bool scale_std)
Standardize the columns of a Matrix by subtracting the mean and scaling to unit variance.
Matrix< double > correlation_matrix(Matrix< double > const &data)
Calculate the correlation Matrix of a given data Matrix.