mat.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_MATRIX_H
#define THEORETICA_MATRIX_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include <array>
#include <vector>
 
#include "../core/error.h"
#include "../core/constants.h"
#include "../core/real_analysis.h"
#include "./algebra.h"
#include "./transform.h"
#include "./vec.h"
 
 
namespace theoretica {
 
 
    template<typename Matrix, typename ReturnType = matrix_element_t<Matrix>&>
    class mat_iterator {
 
        private:
 
            Matrix& matrix;
 
            size_t row;
 
            size_t col;
 
        public:
 
            using iterator_category = std::forward_iterator_tag;
            using value_type = matrix_element_t<Matrix>;
            using pointer = value_type*;
            using reference = value_type&;
 
        
            mat_iterator(
                Matrix& matrix,
                size_t row = 0,
                size_t col = 0)
            : matrix(matrix), row(row), col(col) {}
            
            
            ReturnType operator*() {
                return matrix(row, col);
            }
 
        
            mat_iterator& operator++() {
 
                col++;
 
                if(col == matrix.cols()) {
                    col = 0;
                    row++;
                }
 
                return *this;
            }
 
 
            size_t row_index() {
                return row;
            }
 
 
            size_t col_index() {
                return col;
            }
 
 
            bool operator==(const mat_iterator& other) const {
                return (row == other.row) &&
                    (col == other.col);
            }
 
 
            bool operator!=(const mat_iterator& other) const {
                return !(*this == other);
            }
    };
 
 
    template<typename Type = real, unsigned int N = 0, unsigned int K = 0>
    class mat {
        public:
 
#ifdef THEORETICA_ROW_FIRST
        Type elements[N][K];
#else
        Type elements[K][N];
#endif
 
 
        mat() {
            algebra::mat_zeroes(*this);
        }
 
 
        template<typename Matrix, enable_matrix<Matrix> = true>
        mat(const Matrix& m) {
            algebra::mat_copy(*this, m);
        }
 
 
        template<typename T = Type>
        inline mat(const std::initializer_list<std::initializer_list<T>>& rows) {
 
            if(rows.size() != N) {
                TH_MATH_ERROR("mat::mat", rows.size(), MathError::InvalidArgument);
                algebra::mat_error(*this);
                return;
            }
 
            unsigned int i = 0;
            unsigned int j = 0;
 
            for (const auto& row : rows) {
 
                if (row.size() != K) {
                    TH_MATH_ERROR("mat::mat", rows.size(), MathError::InvalidArgument);
                    algebra::mat_error(*this);
                    return;
                }
 
                for (const auto& x : row) {
                    
                    get(i, j) = x;
                    j = (j + 1) % K;
                }
 
                i++;
            }
        }
 
 
        mat(Type diagonal, unsigned int n, unsigned int k) {
            
            if (n && k)
                resize(n, k);
 
            if (n != rows() && k != cols()) {
                TH_MATH_ERROR("mat::mat()", n, MathError::ImpossibleOperation);
                return;
            }
 
            algebra::mat_zeroes(*this);
            const unsigned int m = min(rows(), cols());
 
            for (unsigned int i = 0; i < m; ++i)
                get(i, i) = diagonal;
        }
 
 
        mat(Type diagonal) {
 
            algebra::mat_zeroes(*this);
            const unsigned int m = min(rows(), cols());
 
            for (unsigned int i = 0; i < m; ++i)
                get(i, i) = diagonal;
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        mat(const Vector& v, unsigned int row_size = N, unsigned int col_size = K) {
 
            if(row_size != N || col_size != K) {
                TH_MATH_ERROR("mat::mat", row_size * col_size, MathError::InvalidArgument);
                algebra::mat_error(*this);
                return;
            }
 
            if(v.size() != size()) {
                TH_MATH_ERROR("mat::mat", v.size(), MathError::InvalidArgument);
                algebra::mat_error(*this);
                return;
            }
 
            for (size_t i = 0; i < rows(); i++)
                for (size_t j = 0; j < cols(); j++)
                    get(i, j) = v[i * cols() + j];
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator=(const Matrix& other) {
            return algebra::mat_copy(*this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K> operator+(const Matrix& other) const {
            mat<Type, N, K> res;
            return algebra::mat_sum(res, *this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K> operator-(const Matrix& other) const {
            mat<Type, N, K> res;
            return algebra::mat_diff(res, *this, other);
        }
 
 
        inline mat<Type, N, K> operator*(Type scalar) const {
            mat<Type, N, K> res;
            return algebra::mat_scalmul(res, scalar, *this);
        }
 
 
        inline friend mat<Type, N, K> operator*(Type a, const mat<Type, N, K>& B) {
            return B * a;
        }
 
 
        template<typename VecType, unsigned int M>
        inline friend vec<VecType, K> operator*(
            const vec<VecType, M>& v, const mat<Type, N, K>& A) {
            return algebra::vec_mat_mul<vec<VecType, K>>(v, A);
        }
 
 
        inline mat<Type, N, K> operator/(Type scalar) const {
 
            mat<Type, N, K> res;
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, MathError::DivByZero);
                return algebra::mat_error(res);
            }
            
            return algebra::mat_scalmul(res, 1.0 / scalar, *this);
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Vector transform(const Vector& v) const {
 
            if(v.size() != cols()) {
                TH_MATH_ERROR("mat::transform", v.size(), MathError::InvalidArgument);
                Vector res;
                res.resize(cols());
                algebra::vec_error(res);
                return res;
            }
 
            return algebra::transform(*this, v);
        }
 
 
        inline vec<Type, N> transform(const vec<Type, K>& v) const {
 
            vec<Type, N> res;
            return algebra::transform(res, *this, v);
        }
 
 
        inline vec<Type, N> operator*(const vec<Type, K>& v) const {
            return transform(v);
        }
 
 
        template<unsigned int M>
        inline mat<Type, N, M> mul(const mat<Type, K, M>& B) const {
            mat<Type, N, M> res;
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix, enable_matrix<Matrix> = true>
        inline Matrix mul(const Matrix& B) const {
 
            Matrix res;
            res.resize(N, B.cols());
 
            if(B.rows() != K) {
                TH_MATH_ERROR("mat::transform", B.rows(), MathError::InvalidArgument);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix, enable_matrix<Matrix> = true>
        inline auto operator*(const Matrix& B) const {
 
            Matrix res;
            res.resize(N, B.cols());
 
            if(B.rows() != K) {
                TH_MATH_ERROR("mat::transform", B.rows(), MathError::InvalidArgument);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator+=(const Matrix& other) {
            return algebra::mat_sum(*this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator-=(const Matrix& other) {
            return algebra::mat_diff(*this, other);
        }
 
 
        inline mat<Type, N, K>& operator*=(Type scalar) {
            return algebra::mat_scalmul(scalar, *this);
        }
 
 
        inline mat<Type, N, K>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, MathError::DivByZero);
                return algebra::mat_error(*this);
            }
 
            return algebra::mat_scalmul(1.0 / scalar, *this);
        }
 
    
        template<typename Matrix>
        inline mat<Type, N, K>& operator*=(const Matrix& B) {
            return (*this = this->operator*(B));
        }
 
 
        inline mat<Type, N, K>& transpose() {
            static_assert(
                N == K, "The matrix must be square to be transposed in place.");
            return algebra::make_transposed(*this);
        }
 
 
        inline Type& get(unsigned int i, unsigned int j) {
 
#ifdef THEORETICA_ROW_FIRST
            return elements[i][j];
#else
            return elements[j][i];
#endif
        }
 
 
        inline const Type& get(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return elements[i][j];
#else
            return elements[j][i];
#endif
        }
 
    
        inline Type& at(unsigned int i, unsigned int j) {
 
            if (i >= rows()) {
                throw std::out_of_range(
                    "The provided row index in mat::at() is out of range"
                );
            }
 
            if (j >= cols()) {
                throw std::out_of_range(
                    "The provided column index in mat::at() is out of range"
                );
            }
 
            return get(i, j);
        }
 
 
        inline const Type& at(unsigned int i, unsigned int j) const {
 
            if (i >= rows()) {
                throw std::out_of_range(
                    "The provided row index in mat::at() is out of range"
                );
            }
 
            if (j >= cols()) {
                throw std::out_of_range(
                    "The provided column index in mat::at() is out of range"
                );
            }
 
            return get(i, j);
        }
 
 
        inline Type& operator()(unsigned int i, unsigned int j) {
            return get(i, j);
        }
 
 
        inline const Type& operator()(unsigned int i, unsigned int j) const {
            return get(i, j);
        }
 
 
        using iterator = mat_iterator<mat<Type, N, K>, Type&>;
 
 
        using const_iterator = mat_iterator<const mat<Type, N, K>, const Type&>;
 
 
        inline auto begin() {
            return iterator(*this, 0, 0);
        }
 
 
        inline auto end() {
            return iterator(*this, rows(), 0);
        }
 
    
        inline auto begin() const {
            return const_iterator(*this, 0, 0);
        }
 
    
        inline auto end() const {
            return const_iterator(*this, rows(), 0);
        }
 
    
        TH_CONSTEXPR inline unsigned int rows() const {
            return N;
        }
 
 
        TH_CONSTEXPR inline unsigned int cols() const {
            return K;
        }
 
 
        inline unsigned int size() const {
            return N * K;
        }
 
 
        inline Type* data() {
            return (Type*) elements;
        }
 
 
        inline const Type* data() const {
            return (const Type*) elements;
        }
 
 
        template<typename Vector = vec<Type, N * K>>
        inline Vector unpack() const {
            
            Vector res;
            res.resize(size());
 
            if (res.size() != size()) {
                TH_MATH_ERROR("mat::unpack", res.size(), MathError::InvalidArgument);
                return algebra::vec_error(res);
            }
 
            for (size_t i = 0; i < rows(); i++)
                for (size_t j = 0; j < cols(); j++)
                    res[i * cols() + j] = get(i, j);
 
            return res;
        }
 
 
        template<typename Matrix>
        inline bool operator==(const Matrix& other) const {
            return algebra::mat_equals(*this, other);
        }
 
 
        template<typename Matrix>
        inline bool operator!=(const Matrix& other) const {
            return !algebra::mat_equals(*this, other);
        }
 
 
        inline mat<Type, N, K>& invert() {
            static_assert(N == K, "The matrix must be square to be invertible.");
            return algebra::invert(*this);
        }
 
 
        inline mat<Type, N, K> resize(unsigned int n, unsigned int k) {
 
            if(rows() != n) {
                TH_MATH_ERROR("mat::resize", n, MathError::InvalidArgument);
                algebra::mat_error(*this);
            } else if(cols() != k) {
                TH_MATH_ERROR("mat::resize", k, MathError::InvalidArgument);
                algebra::mat_error(*this);
            }
 
            return *this;
        }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(
                std::string separator = ", ", bool parenthesis = true) const {
 
                std::stringstream res;
 
                for (unsigned int i = 0; i < rows(); ++i) {
                        
                    if(parenthesis)
                        res << "(";
 
                    for (unsigned int j = 0; j < cols(); ++j) {
                        
                        if(j)
                            res << separator;
 
                        if(abs(get(i, j)) < MACH_EPSILON)
                            res << "0";
                        else
                            res << get(i, j);
                    }
        
                    if(parenthesis)
                        res << ")" << std::endl;
                }
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(
                std::ostream& out, const mat<Type, N, K>& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
 
 
    template<typename Type>
    class mat<Type, 0, 0> {
        public:
 
        std::vector<Type> elements;
 
        size_t row_sz {0};
 
        size_t col_sz {0};
 
 
        mat() : row_sz(0), col_sz(0) {}
 
 
        mat(unsigned int n, unsigned int k) {
            resize(n, k);
            algebra::mat_zeroes(*this);
        }
 
    
        template <typename Matrix, enable_matrix<Matrix>>
        mat(const Matrix& m) {
            algebra::mat_copy(*this, m);
        }
 
    
        template<typename T = Type>
        inline mat(const std::initializer_list<std::initializer_list<T>>& rows) {
 
            unsigned int i = 0;
            unsigned int j = 0;
 
            for (const auto& row : rows) {
 
                for (const auto& x : row) {
                
                    if (!i && !j) {
                        this->resize(rows.size(), row.size());
                    }
                    else if (row.size() != col_sz) {
                        TH_MATH_ERROR("mat::mat", row.size(), MathError::InvalidArgument);
                        algebra::mat_error(*this);
                        return;
                    }
 
                    get(i, j) = x;
                    j = (j + 1) % col_sz;
                }
 
                i++;
            }
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator=(const Matrix& other) {
            resize(other.rows(), other.cols());
            return algebra::mat_copy(*this, other);
        }
 
 
        mat(Type diagonal, unsigned int n, unsigned int k) {
            
            resize(n, k);
            algebra::mat_zeroes(*this);
            const unsigned int m = min(n, k);
 
            for (unsigned int i = 0; i < m; ++i)
                get(i, i) = diagonal;
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        mat(const Vector& v, unsigned int row_size, unsigned int col_size) {
            
            resize(row_size, col_size);
 
            if(v.size() != size()) {
                TH_MATH_ERROR("mat::mat", v.size(), MathError::InvalidArgument);
                algebra::mat_error(*this);
                return;
            }
 
            for (size_t i = 0; i < rows(); i++)
                for (size_t j = 0; j < cols(); j++)
                    get(i, j) = v[i * cols() + j];
        }
 
 
        ~mat() {
            row_sz = 0;
            col_sz = 0;
        }
 
 
        inline void make_zeroes() {
            algebra::mat_zeroes(*this);
        }
 
 
        template<typename Matrix>
        inline mat<Type> operator+(const Matrix& other) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_sum(res, *this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type> operator-(const Matrix& other) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_diff(res, *this, other);
        }
 
 
        inline mat<Type> operator*(Type scalar) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_scalmul(res, scalar, *this);
        }
 
 
        inline friend mat<Type> operator*(Type a, const mat<Type>& B) {
            return B * a;
        }
 
 
        template<typename VecType, unsigned int M>
        inline friend vec<VecType> operator*(
            const vec<VecType, M>& v, const mat<Type, 0, 0>& A) {
            return algebra::vec_mat_mul<vec<VecType>>(v, A);
        }
 
 
        inline mat<Type> operator/(Type scalar) const {
 
            mat<Type> res;
            res.resize(rows(), cols());
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, MathError::DivByZero);
                return algebra::mat_error(res);
            }
            
            return algebra::mat_scalmul(res, 1.0 / scalar, *this);
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Vector transform(const Vector& v) const {
 
            if(v.size() != rows()) {
                TH_MATH_ERROR("mat::transform", v.size(), MathError::InvalidArgument);
                Vector res;
                res.resize(rows());
                return algebra::vec_error(res);
            }
 
            return algebra::transform(*this, v);
        }
 
 
        template<unsigned int N = 0, unsigned int K = 0>
        inline vec<Type, N> transform(const vec<Type, K>& v) const {
 
            vec<Type, N> res;
            res.resize(rows());
            
            return algebra::transform(res, *this, v);
        }
 
 
        template<unsigned int N = 0, unsigned int K = 0>
        inline vec<Type, N> operator*(const vec<Type, K>& v) const {
            return transform(v);
        }
 
 
        // inline vec<Type> operator*(const vec<Type>& v) const {
        //  return transform(v);
        // }
 
 
        inline mat<Type> mul(const mat<Type>& B) const {
 
            mat<Type> res;
            res.resize(rows(), B.cols());
 
            if(B.rows() != cols()) {
                TH_MATH_ERROR("mat::mul", B.rows(), MathError::InvalidArgument);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix, enable_matrix<Matrix> = true>
        inline Matrix mul(const Matrix& B) const {
 
            Matrix res;
            res.resize(rows(), B.cols());
 
            if(B.rows() != cols()) {
                TH_MATH_ERROR("mat::mul", B.rows(), MathError::InvalidArgument);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline auto operator*(const Matrix& B) const {
            return mul(B);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator+=(const Matrix& other) {
            return algebra::mat_sum(*this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator-=(const Matrix& other) {
            return algebra::mat_diff(*this, other);
        }
 
 
        inline mat<Type>& operator*=(Type scalar) {
            return algebra::mat_scalmul(scalar, *this);
        }
 
 
        inline mat<Type>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, MathError::DivByZero);
                return algebra::mat_error(*this);
            }
 
            return algebra::mat_scalmul(1.0 / scalar, *this);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator*=(const Matrix& B) {
            return (*this = this->operator*(B));
        }
 
 
        inline mat<Type>& transpose() {
 
            mat<Type> res;
            res.resize(cols(), rows());
 
            for (size_t i = 0; i < rows(); i++)
                for (size_t j = 0; j < cols(); j++)
                    res(j, i) = get(i, j);
 
            return (*this = res);
        }
 
 
        inline Type& get(unsigned int i, unsigned int j) {
 
#ifdef THEORETICA_ROW_FIRST
            return elements[j + i * col_sz];
#else
            return elements[i + j * row_sz];
#endif
        }
 
 
        inline const Type& get(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return elements[j + i * col_sz];
#else
            return elements[i + j * row_sz];
#endif
        }
 
 
        inline Type& at(unsigned int i, unsigned int j) {
 
            if (i >= rows()) {
                throw std::out_of_range(
                    "The provided row index in mat::at() is out of range"
                );
            }
 
            if (j >= cols()) {
                throw std::out_of_range(
                    "The provided column index in mat::at() is out of range"
                );
            }
 
            return get(i, j);
        }
 
 
        inline const Type& at(unsigned int i, unsigned int j) const {
 
            if (i >= rows()) {
                throw std::out_of_range(
                    "The provided row index in mat::at() is out of range"
                );
            }
 
            if (j >= cols()) {
                throw std::out_of_range(
                    "The provided column index in mat::at() is out of range"
                );
            }
 
            return get(i, j);
        }
 
 
        inline Type& operator()(unsigned int i, unsigned int j) {
            return get(i, j);
        }
 
 
        inline const Type& operator()(unsigned int i, unsigned int j) const {
            return get(i, j);
        }
 
 
        using iterator = mat_iterator<mat<Type, 0, 0>, Type&>;
 
 
        using const_iterator = mat_iterator<const mat<Type, 0, 0>, const Type&>;
 
 
        inline auto begin() {
            return iterator(*this, 0, 0);
        }
 
 
        inline auto end() {
            return iterator(*this, rows(), 0);
        }
 
 
        inline auto begin() const {
            return const_iterator(*this, 0, 0);
        }
 
 
        inline auto end() const {
            return const_iterator(*this, rows(), 0);
        }
 
 
        TH_CONSTEXPR inline unsigned int rows() const {
            return row_sz;
        }
 
 
        TH_CONSTEXPR inline unsigned int cols() const {
            return col_sz;
        }
 
 
        inline unsigned int size() const {
            return rows() * cols();
        }
 
 
        inline Type* data() {
            return elements.data();
        }
 
 
        inline const Type* data() const {
            return elements.data();
        }
 
 
        template<typename Vector = vec<Type>>
        inline Vector unpack() const {
            
            Vector res;
            res.resize(size());
 
            if (res.size() != size()) {
                TH_MATH_ERROR("mat::unpack", res.size(), MathError::InvalidArgument);
                return algebra::vec_error(res);
            }
 
            for (size_t i = 0; i < rows(); i++)
                for (size_t j = 0; j < cols(); j++)
                    res[i * cols() + j] = get(i, j);
 
            return res;
        }
 
 
        template<typename Matrix>
        inline bool operator==(const Matrix& other) const {
            return algebra::mat_equals(*this, other);
        }
 
 
        template<typename Matrix>
        inline bool operator!=(const Matrix& other) const {
            return !algebra::mat_equals(*this, other);
        }
 
 
        inline mat<Type>& invert() {
            return algebra::invert(*this);
        }
 
 
        inline mat<Type>& resize(unsigned int rows, unsigned int cols) {
 
            // Do nothing if the size is already correct
            if (row_sz == rows && col_sz == cols)
                return *this;
 
            std::vector<Type> new_data (uint64_t(rows) * cols);
 
            if (elements.size()) {
 
                size_t min_elements = min(rows, row_sz) * min(cols, col_sz);
 
                // Copy the overlapping elements
                for (unsigned int i = 0; i < min_elements; ++i)
                    new_data[i] = elements[i];
            }
 
            elements = new_data;
            row_sz = rows;
            col_sz = cols;
 
            return *this;
        }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(
                std::string separator = ", ", bool parenthesis = true) const {
 
                std::stringstream res;
 
                for (unsigned int i = 0; i < rows(); ++i) {
                        
                    if(parenthesis)
                        res << "(";
 
                    for (unsigned int j = 0; j < cols(); ++j) {
                        
                        if(j)
                            res << separator;
 
                        if(abs(get(i, j)) < MACH_EPSILON)
                            res << "0";
                        else
                            res << get(i, j);
                    }
        
                    if(parenthesis)
                        res << ")" << std::endl;
                }
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(
                std::ostream& out, const mat<Type>& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
}
 
#endif