matrix.hpp

Go to the documentation of this file.

/*
 *  Cypress -- C++ Neural Associative Memory Simulator
 *  Copyright (C) 2016  Andreas Stöckel
 *
 *  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/>.
 */

#ifndef CYPRESS_UTIL_MATRIX_HPP
#define CYPRESS_UTIL_MATRIX_HPP

#include <algorithm>
#include <cassert>
#include <memory>
#include <ostream>
#include <vector>
#ifndef NDEBUG
#include <sstream>
#include <stdexcept>
#endif

namespace cypress {

enum class MatrixFlags { NONE = 0, ZEROS = 1 };

template <typename T>
class Matrix {
private:
    T *m_buf;

    size_t m_rows, m_cols;

    bool m_destroy = true;

#ifndef NDEBUG

    void check_range(size_t i) const
    {
        if (i > size()) {
            std::stringstream ss;
            ss << "[" << i << "] out of range for matrix of size " << size()
               << std::endl;
            throw std::out_of_range(ss.str());
        }
    }

    void check_range(size_t row, size_t col) const
    {
        if (row >= m_rows || col >= m_cols) {
            std::stringstream ss;
            ss << "[" << row << ", " << col
               << "] out of range for matrix of size " << m_rows << " x "
               << m_cols << std::endl;
            throw std::out_of_range(ss.str());
        }
    }
#else

    void check_range(size_t) const {}
    void check_range(size_t, size_t) const {}
#endif

public:
    Matrix() : m_buf(nullptr), m_rows(0), m_cols(0) {}

    Matrix(std::initializer_list<T> init) : Matrix(init.size(), 1)
    {
        size_t i = 0;
        for (const auto &elem : init) {
            (*this)[i++] = elem;
        }
    }

    template <size_t Rows, size_t Cols>
    Matrix(const std::array<std::array<T, Cols>, Rows> &init)
        : Matrix(Rows, Cols)
    {
        for (size_t i = 0; i < Rows; i++) {
            for (size_t j = 0; j < Cols; j++) {
                (*this)(i, j) = init[i][j];
            }
        }
    }

    Matrix(const std::vector<std::vector<T>> &init)
        : Matrix(init.size(), init[0].size())
    {
        for (size_t i = 1; i < init.size(); i++) {
            assert(init[0].size() == init[i].size());
        }
        for (size_t i = 0; i < init.size(); i++) {
            for (size_t j = 0; j < init[0].size(); j++) {
                (*this)(i, j) = init[i][j];
            }
        }
    }

    Matrix(const std::vector<T> &init) : Matrix(init.size(), 1)
    {
        for (size_t i = 0; i < init.size(); i++) {
            (*this)(i, 0) = init[i];
        }
    }

    Matrix(size_t rows, size_t cols, MatrixFlags flags = MatrixFlags::NONE)
        : m_buf(new T[rows * cols]), m_rows(rows), m_cols(cols)
    {
        if (flags == MatrixFlags::ZEROS) {
            fill(T());
        }
    }

    Matrix(size_t rows, size_t cols, const T *data)
        : m_buf(new T[rows * cols]), m_rows(rows), m_cols(cols)
    {
        std::copy(data, data + (rows * cols), begin());
    }

    Matrix(size_t rows, size_t cols, T *data, bool destroy)
        : m_buf(data), m_rows(rows), m_cols(cols), m_destroy(destroy){};

    Matrix(const Matrix &o)
        : m_buf(new T[o.rows() * o.cols()]), m_rows(o.rows()), m_cols(o.cols())
    {
        std::copy(o.begin(), o.end(), begin());
    }

    Matrix(Matrix &&o) noexcept
        : m_buf(o.m_buf), m_rows(o.rows()), m_cols(o.cols())
    {
        o.m_buf = nullptr;
        o.m_rows = 0;
        o.m_cols = 0;
    }

    Matrix &operator=(const Matrix &o)
    {
        delete[] m_buf;
        m_rows = o.m_rows;
        m_cols = o.m_cols;
        m_buf = new T[o.rows() * o.cols()];
        std::copy(o.begin(), o.end(), begin());
        return *this;
    }

    Matrix &operator=(Matrix &&o) noexcept
    {
        delete[] m_buf;
        m_rows = o.m_rows;
        m_cols = o.m_cols;
        m_buf = o.m_buf;
        o.m_buf = nullptr;
        o.m_rows = 0;
        o.m_cols = 0;
        return *this;
    }

    ~Matrix()
    {
        if (m_destroy) {
            delete[] m_buf;
        }
    }

    operator std::vector<T>() const { return std::vector<T>(begin(), end()); }

    bool operator==(const Matrix<T> &o) const
    {
        return rows() == o.rows() && cols() == o.cols() &&
               std::equal(begin(), end(), o.begin());
    }

    Matrix<T> &fill(const T &val)
    {
        std::fill(begin(), end(), val);
        return *this;
    }

    T *begin(size_t row = 0) { return m_buf + row * cols(); }
    T *end() { return m_buf + size(); }
    const T *begin(size_t row = 0) const { return m_buf + row * cols(); }
    const T *end() const { return m_buf + size(); }
    const T *cbegin(size_t row = 0) const { return m_buf + row * cols(); }
    const T *cend() const { return m_buf + size(); }

    T &operator()(size_t row, size_t col)
    {
        check_range(row, col);
        return *(m_buf + row * m_cols + col);
    }

    const T &operator()(size_t row, size_t col) const
    {
        check_range(row, col);
        return *(m_buf + row * m_cols + col);
    }

    T &operator()(size_t i)
    {
        check_range(i);
        return *(m_buf + i);
    }

    const T &operator()(size_t i) const
    {
        check_range(i);
        return *(m_buf + i);
    }

    T &operator[](size_t i)
    {
        check_range(i);
        return *(m_buf + i);
    }

    const T &operator[](size_t i) const
    {
        check_range(i);
        return *(m_buf + i);
    }

    size_t rows() const { return m_rows; }

    size_t cols() const { return m_cols; }

    size_t size() const { return m_rows * m_cols; }

    void resize(size_t rows, size_t cols)
    {
        if (rows != m_rows || cols != m_cols) {
            delete[] m_buf;
            m_buf = new T[rows * cols];
            m_rows = rows;
            m_cols = cols;
        }
    }

    void reshape(size_t rows, size_t cols)
    {
        if (rows * cols == m_rows * m_cols) {
            m_rows = rows;
            m_cols = cols;
        }
        else {
            resize(rows, cols);
        }
    }

    T *data() { return m_buf; }

    const T *data() const { return m_buf; }

    friend std::ostream &operator<<(std::ostream &os, const Matrix<T> &m)
    {
        T const *d = m.data();
        for (size_t row = 0; row < m.rows(); row++) {
            for (size_t col = 0; col < m.cols(); col++) {
                os << (col == 0 ? "" : ",") << *d;
                d++;
            }
            os << std::endl;
        }
        return os;
    }

    bool empty() const { return size() == 0; }

    Matrix submatrix(size_t row, size_t col) const
    {
        if (m_rows <= 1 || m_cols <= 1) {
            throw std::invalid_argument(
                "For calculation of a submatrix the dimension must be at least "
                "2x2!");
        }
        if (m_rows <= row || m_cols <= col) {
            throw std::invalid_argument(
                "Matrix too small for calculating the submatrix with given row "
                "and column!");
        }
        Matrix<T> subm(m_rows - 1, m_cols - 1);

        size_t subi = 0;
        for (size_t i = 0; i < m_rows; i++) {
            if (i == row) {
                continue;
            }
            int subj = 0;
            for (size_t j = 0; j < m_cols; j++) {
                if (j == col) {
                    continue;
                }
                subm(subi, subj) = *(m_buf + i * m_cols + j);
                subj++;
            }
            subi++;
        }
        return subm;
    }

    T determinant() const
    {
        if (m_rows * m_cols == 0 || m_rows != m_cols) {
            throw std::invalid_argument(
                "Calculation of the determinant only valid for squared "
                "matrices with dimension >=1!");
        }
        size_t dim = m_rows;
        if (dim == 1) {
            return *(m_buf);
        }
        else if (dim == 2) {
            return (*(m_buf)) * (*(m_buf + 3)) -
                   (*(m_buf + 1)) * (*(m_buf + 2));
        }
        Matrix<T> subm(dim - 1, dim - 1);
        T res = 0;
        for (size_t x = 0; x < dim; x++) {
            T tmp = (*(m_buf + x));  // mat(0.x)
            if (tmp) {
                T sign = x % 2 ? -1 : 1;
                res += sign * tmp * submatrix(0, x).determinant();
            }
        }
        return res;
    }

    Matrix adjugate() const
    {
        if (m_rows * m_cols == 0 || m_rows != m_cols) {
            throw std::invalid_argument(
                "Calculation of the determinant only valid for squared "
                "matrices with dimension >=1!");
        }
        size_t dim = m_rows;
        if (dim == 1) {
            return Matrix(*this);
        }
        Matrix<T> res(m_rows, m_cols);
        for (size_t x = 0; x < dim; x++) {
            for (size_t y = 0; y < dim; y++) {
                T sign = (x + y) % 2 ? -1 : 1;
                res(x, y) = sign * submatrix(y, x).determinant();
            }
        }
        return res;
    }

    Matrix<double> inverse() const
    {
        auto det = determinant();
        if (det == 0) {
            throw std::invalid_argument(
                "Could not calculate the inverse, as the determinant is zero!");
        }
        auto adj = adjugate();
        Matrix<double> res(adj.rows(), adj.cols());
        for (size_t i = 0; i < res.size(); i++) {
            res[i] = double(adj[i]) / double(det);
        }
        return res;
    }
};

template <typename T, size_t Rows, size_t Cols>
Matrix<T> make_matrix(const std::array<std::array<T, Cols>, Rows> &init)
{
    return Matrix<T>(init);
}

template <typename T>
Matrix<T> make_matrix(std::initializer_list<T> list)
{
    return Matrix<T>(list);
}

template <typename T>
class Vector : public Matrix<T> {
public:
    Vector() : Matrix<T>(0, 0) {}
    Vector(size_t s, MatrixFlags flags = MatrixFlags::NONE)
        : Matrix<T>(s, 1, flags)
    {
    }
    Vector(std::initializer_list<T> init) : Vector(init.size())
    {
        size_t i = 0;
        for (const auto &elem : init) {
            (*this)[i++] = elem;
        }
    }
    void resize(size_t s) { resize(s, 1); }
};

}  // namespace cypress

#endif /* CYPRESS_UTIL_MATRIX_HPP */