New matrix data model (breaks everything)

1 files changed, 66 insertions, 529 deletions

diff --git a/include/mm/mmmatrix.hpp b/include/mm/mmmatrix.hpp
index facbb48..d653d72 100644
--- a/include/mm/mmmatrix.hpp
+++ b/include/mm/mmmatrix.hpp
@@ -12,580 +12,117 @@
 #pragma once
 
 #include <iostream>
-#include <iomanip>
-#include <cstring>
 #include <cassert>
 #include <initializer_list>
 #include <array>
 #include <memory>
 
-#include "mm/mmiterator.hpp"
 
 namespace mm {
-
-    /* basic grid structure */
+    using index = std::size_t;
 
     template<typename T, std::size_t Rows, std::size_t Cols>
     class basic_matrix;
 
-    /* basic wrapper */
-
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    class matrix; // simple matrix format
-    
     /* specialisations */
-
-    /* specialization of a matrix */
-    template<typename T, std::size_t N>
-    class vector; // by default, set a column vector
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class matrix;
 
     template<typename T, std::size_t N>
     class square_matrix;
 
-    /* specialisation of a square_matrix for a sub-diagonal composed matrix */
-    template<typename T, std::size_t N, signed long ... Diags>
-    class multi_diag_matrix;
+    template<typename T, std::size_t N>
+    class diagonal_matrix;
 }
 
 /*
  * Matrix class, no access methods
  */
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::basic_matrix
-{
-public:
-    using type = T;
-
-    template<typename U, std::size_t ORows, std::size_t OCols>
-    friend class mm::basic_matrix;
-
-    template<typename U, std::size_t ORows, std::size_t OCols>
-    friend class mm::matrix;
-
-    template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
-    friend class mm::iter::basic_iterator;
-
-    template<typename U, std::size_t ON, class Grid>
-    friend class mm::iter::diag_iterator;
-
-    //template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
-    //friend class mm::iter::basic_iterator<T, Rows, Cols, mm::basic_matrix<T, Rows, Cols>>;
-
-    //template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
-    //friend class mm::iter::basic_iterator<typename std::add_const<T>::type, Rows, Cols, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
-
-    basic_matrix();
-
-    // from initializer_list
-    basic_matrix(std::initializer_list<std::initializer_list<T>> l);
-
-    // copyable and movable
-    basic_matrix(const basic_matrix<T, Rows, Cols>& other) = default;
-    basic_matrix(basic_matrix<T, Rows, Cols>&& other) = default;
-
-    // copy from another matrix
-    template<std::size_t ORows, std::size_t OCols>
-    basic_matrix(const basic_matrix<T, ORows, OCols>& other);
-
-    void swap_rows(std::size_t x, std::size_t y);
-    void swap_cols(std::size_t x, std::size_t y);
-
-    // mathematical operations
-    //virtual basic_matrix<T, Cols, Rows> transposed() const;
-    //inline basic_matrix<T, Cols, Rows> td() const { return transposed(); }
-
-protected:
-    template<typename ConstIterator>
-    basic_matrix(ConstIterator begin, ConstIterator end);
-
-private:
-    std::array<T, Rows * Cols> data;
-};
-
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix() {
-    std::fill(data.begin(), data.end(), 0);
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
-    std::initializer_list<std::initializer_list<T>> l
-) {
-    assert(l.size() == Rows);
-    auto data_it = data.begin();
-
-    for (auto&& row : l) {
-        data_it = std::copy(row.begin(), row.end(), data_it);
-    }
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-template<std::size_t ORows, std::size_t OCols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
-    const mm::basic_matrix<T, ORows, OCols>& other
-) {
-    static_assert((ORows <= Rows),
-        "cannot copy a taller matrix into a smaller one"
-    );
-
-    static_assert((OCols <= Cols),
-        "cannot copy a larger matrix into a smaller one"
-    );
-
-    std::fill(data.begin(), data.end(), 0);
-    for (unsigned row = 0; row < Rows; row++)
-        for (unsigned col = 0; col < Cols; col++)
-            this->at(row, col) = other.at(row, col);
-}
-
-/* protected construtor */
-template<typename T, std::size_t Rows, std::size_t Cols>
-template<typename ConstIterator>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
-    ConstIterator begin, ConstIterator end
-) {
-    assert(static_cast<unsigned>(std::distance(begin, end)) >= ((Rows * Cols)));
-    std::copy(begin, end, data.begin());
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-void mm::basic_matrix<T, Rows, Cols>::swap_rows(std::size_t x, std::size_t y) {
-    if (x == y)
-        return;
-
-    for (unsigned col = 0; col < Cols; col++)
-        std::swap(this->at(x, col), this->at(y, col));
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-void mm::basic_matrix<T, Rows, Cols>::swap_cols(std::size_t x, std::size_t y) {
-    if (x == y)
-        return;
-
-    for (unsigned row = 0; row < Rows; row++)
-        std::swap(this->at(row, x), this->at(row, y));
-}
-
-/*
- * Matrix object
- */
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::matrix
-{
-public:
-
-    //template<typename U, std::size_t ORows, std::size_t OCols>
-    using vec_iterator = mm::iter::basic_iterator<T, Rows, Cols, mm::basic_matrix<T, Rows, Cols>>;
-
-    //template<typename U, std::size_t ORows, std::size_t OCols>
-    using const_vec_iterator = mm::iter::basic_iterator<typename std::add_const<T>::type, Rows, Cols, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
-    
-    // default zeros constructor
-    matrix() : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>()), transposed(false) {}
-    
-    // from initializer_list
-    matrix(std::initializer_list<std::initializer_list<T>> l)
-        : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>(l)), transposed(false) {}
-
-    // copyable and movable
-    matrix(const matrix<T, Rows, Cols>& other) // deep copy
-        : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>(*other.M)), transposed(other.transposed) {}
-
-    matrix(basic_matrix<T, Rows, Cols>&& other) // move ptr
-        : M(other.M), transposed(other.transposed)
-    {
-        other.M = nullptr;
-    }
-
-    matrix<T, Rows, Cols> operator=(const basic_matrix<T, Rows, Cols>& other) // deep copy
-    {         
-        *M = *other.M;
-        transposed = other.transposed;
-    }
-
-    /* 
-     * Transposition
-     */
-
-    matrix<T, Rows, Cols>& transpose_d()
-    {
-        transposed = !transposed;
-        return *this;
-    }
-
-    const matrix<T, Rows, Cols> transpose() const
-    {
-        return matrix<T, Rows, Cols>(M, !transposed);
-    }
-
-    inline matrix<T, Rows, Cols>& td()
-    {
-        return transpose();
-    }
-
-    inline matrix<T, Rows, Cols> t() const
-    {
-        return transpose();
-    }
-
-    // strongly transpose
-    matrix<T, Cols, Rows> transpose_cpy() const
-    {
-        matrix<T, Cols, Rows> out(); // copy
-        // TODO
-    }
-
-    /*
-     * Pointer status 
-     */
-
-    bool expired() const
+namespace mm {
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class basic_matrix
     {
-        return M == nullptr;
-    }
-
-    /*
-     * Downcasting conditions
-     */
-
-    /// downcast to square matrix
-    static inline constexpr bool is_square() { return (Rows == Cols); }
-    inline constexpr square_matrix<T, Rows> to_square() const {
-        static_assert(is_square());
-        return static_cast<square_matrix<T, Rows>>(*this);
-    }
-
-    /// downcast to col_vector
-    static inline constexpr bool is_vector() { return (Rows == 1 || Cols == 1); }
-    inline vector<T, Cols> to_vector() const {
-        if constexpr(Cols == 1)
-            return static_cast<vector<T, Rows>>(*this);
-        else if (Rows == 1)
-            return vector<T, Cols>(*this); // copy into column vector
+    public:
+        using type = T;
 
-    }
+        template<typename U, std::size_t ORows, std::size_t OCols>
+        friend class mm::matrix;
 
-    /* Accessors */
+        // copy from another matrix
+        template<std::size_t ORows, std::size_t OCols>
+        matrix(const basic_matrix<T, ORows, OCols>& other);
 
-    virtual T& at(std::size_t row, std::size_t col)
-    {
-        return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
-    }
+        virtual T& at(index row, index col) = 0;
+        virtual const T& at(index row, index col) const = 0;
+    };
 
-    virtual const T& at(std::size_t row, std::size_t col) const
-    {
-        return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
-    }
 
-    std::size_t rows() const {
-        return (transposed) ? Cols : Rows;
-    }
-
-    std::size_t cols() const {
-        return (transposed) ? Rows : Cols;
-    }
-
-    virtual mm::matrix<T, Rows, Cols>::vec_iterator operator[](std::size_t index)
-    {
-        return mm::matrix<T, Rows, Cols>::vec_iterator(*M, index, 0, !transposed);
-    }
+    /* Specializations */
 
-    virtual mm::matrix<T, Rows, Cols>::const_vec_iterator operator[](std::size_t index) const
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    struct matrix : public basic_matrix<T, N>
     {
-        return mm::matrix<T, Rows, Cols>::const_vec_iterator(*M, index, 0, !transposed);
-    }
-
-    /*
-     * Basic matematical operations (dimension indipendent)
-     */
-
-    mm::matrix<T, Rows, Cols>& operator+=(const mm::matrix<T, Rows, Cols>& m) {
-
-        for (unsigned row = 0; row < std::min(rows(), m.rows()); ++row)
-            for (unsigned col = 0; col < std::min(cols(), m.cols()); ++col)
-                at(row, col) += m.at(row, col);
-    
-        return *this;
-    }
-
-    mm::matrix<T, Rows, Cols>& operator-=(const mm::matrix<T, Rows, Cols>& m) {
-
-        for (unsigned row = 0; row < std::min(rows(), m.rows()); ++row)
-            for (unsigned col = 0; col < std::min(cols(), m.cols()); ++col)
-                at(row, col) -= m.at(row, col);
-    
-        return *this;
-    }
-
-    mm::matrix<T, Rows, Cols> operator*=(const T& k) {
-
-        for (unsigned row = 0; row < rows(); ++row)
-            for (auto& x : (*this)[row])
-                x *= k;
-    
-        return *this;
-    }
-
-protected:
-
-    std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> M;
-
-    // shallow construction
-    matrix(std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> grid, bool tr = false) : M(grid), transposed(tr) {}
-    
-private:
-
-    bool transposed;     
-};
-
-/* Basic operator overloading (dimension indipendent) */
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator+(
-    mm::matrix<T, Rows, Cols> a,
-    const mm::matrix<T, Rows, Cols>& b
-) {
-    return a += b;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator-(
-    mm::matrix<T, Rows, Cols> a,
-    const mm::matrix<T, Rows, Cols>& b
-) {
-    return a -= b;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator*(
-    mm::matrix<T, Rows, Cols> a,
-    const T& k
-) {
-    return a *= k;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator*(
-    const T& k,
-    mm::matrix<T, Rows, Cols> a
-) {
-    return a *= k;
-}
-
-// simple multiplication
-template<typename T, std::size_t M, std::size_t P1, std::size_t P2, std::size_t N>
-mm::matrix<T, M, N> operator*(
-    const mm::matrix<T, M, P1>& a,
-    const mm::matrix<T, P2, N>& b
-) {
-    // TODO, adjust asserts for transposed cases
-    static_assert(P1 == P2, "invalid matrix multiplication");
-    assert(a.cols() == b.rows());
-
-    mm::matrix<T, M, N> result;
-    const mm::matrix<T, P2, N> bt = b.t(); // weak transposition
-
-    //npdebug("Calling *")
-
-    for (unsigned row = 0; row < M; row++)
-        for (unsigned col = 0; col < N; col++)
-            result.at(row, col) = a[row] * bt[col]; // scalar product
-
-    return result;
-}
-
-/*
- * Matrix operator <<
- */
-
-template<typename T, std::size_t Rows, std::size_t Cols, unsigned NumW = 3>
-std::ostream& operator<<(std::ostream& os, const mm::matrix<T, Rows, Cols>& m) {
-
-    for (unsigned index = 0; index < m.rows(); index++) {
-        os << "[ ";
-        for (unsigned col = 0; col < m.cols()-1; ++col) {
-            os << std::setw(NumW) << m.at(index, col) << ", ";
+    public:
+        virtual T& at(index row, index col) override {
+            return m_data[row * Cols + col];
         }
-        os << std::setw(NumW) << m.at(index, m.cols()-1) << " ]\n";
-    }
-
-    return os;
-}
-
-/*
- * Vector, TODO better manage column and row
- */
-
-template<typename T, std::size_t N>
-class mm::vector : public mm::matrix<T, N, 1>
-{
-public:
-
-    using mm::matrix<T, N, 1>::matrix;
-
-    vector(std::initializer_list<T> l)
-        : mm::matrix<T, N, 1>(l) {}
-};
 
-template<typename T>
-mm::vector<T, 3> operator^(const mm::vector<T, 3>& v, const mm::vector<T, 3>& w)
-{
-    mm::vector<T, 3> out;
-
-    out[0] = v[1] * w[2] - v[2] * w[2];
-    out[1] = v[2] * w[0] - v[0] * w[2];
-    out[2] = v[0] * w[1] - v[1] * w[0];
-
-    return out;
-}
-
-/*
- * Square matrix
- */
-
-template<typename T, std::size_t N>
-class mm::square_matrix : public mm::matrix<T, N, N>
-{
-public:
-
-    using mm::matrix<T, N, N>::matrix;
-
-    using diag_iterator = mm::iter::diag_iterator<T, N, mm::basic_matrix<T, N, N>>;
-
-    using const_diag_iterator = mm::iter::diag_iterator<typename std::add_const<T>::type, N, typename std::add_const<mm::basic_matrix<T, N, N>>::type>;
-
-    virtual T trace();
-    inline T tr() { return trace(); }
-
-    virtual mm::square_matrix<T, N>::diag_iterator diag_beg(int row = 0)
-    {
-        return diag_iterator(*(this->M), row, 0);
-    }
-
-    virtual mm::square_matrix<T, N>::const_diag_iterator diag_end(int row = 0) const
-    {
-        return const_diag_iterator(*(this->M), row, N);
-    }
-
-    // TODO, determinant
-
-    /// in place inverse
-    // TODO, det != 0
-    // TODO, use gauss jordan for invertible ones
-    //void invert();, TODO, section algorithm
-    
-    /*
-     * Generate the identity
-     */
-
-    static inline constexpr mm::square_matrix<T, N> identity() {
-        mm::square_matrix<T, N> i;
-        for (unsigned row = 0; row < N; row++)
-            for (unsigned col = 0; col < N; col++)
-                i.at(row, col) = (row == col) ? 1 : 0;
-
-        return i;
-    }
-};
-
-template<typename T, std::size_t N>
-T mm::square_matrix<T, N>::trace() 
-{
-    T sum = 0;
-    for (const auto& x : diag_beg())
-        sum += x;
+        virtual const T& at(index row, index col) const override {
+            return at(row, col);
+        }
 
-    return sum;
-}
+    private:
+        std::array<T, Rows * Cols> m_data;
+    };
 
-// TODO, static assert, for all: Diags > -N, Diags < N
-// TODO, force Diags to be ordered
-template<typename T, std::size_t N, signed long ... Diags>
-class mm::multi_diag_matrix : public mm::square_matrix<T, N>
-{
-    T& shared_zero = 0;
 
-public:
-    using mm::square_matrix<T, N>::square_matrix;
+    template<typename T, std::size_t N>
+    struct vector : public matrix<T, 1, N> {};
 
-    // TODO, ordered case: dichotomy search O(log(M))
-    // M = parameter pack size
-    static inline bool constexpr is_in(std::size_t i, std::size_t j)
+    template<typename T, std::size_t N>
+    struct square_matrix : public basic_matrix<T, N>
     {
-        auto t = std::make_tuple(Diags...);
+    public:
+        virtual T& at(index row, index col) override {
+            return m_data[row * N + col];
+        }
 
-        for(unsigned k(0); k < sizeof...(Diags); ++k)
-            if ((i - j) == std::get<k>(t))
-                return true;
+        virtual const T& at(index row, index col) const override {
+            return at(row, col);
+        }
 
-        return false;
-    }
+    private:
+        std::array<T, N*N> m_data;
+    };
 
-    virtual T& at(std::size_t row, std::size_t col) override
+    template<typename T, std::size_t N>
+    struct identity_matrix : public basic_matrix<T, N, N>
     {
-        if (is_in(row, col))
-            return mm::square_matrix<T, N>::at(row, col);
+    public:
+        const T& at(index row, index col) const override {
+            return (row != col) ? static_cast<T>(1) : static_cast<T>(0);
+        }
 
-        shared_zero = 0;
-        return shared_zero;
+    private:
+        T m_useless;
+        T& at(index row, index col) { return m_useless; }
     }
 
-    virtual const T& at(std::size_t row, std::size_t col) const override
+    template<typename T, std::size_t N>
+    struct diagonal_matrix : public basic_matrix<T, N, N>
     {
-        if (is_in(row, col))
-            return mm::square_matrix<T, N>::at(row, col);
-
-        return 0;
-    }
-
-    
-
-    // TODO, implement limited iterators
-};
+    public:
+        T& at(index row, index col) override {
+            n_null_element = static_cast<T>(0);
+            return (row != col) ? m_data[row] : n_null_element;
+        }
 
-/*template<typename T, std::size_t N, signed long ... Diags, std::size_t P, std::size_t M>
-void constexpr diag_mult(const mm::multi_diag_matrix<T, N, Diags...>& a,
-    const mm::matrix<T, P, M>& b, mm::matrix<T, M, N>& result)
-{
-    static_assert(N == P && N == M, "invalid diagonal multiplication");
+        const T& at(index row, index col) const override {
+            return (row != col) ? m_data[row] : static_cast<T>(0);
+        }
 
-    auto d = a.diagonal(Diags);
-    if constexpr (Diags < 0) {
-        for (unsigned k = 0; k < M; ++k) 
-            for (unsigned i = -Diags; i < N; ++i)
-                result.at(i + Diags, k) += d[i + Diags] * b.at(i, k);
-    } else {
-        for (unsigned k = 0; k < M; ++k) 
-            for (unsigned i = Diags; i < N; ++i)
-                result.at(i, k) += d[i - Diags] * b.at(i - Diags, k);
+    private:
+        T m_null_element;
+        std::array<T, N> m_data;
     }
 }
-
-template<typename T, std::size_t N, signed long ... Diags, std::size_t P, std::size_t M>
-mm::matrix<T, M, N> operator*(
-    const mm::multi_diag_matrix<T, N, Diags...>& a,
-    const mm::matrix<T, P, M>& b
-) {
-    static_assert(N == P && N == M, "invalid matrix multiplication");
-    assert(a.cols() == b.rows());
-
-    mm::matrix<T, M, N> result;
-    ((
-    auto d = a.diagonal(Diags);
-    if constexpr (Diags < 0) {
-        for (unsigned k = 0; k < M; ++k) 
-            for (unsigned i = -Diags; i < N; ++i)
-                result.at(i + Diags, k) += d[i + Diags] * b.at(i, k);
-    } else {
-        for (unsigned k = 0; k < M; ++k) 
-            for (unsigned i = Diags; i < N; ++i)
-                result.at(i, k) += d[i - Diags] * b.at(i - Diags, k);
-    }
-    ) ...);
-
-    return result;
-}*/
-