#ifndef MATRIX_H_
#define MATRIX_H_

#include <iostream>
#include <vector>
#include <exception>
#include <functional>

template <typename T, int m, int n>
class matrix
{
    T values[m*n];

    public:
        typedef matrix<T, 1, n> row_t;
        typedef matrix<T, m, 1> column_t;

        matrix()
        {
            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; j++)
                    values[i*n + j] = T();
        }

        matrix(std::initializer_list<std::initializer_list<T>> list)
        {
            int i = 0, j = 0;
            for (const auto& row : list) {
                j = 0;
                for (const auto& cell : row) {
                    set(i, j++, cell);
                }
                i++;
            }
        }

        T get(int i, int j) const
        {
            return values[i*n + j];
        }

        void set(int i, int j, T value)
        {
            values[i*n + j] = value;
        }

        template <int p> matrix<T, m, p> operator* (const matrix<T, n, p>& rhs) const
        {
            matrix<T, m, p> result;

            for (int i = 0; i < m; ++i) {
                for(int j = 0; j < p; j++) {
                    T accumulator = 0;
                    for(int k = 0; k < n; k++) 
                        accumulator += this->get(i, k) * rhs.get(k, j);

                    result.set(i, j, accumulator);
                }
            }

            return result;
        }

        matrix<T, m, n> operator* (const T& rhs) const
        {
            matrix <T, m, n> result;

            for (int i = 0; i < m; ++i)
                for (int j = 0; j < m; ++j) 
                    result.set(i, j, rhs * get(i, j));

            return result;
        }

        matrix<T, n, m> transpose() const
        {
            matrix<T, n, m> result;

            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j) 
                    result.set(j, i, get(i, j));


            return result;
        }

        matrix<T, m, n> operator- () const
        {
            matrix<T, m, n> result;

            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    result.set(i, j, -get(i, j));

            return result;
        }

        matrix<T, m, n> operator- (const matrix<T, m, n>& rhs) const
        {
            return *this + (-rhs);
        }

        matrix<T, m, n> operator- (const T& value) const
        {
            return *this + (-value);
        }

        matrix<T, m, n> operator+ (const matrix<T, m, n>& rhs) const
        {
            matrix<T, m, n> result;

            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    result.set(i, j, this->get(i, j) + rhs.get(i, j));

            return result;
        }

        matrix<T, m, n> operator+ (const T& value) const
        {
            matrix <T, m, n> result;

            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    result.set(this->get(i, j) + value);

            return result;
        }

        column_t column(std::size_t j)
        {
            column_t result;

            for (int i = 0; i < m; ++i) {
                result.set(i, 0, get(i, j));
            }
            
            return result;
        }

        row_t row(std::size_t i)
        {
            row_t result;

            for (int j = 0; j < n; ++j) {
                result.set(0, j, get(i, j));
            }
            
            return result;
        }

        friend std::ostream& operator<< (std::ostream& stream, const matrix<T, m, n>& matrix)
        {
            for (int i = 0; i < m; ++i) {
                for (int j = 0; j < n; j++) {
                    stream << matrix.get(i, j) << " ";
                }
                stream << std::endl;
            }

            return stream;
        }

        std::ostream& save(std::ostream& stream)
        {
            char start[2] = { char(0xDE), char(0xAD) };
            char end[2]   = { char(0xBE), char(0xEF) };

            T temp;
            stream.write(start, 2);
            for (int i = 0; i < m; ++i) {
                for (int j = 0; j < n; j++) {
                    temp = get(i, j);
                    stream.write((char*)(&temp), sizeof(T));
                }
            }
            stream.write(end, 2);

            return stream;
        }

        void load(std::istream& stream)
        {
            char buf[2];

            stream.read(buf, 2);
            if (buf[0] != char(0xDE) || buf[1] != char(0xAD)) {
                std::cerr << "Bledny prefix" << std::endl;
                return;
            }

            T temp;
            for (int i = 0; i < m; ++i) {
                for (int j = 0; j < n; j++) {
                    stream.read((char*)(&temp), sizeof(T));
                    set(i, j, temp);
                }
            }

            stream.read(buf, 2);
            if (buf[0] != char(0xBE) || buf[1] != char(0xEF))  {
                std::cerr << "Bledny suffix" << std::endl;
                return;
            }
        }
};

template <typename T, int m, int n> matrix<T, m, n> operator*(const T& lhs, const matrix<T, m, n>& rhs) {
    return rhs * lhs;
}

template <typename T, int n> using vector = matrix<T, n, 1>;

template <typename T, int m, int n> void fill(matrix<T, m, n> &mat, std::function<T(const int&, const int&)> filler)
{
    for (int i = 0; i < m; i++)
        for (int j = 0; j < n; j++)
            mat.set(i, j, filler(i, j));
}

template <typename T, int n, int m> 
vector<T, m+n> concat(const vector<T, n> &a, const vector<T, m> &b)  
{
    vector<T, m+n> result;

    for (int i = 0; i < n; i++)
        result.set(i, 0, a.get(i, 0));
    for (int i = 0; i < m; i++)
        result.set(i+n, 0, b.get(i, 0)); // YOU LIITLE FUCKING BASTARD

    return result;
}

template <typename T, int n>
vector<T, n> normalize(const vector<T, n> &vec)
{
    T accumulator = T();

    for (int i = 0; i < n; i++) {
        T temp = vec.get(i, 0);
        accumulator += abs(temp);
    }

    accumulator /= n;

    if (!accumulator) return vec;

    std::function<T(const T&)> normalizer = [accumulator](const T& item) { return item / accumulator; };
    return map(vec, normalizer);
}

template <typename T, int m, int n> 
matrix<T, m, n> map(const matrix<T, m, n>& matrix, std::function<T(const T&)> func)
{
    ::matrix<T, m, n> result;

    for (int i = 0; i < m; i++)
        for (int j = 0; j < n; j++)
            result.set(i, j, func(matrix.get(i, j)));

    return result;
}

template <typename T, int m, int n> 
matrix<T, m, n> combine(const matrix<T, m, n>& a, const matrix<T, m, n>& b, std::function<T(const T&, const T&)> func)
{
    ::matrix<T, m, n> result;

    for (int i = 0; i < m; i++)
        for (int j = 0; j < n; j++)
            result.set(i, j, func(a.get(i, j), b.get(i, j)));

    return result;
}

template <typename T, int m, int n> 
matrix<T, m, n> combine(const matrix<T, m, n>& a, const matrix<T, m, n>& b, T c, T d)
{
    return a*c + d*b;
}

#endif