utils/containers/matrix/operators.hpp

Go to the documentation of this file.

#ifndef GENESIS_UTILS_CONTAINERS_MATRIX_OPERATORS_H_
#define GENESIS_UTILS_CONTAINERS_MATRIX_OPERATORS_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2018 Lucas Czech and HITS gGmbH
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
    Contact:
    Lucas Czech <lucas.czech@h-its.org>
    Exelixis Lab, Heidelberg Institute for Theoretical Studies
    Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
*/
 
#include <cassert>
#include <ostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
#include "genesis/utils/containers/matrix.hpp"
#include "genesis/utils/core/algorithm.hpp"
 
namespace genesis {
namespace utils {
 
// =================================================================================================
//     Helpful Functions
// =================================================================================================
 
std::pair<size_t, size_t> triangular_indices( size_t k, size_t n );
 
size_t triangular_index( size_t i, size_t j, size_t n );
 
size_t triangular_size( size_t n );
 
// =================================================================================================
//     General Matrix Operators
// =================================================================================================
 
template <typename T>
Matrix<T> transpose( Matrix<T> const& data )
{
    auto res = Matrix<T>( data.cols(), data.rows() );
    for( size_t r = 0; r < data.rows(); ++r ) {
        for( size_t c = 0; c < data.cols(); ++c ) {
            res( c, r ) = data( r, c );
        }
    }
    return res;
}
 
template <typename T>
bool is_square( Matrix<T> const& data )
{
    return data.rows() == data.cols();
}
 
template <typename T>
bool is_symmetric( Matrix<T> const& data )
{
    if( data.rows() != data.cols() ) {
        return false;
    }
 
    // We only need to check the upper triangle, and compare it to the lower triangle.
    // Also, we use early stopping. Of course we do. Who wouldn't.
    for( size_t i = 0; i < data.rows(); ++i ) {
        for( size_t j = i + 1; j < data.cols(); ++j ) {
            if( data( i, j ) != data( j, i ) ) {
                return false;
            }
        }
    }
 
    return true;
}
 
template <typename T>
std::ostream& operator<< (std::ostream& os, const Matrix<T>& matrix)
{
    for (size_t i = 0; i < matrix.rows(); ++i) {
        for (size_t j = 0; j < matrix.cols(); ++j) {
            os << matrix(i, j);
            if (j < matrix.cols() - 1) {
                os << " ";
            }
        }
        os << "\n";
    }
    return os;
}
 
template<>
std::ostream& operator<< (std::ostream& os, const Matrix<signed char>& matrix);
 
template<>
std::ostream& operator<< (std::ostream& os, const Matrix<unsigned char>& matrix);
 
template <typename T>
void print( std::ostream& out, Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
{
    // If the user does not want limits, or uses wrong ones, just use everything!
    if( rows == 0 || rows >= matrix.rows() ) {
        rows = matrix.rows();
    }
    if( cols == 0 || cols >= matrix.cols() ) {
        cols = matrix.cols();
    }
 
    // Print as many rows and cols as wanted.
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < cols; ++j) {
            out << matrix(i, j);
            if (j < matrix.cols() - 1) {
                out << " ";
            }
        }
        if( cols < matrix.cols() ) {
            out << " ...";
        }
        out << "\n";
    }
    if( rows < matrix.rows() ) {
        out << "...\n";
    }
}
 
template <typename T>
std::string print( Matrix<T> const& matrix, size_t rows = 10, size_t cols = 10 )
{
    std::ostringstream out;
    print( out, matrix, rows, cols );
    return out.str();
}
 
// =================================================================================================
//     Swapping
// =================================================================================================
 
template <typename T>
void swap_rows( Matrix<T>& data, size_t row_a, size_t row_b )
{
    if( row_a >= data.rows() || row_b >= data.rows() ) {
        throw std::invalid_argument( "Invalid row index for swap_rows()." );
    }
 
    using std::swap;
    for( size_t c = 0; c < data.cols(); ++c ) {
        swap( data( row_a, c ), data( row_b, c ) );
    }
}
 
template <typename T>
void swap_cols( Matrix<T>& data, size_t col_a, size_t col_b )
{
    if( col_a >= data.rows() || col_b >= data.rows() ) {
        throw std::invalid_argument( "Invalid column index for swap_cols()." );
    }
 
    using std::swap;
    for( size_t r = 0; r < data.rows(); ++r ) {
        swap( data( r, col_a ), data( r, col_b ) );
    }
}
 
} // namespace utils
} // namespace genesis
 
#endif // include guard