#pragma once

#include <vector>
#include <exception>

#include <cassert>



namespace lm
{


class invalid_product_exception : std::exception
{};


//! 可変サイズ行列.
/*
 size_x == 3 , size_y == 2 とすると
 [ (0,0) , (1,0) , (2,0) ]
 [ (0,1) , (1,1) , (2,1) ]

 通しインデックスは
 [ 0 , 1 , 2 ]
 [ 3 , 4 , 5 ]
*/ 
template<typename T>
class matrix
{
public:
	matrix(void);
	matrix(size_t num_x, size_t num_y);
	matrix(size_t num_x, size_t num_y, const T& initial_value);
	matrix(const matrix<T>& m);

	T* v(void);
	const T* v(void) const;

	T& at(size_t idx_v);
	const T& at(size_t idx_v) const;
	T& operator()(size_t idx_x, size_t idx_y);
	const T& operator()(size_t idx_x, size_t idx_y) const;
	T& operator[](size_t idx_v);
	const T& operator[](size_t idx_v) const;

	T& at(size_t idx_x, size_t idx_y);
	const T& at(size_t idx_x, size_t idx_y) const;
	T& operator()(size_t idx_v);
	const T& operator()(size_t idx_v) const;

	void resize(size_t num_x, size_t num_y);
	void resize(size_t num_x, size_t num_y, const T& initial_value);

	void clear(void);

	size_t size(void) const;
	size_t size_x(void) const;
	size_t size_y(void) const;

	bool empty(void) const;

	void set(const matrix<T>& m);
	matrix<T>& operator=(const matrix<T>& m);

	bool equals(const matrix<T>& m);
	bool operator==(const matrix<T>& m);
	bool operator!=(const matrix<T>& m);

	typedef typename std::vector<T>::iterator               iterator;
	typedef typename std::vector<T>::const_iterator         const_iterator;
	typedef typename std::vector<T>::reverse_iterator       reverse_iterator;
	typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
	iterator               begin  (void);
	iterator               end    (void);
	reverse_iterator       rbegin (void);
	reverse_iterator       rend   (void);
	const_iterator         begin  (void) const;
	const_iterator         end    (void) const;
	const_reverse_iterator rbegin (void) const;
	const_reverse_iterator rend   (void) const;

	void swap(matrix<T>& m);

protected:
	size_t m_size_x;
	size_t m_size_y;
	std::vector<T> buf;
};


// global method 

template<typename T> inline
std::vector<T> dot(const matrix<T>& _m, const std::vector<T>& _v);
template<typename T> inline
std::vector<T> dot(const std::vector<T>& _v, const matrix<T>& _m);

template<typename T> inline
std::vector<T> operator*(const matrix<T>& _m, const std::vector<T>& _v);
template<typename T> inline
std::vector<T> operator*(const std::vector<T>& _v, const matrix<T>& _m);

template<typename T> inline
matrix<T> dot(const matrix<T>& _m0, const matrix<T>& _m1);
template<typename T> inline
matrix<T> operator*(const matrix<T>& _m0, const matrix<T>& _m1);


// implement

template<typename T> inline
matrix<T>::matrix(void)
	: m_size_x( 0 )
	, m_size_y( 0 )
{}

template<typename T> inline
matrix<T>::matrix(size_t num_x, size_t num_y)
	: m_size_x( num_x )
	, m_size_y( num_y )
{
	buf.resize( m_size_x * m_size_y );
}

template<typename T> inline
matrix<T>::matrix(size_t num_x, size_t num_y, const T& initial_value)
	: m_size_x( num_x )
	, m_size_y( num_y )
{
	buf.resize( m_size_x * m_size_y );
	std::fill( buf.begin() , buf.end() , initial_value );
}

template<typename T> inline
matrix<T>::matrix(const matrix<T>& m)
{
	buf = m.buf;
	m_size_x = m.size_x();
	m_size_y = m.size_y();
}


template<typename T> inline
T* matrix<T>::v()
{
	if( buf.empty() ) return NULL    ;
	else              return &buf[0] ;
}
template<typename T> inline
const T* matrix<T>::v() const
{
	if( buf.empty() ) return NULL    ;
	else              return &buf[0] ;
}


template<typename T> inline
T& matrix<T>::at(size_t idx_v)
{
	return buf[idx_v];
}
template<typename T> inline
const T& matrix<T>::at(size_t idx_v) const
{
	return buf[idx_v];
}
template<typename T> inline
T& matrix<T>::operator()(size_t idx_v)
{
	return at(idx_v);
}
template<typename T> inline
const T& matrix<T>::operator()(size_t idx_v) const
{
	return at(idx_v);
}
template<typename T> inline
T& matrix<T>::operator[](size_t idx_v)
{
	return at( idx_v );
}
template<typename T> inline
const T& matrix<T>::operator[](size_t idx_v) const
{
	return at( idx_v );
}


template<typename T> inline
T& matrix<T>::at(size_t idx_x, size_t idx_y)
{
	// Matrix subscript out of range
	assert( idx_x < m_size_x && idx_y < m_size_y );

	return buf[ idx_x + m_size_x * idx_y ];
}
template<typename T> inline
const T& matrix<T>::at(size_t idx_x, size_t idx_y) const
{
	// Matrix subscript out of range
	assert( idx_x < m_size_x && idx_y < m_size_y );

	return buf[ idx_x + m_size_x * idx_y ];
}
template<typename T> inline
T& matrix<T>::operator()(size_t idx_x,size_t idx_y)
{
	return at( idx_x , idx_y );
}
template<typename T> inline
const T& matrix<T>::operator()(size_t idx_x, size_t idx_y) const
{
	return at( idx_x , idx_y );
}


template<typename T> inline
void matrix<T>::resize(size_t num_x, size_t num_y)
{
	m_size_x = num_x ;
	m_size_y = num_y ;
	buf.resize( m_size_x * m_size_y );
}
template<typename T> inline
void matrix<T>::resize(size_t num_x, size_t num_y, const T& initial_value)
{
	m_size_x = num_x ;
	m_size_y = num_y ;
	buf.resize( m_size_x * m_size_y );
	std::fill( buf.begin() , buf.end() , initial_value );
}


template<typename T> inline
void matrix<T>::clear(void)
{
	m_size_x = 0;
	m_size_y = 0;
	buf.clear();
}


template<typename T> inline
size_t matrix<T>::size(void) const
{
	return buf.size();
}
template<typename T> inline
size_t matrix<T>::size_x(void) const
{
	return m_size_x;
}
template<typename T> inline
size_t matrix<T>::size_y(void) const
{
	return m_size_y;
}

template<typename T> inline
bool matrix<T>::empty(void) const
{
	return buf.empty();
}


template<typename T> inline
void matrix<T>::set(const matrix<T>& m)
{
	buf = m.buf;
	m_size_x = m.size_x() ;   m_size_y = m.size_y() ;
}
template<typename T> inline
matrix<T>& matrix<T>::operator=(const matrix<T>& m)
{
	set( m ) ;
	return (*this);
}


template<typename T> inline
bool matrix<T>::equals(const matrix<T>& m)
{
	if( m_size_x != m.m_size_x || m_size_y != m.m_size_y ) return false;
	for(size_t i=0;i<m.size();i++) if( (*this)(i) != m(i) ) return false;
	return true;
}
template<typename T> inline
bool matrix<T>::operator==(const matrix<T>& m)
{
	return equals( m );
}
template<typename T> inline
bool matrix<T>::operator!=(const matrix<T>& m)
{
	return !equals( m );
}


template<typename T> inline
typename matrix<T>::iterator matrix<T>::begin(void)
{
	return buf.begin();
}
template<typename T> inline
typename matrix<T>::iterator matrix<T>::end(void)
{
	return buf.end();
}
template<typename T> inline
typename matrix<T>::reverse_iterator matrix<T>::rbegin(void)
{
	return buf.rbegin();
}
template<typename T> inline
typename matrix<T>::reverse_iterator matrix<T>::rend(void)
{
	return buf.rend();
}
template<typename T> inline
typename matrix<T>::const_iterator matrix<T>::begin(void) const
{
	return buf.begin();
}
template<typename T> inline
typename matrix<T>::const_iterator matrix<T>::end(void) const
{
	return buf.end();
}
template<typename T> inline
typename matrix<T>::const_reverse_iterator matrix<T>::rbegin(void) const
{
	return buf.rbegin();
}
template<typename T> inline
typename matrix<T>::const_reverse_iterator matrix<T>::rend(void) const
{
	return buf.rend();
}

template<typename T> inline
void matrix<T>::swap(matrix<T>& m)
{
	std::swap(m_size_x, m.m_size_x);
	std::swap(m_size_y, m.m_size_y);
	buf.swap(m.buf);
}


template<typename T> inline
std::vector<T> dot(const matrix<T>& _m, const std::vector<T>& _v)
{
	return _m * _v;
}

template<typename T> inline
std::vector<T> dot(const std::vector<T>& _v, const matrix<T>& _m)
{
	return _v * _m;
}

//! 行列とベクトルの積
//! m.size_x() == 3 , m.size_y() == 5 の場合
//! [r0]   [m00 m10 m20]
//! [r1]   [m01 m11 m21]   [v0]
//! [r2] = [m02 m12 m22] * [v1]
//! [r3]   [m03 m13 m23]   [v2]
//! [r4]   [m04 m14 m24]
//!
//! ※ mij = m( i , j ) = m.at( i , j )
template<typename T> inline
std::vector<T> operator*(const matrix<T>& _m, const std::vector<T>& _v)
{
	assert( _m.size_x() == _v.size() );
	if( _m.size_x() != _v.size() )
		throw invalid_product_exception();

	std::vector<T> ret( _m.size_y() );
	for( size_t i = 0 ; i < _m.size_y() ; ++i )
	{
		T& val = ret[i];

		for( size_t j = 0 ; j < _m.size_x() ; ++j )
		{
			val += _m( j , i ) * _v[j];
		}
	}

	return ret;
}

//! 行列とベクトルの積
//! m.size_x() == 5 , m.size_y() == 3 の場合
//!                                 [m00 m10 m20 m30 m40]
//! [r0 r1 r2 r3 r4] = [v0 v1 v2] * [m01 m11 m21 m31 m41]
//!                                 [m02 m12 m22 m32 m42]
//!
//! ※ mij = m( i , j ) = m.at( i , j )
template<typename T> inline
std::vector<T> operator*(const std::vector<T>& _v, const matrix<T>& _m)
{
	assert( _m.size_y() == _v.size() );
	if( _m.size_y() != _v.size() )
		throw invalid_product_exception();

	std::vector<T> ret( _m.size_x() );
	for( size_t i = 0 ; i < _m.size_x() ; ++i )
	{
		T& val = ret[i];

		for( size_t j = 0 ; j < _m.size_y() ; ++j )
		{
			val += _v[j] * _m( i , j );
		}
	}

	return ret;
}

template<typename T> inline
matrix<T> dot(const matrix<T>& _m0, const matrix<T>& _m1)
{
	return _m0 * _m1;
}

//! 行列同士の積
//! A = M * N =>
//! [a00 a10]   [m00 m10 m 20]
//! [a01 a11]   [m01 m11 m 21]   [n00 n10]
//! [a02 a12] = [m02 m12 m 22] * [n01 n11]
//! [a03 a13]   [m03 m13 m 23]   [n02 n12]
//! [a04 a14]   [m04 m14 m 24]
template<typename T> inline
matrix<T> operator*(const matrix<T>& _m0, const matrix<T>& _m1)
{
	assert( _m0.size_x() == _m1.size_y() );
	if( _m0.size_x() != _m1.size_y() )
		throw invalid_product_exception();

	matrix<T> m( _m1.size_x() , _m0.size_y() );
	for( size_t i = 0 ; i < _m1.size_x() ; ++i )
	{
		for( size_t j = 0 ; j < _m0.size_y() ; ++j )
		{
			T& val = m( i , j );

			for( size_t k = 1 ; k < _m0.size_x() ; ++k )
			{
				val += _m0( k , j ) * _m1( i , k );
			}
		}
	}

	return m;
}


}
