CRNMatrix.h

Go to the documentation of this file.

/* Copyright 2014-2016 INSA-Lyon, Université Paris Descartes, ENS-Lyon
 * 
 * This file is part of libcrn.
 * 
 * libcrn is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * libcrn 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 Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with libcrn.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * file: CRNMatrix.h
 * \author Yann LEYDIER
 */

#ifndef CRNMATRIX_HEADER
#define CRNMATRIX_HEADER

#include <CRNObject.h>
#include <CRNException.h>
#include <CRNMath/CRNMath.h>
#include <CRNString.h>
#include <vector>

namespace crn
{
    /****************************************************************************/
    template<typename T> class Matrix: public Object
    {
        public:
            Matrix(size_t nrow, size_t ncol, const T &value = T(0)) : rows(nrow), cols(ncol)
            {
                if ((rows == 0) || (cols == 0))
                    throw ExceptionDomain("Matrix::Matrix(size_t nrow, size_t ncol, double value): null or negative row or column dimensions");
                datatype(rows * cols, value).swap(data);
            }

            Matrix(const std::vector<std::vector<T>> &m)
            {
                if (m.empty())
                    throw ExceptionDimension("Matrix::Matrix(m): empty matrix.");
                if (m.front().empty())
                    throw ExceptionDimension("Matrix::Matrix(m): empty matrix.");
                auto r = m.size();
                auto c = m.front().size();
                for (const auto &line : m)
                    if (line.size() != c)
                        throw ExceptionInvalidArgument("Matrix::Matrix(m): the argument is not a matrix.");
                rows = r;
                cols = c;
                datatype(rows * cols).swap(data);
                for (size_t r = 0; r < rows; ++r)
                    std::copy(m[r].begin(), m[r].end(), data.begin() + r * cols);
            }

            Matrix(const std::vector<T> &vect, Orientation ori = Orientation::VERTICAL)
            {
                if (vect.empty())
                    throw ExceptionDimension("Matrix::Matrix(vect, orientation): empty matrix.");
                if (ori == Orientation::VERTICAL)
                {
                    rows = vect.size();
                    cols = 1;
                    data = vect;
                }
                else if (ori == Orientation::HORIZONTAL)
                {
                    rows = 1;
                    cols = vect.size();
                    data = vect;
                }
                else
                    throw ExceptionInvalidArgument("Matrix::Matrix(vect, orientation): invalid orientation.");
            }

            Matrix(std::vector<T> &&vect, Orientation ori = Orientation::VERTICAL)
            {
                if (vect.empty())
                    throw ExceptionDimension("Matrix::Matrix(vect, orientation): empty matrix.");
                if (ori == Orientation::VERTICAL)
                {
                    rows = vect.size();
                    cols = 1;
                    data = std::move(vect);
                }
                else if (ori == Orientation::HORIZONTAL)
                {
                    rows = 1;
                    cols = vect.size();
                    data = std::move(vect);
                }
                else
                    throw ExceptionInvalidArgument("Matrix::Matrix(vect, orientation): invalid orientation.");
            }

            virtual ~Matrix() override = default;
            Matrix(const Matrix &) = default;
            template<typename Y> Matrix(const Matrix<Y> &m) : rows(m.GetRows()), cols(m.GetCols())
            {
                datatype(rows * cols).swap(data);
                std::copy(m.Std().begin(), m.Std().end(), data.begin());
            }
            Matrix(Matrix &&) = default;
            Matrix& operator=(const Matrix &) = default;
            Matrix& operator=(Matrix &&) = default;

            inline void Swap(Matrix &m)
            {
                std::swap(rows, m.rows);
                std::swap(cols, m.cols);
                std::swap(data, m.data);
            }

            inline size_t GetRows() const noexcept { return rows; }

            inline size_t GetCols() const noexcept { return cols; }

            inline const T& At(size_t pos) const noexcept { return data[pos]; }
            inline T& At(size_t pos) noexcept { return data[pos]; }
            inline const T& At(size_t r, size_t c) const noexcept { return data[r * cols + c]; }
            inline T& At(size_t r, size_t c) noexcept { return data[r * cols + c]; }

            void SetAll(const T &v)
            {
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        At(r, c) = v;
            }

            void Negative()
            {
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        At(r, c) = -At(r, c);
            }

            inline void SetRow(size_t r, const std::vector<T> &row)
            {
                if (r > rows)
                    throw ExceptionDomain("Matrix::SetRow(): index out of range.");
                if (row.size() != cols)
                    throw ExceptionDimension("Matrix::SetRow(): argument is not of the right size.");
                for (size_t c = 0; c < cols; ++c)
                    At(r, c) = row[c];
            }

            inline void SetColum(size_t c, const std::vector<T> &col)
            {
                if (c > cols)
                    throw ExceptionDomain("Matrix::SetColumn(): index out of range.");
                if (col.size() != rows)
                    throw ExceptionDimension("Matrix::SetColumn(): argument is not of the right size.");
                for (size_t r = 0; r < rows; ++r)
                    At(r, c) = col[r];
            }

            inline void IncreaseElement(size_t r, size_t c, const T &delta)
            {
                if (!areValidIndexes(r, c))
                    throw ExceptionDomain("Matrix::IncreaseElement(): index out of range.");
                At(r, c) += delta;
            }

            void IncreaseRow(size_t r, const T &v)
            {
                if (r >= rows)
                    throw ExceptionDomain("Matrix::IncreaseRow(): column index out of range");
                for (size_t c = 0; c < cols; ++c)
                    At(r, c) += v;
            }

            void IncreaseColumn(size_t c, const T &v)
            {
                if (c >= cols)
                    throw ExceptionDomain("Matrix::IncreaseColumn(): column index out of range");
                for (size_t r = 0; r < rows; ++r)
                    At(r, c) += v;
            }

            Matrix& operator+=(const Matrix &m)
            {
                if ((m.GetRows() != rows) || (m.GetCols() != cols))
                    throw ExceptionDimension("Matrix::+=(): incompatible dimensions");
                for (auto i : Range(data))
                    At(i) += m.At(i);
                return *this;
            }
            Matrix& operator-=(const Matrix &m)
            {
                if ((m.GetRows() != rows) || (m.GetCols() != cols))
                    throw ExceptionDimension("Matrix::+=(): incompatible dimensions");
                for (auto i : Range(data))
                    At(i) -= m.At(i);
                return *this;
            }
            Matrix& operator*=(const Matrix &m)
            { 
                if (cols != m.GetRows())
                    throw ExceptionDimension("Matrix::*=(): incompatible dimensions");
                auto p = datatype(rows * m.GetCols());
                for (size_t r = 0; r < rows; r++)
                    for (size_t c = 0; c < m.GetCols(); c++)
                        for (size_t k = 0; k < cols; k++)
                            p[r * m.GetCols() + c] += At(r, k) * m[k][c];                           
                cols = m.GetCols();
                data.swap(p);
                return *this;
            }
            bool operator==(const Matrix &m) const
            {
                if ((m.GetRows() != rows) || (m.GetCols() != cols))
                    return false;
                return std::equal(data.begin(), data.end(), m.data.begin());
            }
            inline bool operator!=(const Matrix &m) const
            { return !(*this == m); }

            Matrix& operator*=(double d) {  for (auto &v : data) v = T(v * d); return *this; }

            Matrix& operator/=(double d) { for (auto &v : data) v = T(v / d); return *this; }

            void MultRow(size_t r, double v)
            {
                if (r >= rows)
                    throw ExceptionDomain("Matrix::ScaleRow(): row index out of range");
                for (size_t c = 0; c < cols; ++c)
                    At(r, c) = T(v * At(r, c));
            }

            void MultColumn(size_t c, double v)
            {
                if (c >= cols)
                    throw ExceptionDomain("Matrix::ScaleColumn(): column index out of range");
                for (size_t r = 0; r < rows; ++r)
                    At(r, c) = T(v * At(r, c));
            }

            void CenterColumns()
            {
                auto means = std::vector<T>(cols, T(0));
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        means[c] += At(r, c);
                for (size_t c = 0; c < cols; ++c)
                    means[c] = means[c] / T(rows);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        At(r, c) -= means[c];
            }

            void ReduceColumns(bool centered = true)
            {
                auto means = std::vector<T>(cols, T(0));
                if (!centered)
                {
                    for (size_t r = 0; r < rows; ++r)
                        for (size_t c = 0; c < cols; ++c)
                            means[c] += At(r, c);
                    for (size_t c = 0; c < cols; ++c)
                        means[c] = means[c] / T(rows);
                }
                auto devs = std::vector<T>(cols, T(0));
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        devs[c] += Sqr(At(r, c) - means[c]);
                for (size_t c = 0; c < cols; ++c)
                    devs[c] = sqrt(devs[c] / T(rows));
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        At(r, c) /= devs[c];
            }

            void SwapRows(size_t r1, size_t r2)
            {
                if ((r1 >= rows) || (r2 >= rows))
                    throw ExceptionDomain("Matrix::SwapRows(): index out of range");
                if (r1 != r2)
                    for (size_t c = 0; c < cols; ++c)
                        std::swap(At(r1, c), At(r2, c));
            }

            void SwapColumns(size_t c1, size_t c2)
            {
                if ((c1 >= cols) || (c2 >= cols))
                    throw ExceptionDomain("Matrix::SwapColumns(): index out of range");
                if (c1 != c2)
                    for (size_t r = 0; r < rows; ++r)
                        std::swap(At(r, c1), At(r, c2));
            }

            void NormalizeForConvolution()
            {
                auto sum = T(0);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        sum += Abs(At(r, c));
                if (sum != T(0))
                    for (size_t r = 0; r < rows; ++r)
                        for (size_t c = 0; c < cols; ++c)
                            At(r, c) /= sum;
            }

            typename TypeInfo<T>::SumType CumulateCells() const
            {
                auto sum = typename TypeInfo<T>::SumType(0);
                for (size_t r = 0; r < rows; r++)
                    for (size_t c = 0; c < cols; c++)
                        sum += At(r, c);
                return sum;
            }

            T GetMin() const
            {
                return *std::min_element(data.begin(), data.end());
            }

            T GetMax() const
            {
                return *std::max_element(data.begin(), data.end());
            }

            std::pair<size_t, size_t> Argmin() const
            {
                auto pos = std::min_element(data.begin(), data.end()) - data.begin();
                return std::make_pair(pos / cols, pos % cols);
            }

            std::pair<size_t, size_t> Argmax() const
            {
                auto pos = std::max_element(data.begin(), data.end()) - data.begin();
                return std::make_pair(pos / cols, pos % cols);
            }

            Matrix MakeTranspose()
            {
                auto out = Matrix(cols, rows);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        out[c][r] = At(r, c);
                return out;
            }

            virtual Matrix& Transpose()
            {
                datatype newdata(cols * rows);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        newdata[c * rows + r] = At(r, c);
                data.swap(newdata);
                std::swap(rows, cols);
                return *this;
            }

            Matrix MakeRowAsColumn(size_t r) const
            {
                if (r >= GetRows())
                    throw ExceptionDomain("Matrix::MakeRowAsColumn(): index out of range");
                return Matrix(datatype(data.begin() + r * cols, data.begin() + (r + 1) * cols), Orientation::VERTICAL);
            }

            Matrix MakeRow(size_t r) const
            {   
                if (r >= GetRows())
                    throw ExceptionDomain("Matrix::MakeRow(): row index out of range");
                return Matrix(datatype(data.begin() + r * cols, data.begin() + (r + 1) * cols), Orientation::VERTICAL);
            }

            Matrix MakeColumn(size_t k) const
            {
                if (k >= cols)
                    throw ExceptionDomain("Matrix::MakeColumn(): column index out of range");
                auto column = Matrix(rows, 1, 0.0);
                for (size_t r = 0; r < rows; r++)
                    column.At(r, 0) = At(r, k);
                return column;
            }

            size_t ArgmaxInRow(size_t r) const
            {
                if (r >= rows)
                    throw ExceptionDomain("Matrix::ArgmaxInRow(): index out of range");
                return std::max_element(data.begin() + r * cols, data.begin() + (r + 1) * cols) - (data.begin() + r * cols);
            }

            size_t ArgmaxInColumn(size_t c) const
            {
                if (c >= cols)
                    throw ExceptionDomain("Matrix::ArgmaxInColumn(): index out of range");
                auto max = At(0, c);
                auto index = size_t(0);
                for (size_t r = 0; r < rows; r++)
                {
                    const auto &mrc = At(r, c);
                    if (mrc > max)
                    {
                        max = mrc;
                        index = r;
                    }
                }
                return index;
            }

            size_t CountNullCellsInRow(size_t r) const
            {
                if (r >= rows)
                    throw ExceptionDomain("Matrix::CountNullCellsInRow(): row index out of range");
                return std::count(data.begin() + r * cols, data.begin() + (r + 1) * cols, T(0));
            }

            size_t CountNullCellsInColumn(size_t c) const
            {
                if (c >= rows)
                    throw ExceptionDomain("Matrix::CountNullCellsInColumn(): column index out of range");
                auto count = size_t(0);
                for (size_t r = 0; r < rows; ++r)
                    if (At(r, c) == T(0))
                        count += 1;
                return count;
            }

            size_t CountNullCells() const
            {
                return std::count(data.begin(), data.end(), T(0));
            }

            Matrix MakeColumnMeans() const
            {
                auto m = Matrix(1, cols);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        m.At(c) += At(r, c);
                m *= 1.0 / double(rows);
                return m;
            }

            Matrix MakeColumnDeviations(const Matrix &means) const
            {
                auto m = Matrix(1, cols);
                for (size_t r = 0; r < rows; ++r)
                    for (size_t c = 0; c < cols; ++c)
                        m.At(c) += Sqr(At(r, c) - means[c]);
                m *= 1.0 / double(rows);
                return m;
            }

            inline Matrix MakeColumnDeviations(bool zero_means = false) const
            {
                if (zero_means)
                    return MakeColumnDeviations(Matrix{1, cols});
                else
                    return MakeColumnDeviations(MakeColumnMeans());
            }

            Matrix MakeCovariance() const
            {
                auto cov = Matrix(cols, cols);
                for (size_t i = 0; i < cols; ++i)
                    for (size_t j = i; j < cols; ++j)
                    {
                        auto product = T(0);
                        for (size_t k = 0; k < rows; ++k)
                            product += At(k, i) * At(k, j);
                        // Each pattern in supposed equiprobable
                        // and thus weighted with 1/nbPatterns.
                        product /= T(rows);
                        cov[i][j] = product;
                        if (i != j)
                            cov[j][i] = product;
                    }
                return cov;
            }

            inline T* operator[](size_t r) { return data.data() + r * cols; }

            inline const T* operator[](size_t r) const { return data.data() + r * cols; }

            inline const std::vector<T>& Std() const & noexcept { return data; }
            inline std::vector<T> Std() && { return std::move(data); }

            template<typename IMG> IMG ToImage() const
            {
                auto img = IMG{cols, rows};
                std::copy(data.begin(), data.end(), img.Begin());
                return img;
            }

            String ToString() const
            {
                String s;
                for (size_t r = 0; r < rows; r++)
                {
                    for (size_t c = 0; c < cols; c++)
                    {
                        if (c) s += U" ";
                        s += At(r, c);
                    }
                    s += "\n";
                }
                return s;
            }

        protected:
            inline bool areValidIndexes(size_t r, size_t c) const { return (r < rows) && (c < cols); }

            using datatype = std::vector<T>; 
            datatype data; 
            size_t rows, cols; 
    };

    template<typename T> Matrix<T> operator+(const Matrix<T> &m1, const Matrix<T> &m2)
    { return Matrix<T>(m1) += m2; }
    template<typename T> Matrix<T> operator+(Matrix<T> &&m1, const Matrix<T> &m2)
    { return std::move(m1 += m2); }
    template<typename T> Matrix<T> operator+(const Matrix<T> &m1, Matrix<T> &&m2)
    { return std::move(m2 += m1); }

    template<typename T> Matrix<T> operator-(const Matrix<T> &m1, const Matrix<T> &m2)
    { return Matrix<T>(m1) -= m2; }
    template<typename T> Matrix<T> operator-(Matrix<T> &&m1, const Matrix<T> &m2)
    { return std::move(m1 -= m2); }

    template<typename T> Matrix<T> operator*(const Matrix<T> &m, double d)
    { return std::move(Matrix<T>(m) *= d); }
    template<typename T> Matrix<T> operator*(Matrix<T> &&m, double d)
    { return std::move(m *= d); }
    template<typename T> Matrix<T> operator*(double d, const Matrix<T> &m)
    { return Matrix<T>(m) *= d; }
    template<typename T> Matrix<T> operator*(double d, Matrix<T> &&m)
    { return std::move(m *= d); }
    
    template<typename T> Matrix<T> operator/(const Matrix<T> &m, double d)
    { return std::move(Matrix<T>(m) /= d); }
    template<typename T> Matrix<T> operator/(Matrix<T> &&m, double d)
    { return std::move(m /= d); }
    template<typename T> Matrix<T> operator/(double d, const Matrix<T> &m)
    { return Matrix<T>(m) /= d; }
    template<typename T> Matrix<T> operator/(double d, Matrix<T> &&m)
    { return std::move(m /= d); }   
    

    template<typename T> Matrix<T> operator*(const Matrix<T> &m1, const Matrix<T> &m2)
    { return Matrix<T>(m1) *= m2; }
    template<typename T> Matrix<T> operator*(Matrix<T> &&m1, const Matrix<T> &m2)
    { return std::move(m1 *= m2); }

    template<typename I> struct TypeInfo<Matrix<I>>
    {
        using SumType = Matrix<typename TypeInfo<I>::SumType>;
        using DiffType = Matrix<typename TypeInfo<I>::DiffType>;
        using DecimalType = Matrix<typename TypeInfo<I>::DecimalType>;
    };

}

namespace std
{
    template<typename T> void swap(crn::Matrix<T> &m1, crn::Matrix<T> &m2)
    {
        m1.Swap(m2);
    }
}
#endif