#ifndef MATRIX_HPP_
#define MATRIX_HPP_

#include <utility>
#include <vector>
#include <sstream>
#include <stdexcept>

/**
 * @class matrix
 */
template < typename T, typename Allocator = std::allocator<T> >
class Matrix
{
	// Enums
public:
	enum StorageLayout
	{
		RowMajorLayout = 0,
		ColumnMajorLayout
	};

	// Types
public:
	typedef T ValueType;
	typedef Allocator AllocatorType;

public:

	/**
	 * createLayoutStrides
	 * @brief Create strides for a given layout.
	 * @param rows Number of rows in the matrix.
	 * @param cols Number of columns in the matrix.
	 * @param layout Layout of the matrix (row-major or column-major).
	 * @return Stride pair.
	 */
	static std::pair<int, int> createLayoutStrides(
		const int& rows, const int& cols,
		const int& layout)
	{
		// default is row-major (i.e. for layout != ColumnMajorLayout)
		std::pair<int, int> stride(cols, 1);
		if(layout == ColumnMajorLayout)
			stride = std::make_pair(1, rows);
		return stride;
	}

	/**
	 * Constructor
	 */
	Matrix(const AllocatorType& alloc = AllocatorType())
		: rows_(0), cols_(0), layout_(RowMajorLayout),
		  strides_(createLayoutStrides(rows_, cols_, layout_)),
		  data_(alloc), ptr_(0), ownData_(true)
	{};

	Matrix(const int& rows, const int& cols,
	       const ValueType& val, const int& layout = RowMajorLayout, 
	       const AllocatorType alloc = AllocatorType())
		: rows_(rows), cols_(cols), layout_(layout),
		  strides_(createLayoutStrides(rows_, cols_, layout_)),
		  data_(rows_*cols_, val, alloc), ptr_(&data_[0]), ownData_(true)
	{};

	Matrix(const int& rows, const int& cols,
	       const AllocatorType alloc = AllocatorType())
		: rows_(rows), cols_(cols), layout_(RowMajorLayout),
		  strides_(createLayoutStrides(rows_, cols_, layout_)),
		  data_(rows_*cols_, ValueType(), alloc), ptr_(&data_[0]), ownData_(true)
	{};

	Matrix(ValueType* ptr, const int& rows, const int& cols,
	       const int& layout = RowMajorLayout, 
	       const AllocatorType alloc = AllocatorType())
		: rows_(rows), cols_(cols), layout_(layout),
		  strides_(createLayoutStrides(rows_, cols_, layout_)),
		  data_(), ptr_(ptr), ownData_(false)
	{}

	Matrix(const Matrix<ValueType, AllocatorType>& other)
	: rows_(other.rows()), cols_(other.columns()), layout_(other.layout_),
	  strides_(createLayoutStrides(rows_, cols_, layout_)),
	  data_(other.data_), ptr_(&data_[0]), ownData_(true)
	{};

	/**
	 * Destructor
	 */
	virtual ~Matrix(){};

	/**
	 * setLayout
	 * @brief Set the storage layout of the matrix (row-major or column-major)
	 * @param layout Layout to be set.
	 */
	void setLayout(const int& layout)
	{
		layout_ = layout;
		strides_ = createLayoutStrides(rows_, cols_, layout_);
	}

	/**
	 * resize
	 * @brief Resize matrix, if the underlying data is owned by the matrix.
	 * @param rows New number of rows.
	 * @param cols New number of columns.
	 * @return 0 if successful, -1 otherwise.
	 */
	int resize(const int& rows, const int& cols)
	{
		// If matrix does not own the underlying data,
		// we cannot resize the matrix.
		if(!ownData_)
			return -1;

		// check for valid arguments
		if(rows <= 0 || cols <= 0)
			return -1;

		// Resize matrix
		if(rows_ != rows || cols_ != cols)
		{
			rows_ = rows;
			cols_ = cols;
			strides_ = createLayoutStrides(rows_, cols_, layout_);
			data_.resize(rows*cols);
			ptr_ = &data_[0];
		}
		return 0;
	}

	/**
	 * fill
	 * @brief Fill the entire array with given value.
	 * @param val Value to be assigned to all elements in the matrix.
	 */
	void fill(const ValueType& val)
	{
		std::fill(ptr_, ptr_+rows_*cols_, val);
	}

	/**
	 * copy
	 * @brief Copy data in a range.
	 * @param first Iterator pointing to the start of data.
	 * @param last Iterator beyond the end of data.
	 */
	template <class InputIterator>
	void copy(InputIterator first, InputIterator last)
	{
		std::copy(first, last, ptr_);
	}

	/**
	 * equal
	 * @brief Compare matrix with another for its size and contents.
	 * @param other Matrix to be compared with.
	 * @return true if equal, false otherwise.
	 */
	bool equal(const Matrix& other) const
	{
		return rows_ == other.rows_ &&
			cols_ == other.cols_ &&
			std::equal(ptr_, ptr_+rows_*cols_, other.ptr_);
	}

	/**
	 * rows
	 * @brief Get number of rows in the matrix.
	 * @return Number of rows.
	 */
	int rows() const { return rows_; };

	/**
	 * columns
	 * @brief Get number of columns in the matrix.
	 * @return Number of columns.
	 */
	int columns() const { return cols_; };

	/**
	 * data
	 * @brief Get access to underlying data array.
	 * @return const pointer to matrix data array.
	 */
	const ValueType* data() const { return ptr_; };

	/**
	 * data
	 * @brief Get access to underlying data array.
	 * @return Pointer to matrix data array.
	 */
	ValueType* data() { return ptr_; };

	/**
	 * operator()
	 * @brief Access element in matrix by linear index (does not throw).
	 * @param linear Index of the element.
	 * @return const reference to element at linear index.
	 */
	const ValueType& operator()(const int& linear) const
	{ return ptr_[linear]; };

	/**
	 * operator()
	 * @brief Access element in matrix by linear index (does not throw).
	 * @param linear Index of the element.
	 * @return Reference to element at linear index.
	 */
	ValueType& operator()(const int& linear)
	{ return ptr_[linear]; };

	/**
	 * at
	 * @brief Access element in matrix by linear index. This function
	 * checks for range violations and would throw std::out_of_range
	 * in case of an invalid index.
	 * @param index Index of the element.
	 * @return Reference to element at linear index.
	 */
	const ValueType& at(const int& index) const throw(std::out_of_range)
	{
		if(index < 0 || index >= rows_*cols_)
		{
			throw std::out_of_range("[Matrix] linear index out of range");
		}
		return ptr_[index];
	};

	/**
	 * at
	 * @brief Access element in matrix by linear index. This function
	 * checks for range violations and would throw std::out_of_range
	 * in case of an invalid index.
	 * @param index Index of the element.
	 * @return Reference to element at linear index.
	 */
	ValueType& at(const int& index) throw(std::out_of_range)
	{
		if(index < 0 || index >= rows_*cols_)
		{
			throw std::out_of_range("[Matrix] linear index out of range");
		}
		return ptr_[index];
	};

	/**
	 * operator()
	 * @brief Access element in matrix by matrix index (does not throw).
	 * @param row Row of the element.
	 * @param col Column of the element.
	 * @return const reference to element at (row, col).
	 */
	const ValueType& operator()(const int& row, const int& col) const
	{ return operator()(linearIndex(row, col)); };

	/**
	 * operator()
	 * @brief Access element in matrix by matrix index (does not throw).
	 * @param row Row of the element.
	 * @param col Column of the element.
	 * @return Reference to element at (row, col).
	 */
	ValueType& operator()(const int& row, const int& col)
	{ return operator()(linearIndex(row, col)); };

	/**
	 * at
	 * @brief Access element in matrix by matrix index. This function
	 * checks for range violations and would throw std::out_of_range
	 * in case of an invalid index.
	 * @param row Row of the element.
	 * @param col Column of the element.
	 * @return const reference to element at (row, col).
	 */
	const ValueType& at(const int& row, const int& col) const throw(std::out_of_range)
	{
		if(row < 0 || row >= rows_ || col < 0 || col >= cols_)
		{
			throw std::out_of_range("[Matrix] matrix index out of range");
		}
		return operator()(linearIndex(row, col));
	};

	/**
	 * at
	 * @brief Access element in matrix by matrix index. This function
	 * checks for range violations and would throw std::out_of_range
	 * in case of an invalid index.
	 * @param row Row of the element.
	 * @param col Column of the element.
	 * @return Reference to element at (row, col).
	 */
	ValueType& at(const int& row, const int& col) throw(std::out_of_range)
	{
		if(row < 0 || row >= rows_ || col < 0 || col >= cols_)
		{
			throw std::out_of_range("[Matrix] matrix index out of range");
		}
		return operator()(linearIndex(row, col));
	};

private:

	/**
	 * linearIndex
	 * @brief Compute linear index from row-column.
	 * @param row Row of the element.
	 * @param col Column of the element.
	 * @return Linear index corresponding to (row, column).
	 */
	int linearIndex(const int& row, const int& col) const
	{
		return row*strides_.first + col*strides_.second;
	}

	// Implementation
protected:
	int rows_;
	int cols_;
	int layout_;
	std::pair<int, int> strides_;
	std::vector<ValueType, AllocatorType> data_;
	ValueType* ptr_;
	bool ownData_;
};

#endif /* MATRIX_HPP_ */