#ifndef GELIB_d_matrix_H_
#define GELIB_d_matrix_H_

#include "../1_dimensional/d_vector.h"
#include <algorithm>
#include <iterator>

namespace gelib_impl
{

	class d_matrix : public std::vector<double>
	{
		size_t m_rows;
		size_t m_cols;

	public:
		void elem_op(std::function< double(double, double) > op, const d_matrix& w, d_matrix& z) const;

		static void vector_op(const std::function<void(double&, double&)>& op, std::vector<double>::iterator mat_begin, d_vector& w, size_t total_rows, size_t total_cols, int row, int col);

		d_matrix operator+(const d_matrix& w) const;

		d_matrix operator-(const d_matrix& w) const;

		d_matrix operator*(const d_matrix& w) const;

		/*
		* replace elements at (row, col) with those of w
		* must have one coordinate == -1 to indicate dimension to insert
		*/
		d_matrix& insert_vector(d_vector& w, int row, int col);

		/*
		* replace elements at (row, col) with those of w
		* must have one coordinate == -1 to indicate dimension to insert
		*/
		static void insert_vector(std::vector<double>::iterator mat_begin, d_vector& w, size_t total_rows, size_t total_cols, int row, int col);

		/*
		* copy of vector at (row, col)
		* must have one coordinate == -1 to indicate dimension to extract
		*/
		d_vector extract_vector(int row, int col) const;

		/*
		* copy of vector at (row, col)
		* must have one coordinate == -1 to indicate dimension to extract
		*/
		void extract_vector(d_vector& w, int row, int col) const;

		static d_vector extract_vector(std::vector<double>::const_iterator const& mat_begin, size_t total_rows, size_t total_cols, int row, int col);

		static void extract_vector(d_vector& w, std::vector<double>::const_iterator const& mat_begin, size_t total_rows, size_t total_cols, int row, int col);

		const size_t rows() const;

		const size_t cols() const;

		/*
		* total number of elements in the matrix
		*/
		const size_t no_elements() const;

		/*
		* convert to vector of vectors matrix
		*/
		std::vector<std::vector<double>> to_vec_of_vec();

		/*
		* allow construction for iterator like vector
		*/
		d_matrix(d_matrix::const_iterator it1, d_matrix::const_iterator it2, size_t rows, size_t cols);

		d_matrix(std::initializer_list<double> vals, size_t rows, size_t cols);
		/*
		* create from vector of vectors
		*/
		d_matrix(std::vector<std::vector<double>>& mat);

		d_matrix(double* p1, double* p2, size_t rows, size_t cols);

		d_matrix(size_t rows, size_t cols);

		d_matrix();

		d_matrix& fix_dimensions(size_t rows, size_t cols);

		size_t index(size_t row, size_t column) const;

		static size_t index(size_t total_cols, size_t row, size_t column);

		/* load this vector into the matrix as if double-layered std::vector*/
		void push_back(d_vector& v);

		double at(size_t row, size_t column) const;

		double operator()(size_t row, size_t column) const;

		void resize(size_t rows, size_t columns);

		void transpose();

		// column vector times matrix
		void mat_multiply_cr(d_vector& v, d_vector& out) const;
		// row vector times matrix
		void mat_multiply_rc(d_vector& v, d_vector& out) const;

	};

}

#endif /* D_MATRIX_H_ */
