Matrix.h

Go to the documentation of this file.

#ifndef MATRIX_H
#define MATRIX_H
 
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <vector>
#include <array>
#include <memory>
#include <stdexcept>
#include <algorithm>
 
#include "../special/special_generic.h"
 
// -----------------------------------------------------------------------------
// RAII wrappers around GSL resources
// -----------------------------------------------------------------------------
 
using Vector = std::vector<double>;
 
using gsl_matrix_sptr = std::unique_ptr<gsl_matrix, decltype(&gsl_matrix_free)>;
 
using gsl_vector_sptr = std::unique_ptr<gsl_vector, decltype(&gsl_vector_free)>;
 
using gsl_permutation_sptr = std::unique_ptr<gsl_permutation, decltype(&gsl_permutation_free)>;
 
using gsl_eigen_workspace_sptr = std::unique_ptr<gsl_eigen_symmv_workspace, decltype(&gsl_eigen_symmv_free)>;
 
inline gsl_matrix_sptr make_gsl_matrix(size_t rows, size_t cols) {
    return gsl_matrix_sptr(gsl_matrix_alloc(rows, cols), &gsl_matrix_free);
}
 
inline gsl_vector_sptr make_gsl_vector(size_t size) {
    return gsl_vector_sptr(gsl_vector_alloc(size), &gsl_vector_free);
}
 
inline gsl_permutation_sptr make_gsl_permutation(size_t size) {
    return gsl_permutation_sptr(gsl_permutation_alloc(size), &gsl_permutation_free);
}
 
inline gsl_eigen_workspace_sptr make_eigen_workspace(size_t n) {
    return { gsl_eigen_symmv_alloc(n), &gsl_eigen_symmv_free };
}
 
 
// -------------------------------------------------------
 
struct EigenSystem; 
 
struct SignedLogDet {
    double logdet;  
    int sign;       
};
 
class RealMatrix {
public:
    RealMatrix() : rows_(0), cols_(0) {}
 
    RealMatrix(std::vector<double> data_, std::size_t rows, std::size_t cols);
 
    RealMatrix(std::vector<std::vector<double>> data_);
 
    RealMatrix(std::size_t rows, std::size_t cols);
 
     // -------------------------------------------------------------------------
    // Element access / shape
    // -------------------------------------------------------------------------
 
    double& at(size_t i, size_t j);
 
    double& unchecked_at(size_t i, size_t j);
 
    const double& at(size_t i, size_t j) const;
 
    const double& unchecked_at(size_t i, size_t j) const;
 
    std::size_t rows() const;
 
    std::size_t cols() const;
 
    // -------------------------------------------------------------------------
    // Row / column manipulation
    // -------------------------------------------------------------------------
 
    void remove_row(std::size_t row_idx);
 
    void remove_column(std::size_t col_idx);
 
    void remove_row_and_column(std::size_t dim_idx);
 
    // -------------------------------------------------------------------------
    // GSL support
    // -------------------------------------------------------------------------
 
    static RealMatrix from_gsl_copy(const gsl_matrix* A);
 
    gsl_matrix_sptr to_gsl_matrix() const;
 
    // -------------------------------------------------------------------------
    // Linear algebra
    // -------------------------------------------------------------------------
 
    bool is_symmetric() const;
 
    RealMatrix transpose() const;
 
    EigenSystem eig() const;
 
    SignedLogDet slogdet() const;
 
    RealMatrix inv() const;
 
    // -------------------------------------------------------------------------
    // Arithmetic
    // -------------------------------------------------------------------------
 
    RealMatrix operator-() const;
 
    RealMatrix& operator+=(const RealMatrix& rhs);
 
    RealMatrix& operator-=(const RealMatrix& rhs);
 
    RealMatrix& operator*=(const RealMatrix& rhs);
 
    RealMatrix& operator*=(double scalar);
 
    RealMatrix& operator/=(double scalar);
 
    friend RealMatrix operator+(RealMatrix lhs, const RealMatrix& rhs) {
        lhs += rhs;
        return lhs;
    }
 
    friend RealMatrix operator-(RealMatrix lhs, const RealMatrix& rhs) {
        lhs -= rhs;
        return lhs;
    }
 
    friend RealMatrix operator*(RealMatrix lhs, const RealMatrix& rhs) {
        lhs *= rhs;
        return lhs;
    }
 
    friend RealMatrix operator*(RealMatrix lhs, double scalar) {
        lhs *= scalar;
        return lhs;
    }
 
    friend RealMatrix operator*(double scalar, RealMatrix rhs) {
        rhs *= scalar;
        return rhs;
    }
 
    friend RealMatrix operator/(RealMatrix lhs, double scalar) {
        lhs /= scalar;
        return lhs;
    }
 
    friend std::ostream& operator<<(std::ostream& os, RealMatrix A);
    
private:
    std::vector<double> data;   
    std::size_t rows_, cols_;   
};
 
struct EigenSystem {
    RealMatrix D;   
    RealMatrix P;   
};
 
RealMatrix eye(std::size_t n);
 
RealMatrix diag(const gsl_vector* X);
 
RealMatrix nearest_psd(RealMatrix R, double thr = 1e-12);
 
RealMatrix cholesky_L(RealMatrix R);
 
RealMatrix block_diag(const std::vector<RealMatrix>& blocks);
 
template<typename... Ms>
RealMatrix block_diag(const Ms&... ms) {
    return block_diag(std::vector<RealMatrix>{ms...});
}
 
#endif // MATRIX_H