#ifndef __MATRIX_AMV
	#define __MATRIX_AMV

#pragma warning(disable:4786)

#include <iostream>
#include <strstream>
#include <cstdlib>
#include <new>
#include "valarray.h"
#include "bitarray.h"
#include "exception.h"
#include "function.h"
#include "deque.h"
#include "binaryTree.h"
#include "expression.h"
#include "amvdefs.h"

namespace amv	{

template <class type>
class matrix;

template <class type>
matrix<type> GaussJordan(const matrix<type> &m);

template <class type>
matrix<type> Gauss(const matrix<type> &m, size_t *pSwaps);

template <class type>
matrix<type> GaussX(const matrix<type> &m);

template <class type>
matrix<type> inverse(const matrix<type> &m);

template <class type>
matrix<type> identity(const matrix<type> &m);

template <class type>
matrix<type> transposed(const matrix<type> &m);

template <class type>
matrix<type> determinant(const matrix<type> &m);

template <class type>
matrix<type> plan(const matrix<type> &m);

template <class type>
matrix<type> mul_diag(const matrix<type> &m);

template <class type>
matrix<type> sin(const matrix<type> &m);

template <class type>
matrix<type> cos(const matrix<type> &m);

template <class type>
matrix<type> tan(const matrix<type> &m);

template <class type>
ostream &operator<<(ostream &os, const matrix<type> &m);

template <class type>
istream &operator>>(istream &is, matrix<type> &m);

template <class type>
class matrix	{
	size_t nRows,nColumns;
	valarray<type> *mat;
	ostream *pOstream;
	matrix<type> calculateInverse();
	type calculateDeterminant(bitarray &C, bitarray &R);
	bool arrangeRows(size_t beg=0);
	void deleteOtherItemsInColumn(size_t k);
	void deleteDownItemsInColumn(size_t k);
	public:
		class invalid_operator: public exception	{
			public:
				invalid_operator(const char *message="Error, operador no válido"):
				exception(message)
				{
				}

		};

		matrix(size_t numberRows=1, size_t numberColumns=1, ostream *pOs=0);
		matrix(const matrix<type> &m);
		matrix(const type &t, ostream *pOs=0);
		~matrix();
		void swapRows(size_t r1, size_t r2);
		valarray<type> &operator[](size_t subscript);
		const valarray<type> &operator[](size_t subscript) const;
		matrix<type> &operator=(const matrix<type> &m);
		matrix<type> &operator=(const type &value);
		size_t rows() const;
		size_t columns() const;
		matrix<type> operator+(const matrix<type> &m2);
		matrix<type> operator-(const matrix<type> &m2);
		matrix<type> operator*(const matrix<type> &m2);
		matrix<type> operator*(const type &factor);
		matrix<type> operator/(const matrix<type> &m2) { return *this; }
		matrix<type> operator^(const matrix<type> &m2) { return *this; }
		bool operator==(const matrix<type> &m2)
			{
			if(nRows==m2.nRows && nColumns==m2.nColumns)
				{
				for(size_t i=0; i<nRows; ++i)
					if(!((*this)[i]==m2[i]))
						return false;
				}
			else
				return false;

			return true;
			}
		
		static matrix<type> identity(size_t dimension);
		
		friend matrix<type> GaussJordan(const matrix<type> &m);
		friend matrix<type> Gauss(const matrix<type> &m, size_t *pSwaps);
		friend matrix<type> GaussX(const matrix<type> &m);
		friend matrix<type> inverse(const matrix<type> &m);
		friend matrix<type> identity(const matrix<type> &m);
		friend matrix<type> transposed(const matrix<type> &m);
		friend matrix<type> determinant(const matrix<type> &m);
		friend matrix<type> plan(const matrix<type> &m);
		friend matrix<type> mul_diag(const matrix<type> &m);
		friend matrix<type> sin(const matrix<type> &m);
		friend matrix<type> cos(const matrix<type> &m);
		friend matrix<type> tan(const matrix<type> &m);

		friend ostream &operator<<(ostream &os, const matrix<type> &m);
		friend istream &operator>>(istream &is, matrix<type> &m);

//agregados 15 oct 2003
		bool scalar; //permite saber si la matriz consta de un sólo elemento, equiv. a un escalar
	bool isScalar() const
		{
		return scalar;
		}
	void setScalar()
		{
		if(nRows==1 && nColumns==1)
			scalar=true;
		else
			scalar=false;
		}
	void setPOstream(ostream *pOs)
		{
		pOstream=pOs;
		}

	matrix<type> apply(type fn(type)) const
		{
		matrix temp(nRows,nColumns);
		for(size_t i=0; i<nRows; ++i)
			temp[i]=mat[i].apply(fn);

		return temp;
		}

	matrix<type> apply(type fn(const type &)) const
		{
		matrix temp(nRows,nColumns);
		for(size_t i=0; i<nRows; ++i)
			temp[i]=mat[i].apply(fn);

		return temp;
		}

	matrix<type> &apply_itself(type fn(type))
		{
		for(size_t i=0; i<nRows; ++i)
			mat[i]=mat[i].apply(fn);

		return *this;
		}

	matrix<type> &apply_itself(type fn(const type &))
		{
		for(size_t i=0; i<nRows; ++i)
			mat[i]=mat[i].apply(fn);

		return *this;
		}

	static binaryTree<function<matrix<type> > > getFunctionsTree(const matrix<type> &m);
//
};

	template <class type>
	binaryTree<function<matrix<type> > > matrix<type>::getFunctionsTree(const matrix<type> &m)
		{
		type xxx;
		matrix<type> m1;
		function<matrix<type> > f1;
		binaryTree<function<matrix<type> > > t;

		t.insert(function<matrix<type> > ("gauss",GaussX)),
		t.insert(function<matrix<type> > ("gj",GaussJordan)),
		t.insert(function<matrix<type> > ("inv",inverse)),
		t.insert(function<matrix<type> > ("trans",transposed)),
		t.insert(function<matrix<type> > ("det",determinant)),
		t.insert(function<matrix<type> > ("tr",plan)),
		t.insert(function<matrix<type> > ("muld",mul_diag)),
		t.insert(function<matrix<type> > ("sin",sin)),
		t.insert(function<matrix<type> > ("cos",cos)),
		t.insert(function<matrix<type> > ("tan",tan));

		return t;
		}

template <class type>
matrix<type> matrix<type>::calculateInverse()
{
size_t i,j;
matrix<type> temp(nRows,nColumns*2), inverseMatrix(nRows,nColumns);
for(i=0; i<nRows; i++)
	for(j=0; j<nColumns; j++)
		temp[i][j]=mat[i][j];
for(i=0; i<nRows; i++)
	temp[i][nColumns+i]=type(true);

for(j=0; j<nColumns; j++)
	{
	temp.arrangeRows(j);
	if(!(temp[j][j]==type()))
		{
		temp.deleteOtherItemsInColumn(j);
		}
	}

if(pOstream)
	*pOstream<<endl;
for(i=0; i<nRows; i++)
	for(j=0; j<nColumns; j++)
		inverseMatrix[i][j]=temp[i][j+nColumns];
return inverseMatrix;
}

template <class type>
type matrix<type>::calculateDeterminant(bitarray &C, bitarray &R)
{
type result=type();
size_t nRow=0, nItems=C.size();
if((R.count(true) == 1) && (C.count(true)==1))
	{
	result=mat[R.find(true)][C.find(true)];
	}
else
	{
	size_t nCols=0;
	bitarray RR(R);
	nRow = R.find(true);
	RR[nRow] =false;
	bitarray CC(C);
	for(size_t j = 0; j < nItems; j++)
		{
		if(CC[j])
			{
			CC[j]=false;
			nCols++;
			if(nCols%2==0)
				result+=-mat[nRow][j]*calculateDeterminant(CC,RR);
			else
				result+=mat[nRow][j]*calculateDeterminant(CC,RR);
			CC[j]=true;
			}
		}
	RR[nRow]=true;
	}
return result;
}

template <class type>
bool matrix<type>::arrangeRows(size_t beg)
{
size_t k;

size_t rowsLimit=nRows-1;
bool swap=false;


if(mat[beg][beg]==type())
	{
	for(k=(beg+1); k<nRows && !swap; ++k)
		{
		if(!(mat[k][beg]==type()))
			{
			swapRows(beg,k);
				if(pOstream)
				{
				*pOstream<<'r'<<beg<<"-><-r"<<k<<endl;
				*pOstream<<*this;
				}

			swap=true;
			}
		}
	}

return swap;
}

template <class type>
void matrix<type>::deleteOtherItemsInColumn(size_t k)
{
size_t i;
type divisor=mat[k][k], pivot;

mat[k]/=divisor;
if(pOstream)
	{
	*pOstream<<'r'<<k<<'/'<<divisor<<endl;
	*pOstream<<*this;
	}

for(i=0; i<nRows; i++)
	{
	if((i!=k) && (!(mat[i][k]==type())))
		{
		pivot=-mat[i][k];
		divisor=pivot;
		mat[k]*=pivot;

		if(pOstream)
			{
			*pOstream<<'r'<<k<<'*'<<pivot<<endl;
			*pOstream<<*this;
			mat[i]+=mat[k];
			*pOstream<<'r'<<i<<"+r"<<k<<endl;
			*pOstream<<*this;
			mat[k]/=divisor;
			*pOstream<<'r'<<k<<'/'<<divisor<<endl;
			*pOstream<<*this;
			}
		else

			{
			mat[i]+=mat[k];
			mat[k]/=divisor;
			}
		}
	}
}

template <class type>
void matrix<type>::deleteDownItemsInColumn(size_t k)
{
size_t i;
type reciprocal, pivot;

for(i=k+1; i<nRows; i++)
	{
	if(!(mat[i][k]==type()))
		{
		pivot=-mat[i][k];
		reciprocal=mat[k][k]/pivot;
		mat[k]*=pivot/mat[k][k];

		if(pOstream)
			{
			*pOstream<<'r'<<k<<'*'<<pivot<<endl;
			*pOstream<<*this;
			mat[i]+=mat[k];
			*pOstream<<'r'<<i<<"+r"<<k<<endl;
			*pOstream<<*this;
			mat[k]*=reciprocal;
			*pOstream<<'r'<<k<<'*'<<reciprocal<<endl;
			*pOstream<<*this;
			}
		else

			{
			mat[i]+=mat[k];
			mat[k]*=reciprocal;
			}
		}
	}
}


template <class type>
matrix<type>::matrix(size_t numberRows, size_t numberColumns, ostream *pOs):
nRows(numberRows), nColumns(numberColumns), pOstream(pOs)
{
size_t i;
mat=new valarray<type>[nRows];
for(i=0; i<nRows; i++)
	{
	mat[i].resize(nColumns);
	}
setScalar();
}

template <class type>
matrix<type>::matrix(const matrix<type> &m):
nRows(m.nRows), nColumns(m.nColumns), pOstream(m.pOstream), scalar(m.scalar)
{
size_t i;

mat=new valarray<type>[nRows];
for(i=0; i<nRows; i++)
	{
	mat[i].resize(nColumns);
	mat[i]=m[i];
	}
}

template <class type>
matrix<type>::matrix(const type &t, ostream *pOs):
nRows(1), nColumns(1), pOstream(pOs), scalar(true)
{
mat=new valarray<type>[1];
mat[0].resize(1);
mat[0][0]=t;

}

template <class type>
matrix<type>::~matrix()
{
delete [] mat;
}

template <class type>
void matrix<type>::swapRows(size_t r1, size_t r2)
{
valarray<type> temp(mat[r1]);

mat[r1]=mat[r2];
mat[r2]=temp;
}

template <class type>
valarray<type> &matrix<type>::operator[](size_t subscript)
{
return mat[subscript];
}

template <class type>
const valarray<type> &matrix<type>::operator[](size_t subscript) const
{
return mat[subscript];
}

template <class type>
matrix<type> &matrix<type>::operator=(const matrix<type> &m)
{
size_t i;

delete [] mat;

nRows=m.nRows;
nColumns=m.nColumns;
mat=new valarray<type>[nRows];
for(i=0; i<nRows; ++i)
	{
	mat[i].resize(nColumns);
	mat[i]=m[i];
	}
scalar=m.scalar;

return *this;
}

template <class type>
matrix<type> &matrix<type>::operator=(const type &value)
{
size_t i;
for(i=0; i<nRows; ++i)
	mat[i]=value;

return *this;
}

template <class type>
size_t matrix<type>::rows() const
{
return nRows;
}

template <class type>
size_t matrix<type>::columns() const
{
return nColumns;
}

template <class type>
matrix<type> matrix<type>::operator+(const matrix<type> &m2)
{
if(scalar && m2.scalar)
	return matrix<type>(mat[0][0]+m2.mat[0][0]);

if(nRows!=m2.nRows || nColumns!=m2.nColumns)
	throw invalid_operator();

size_t i;
matrix<type> temp(*this);
for(i=0; i<nRows; i++)
	temp[i]+=m2[i];
	
return temp;
}

template <class type>
matrix<type> matrix<type>::operator-(const matrix<type> &m2)
{
if(scalar && m2.scalar)
	return matrix<type>(mat[0][0]-m2.mat[0][0]);

if(nRows!=m2.nRows || nColumns!=m2.nColumns)
	throw invalid_operator();

size_t i;
matrix<type> temp(*this);
for(i=0; i<nRows; i++)
		temp[i]-=m2[i];
		
return temp;
}

template <class type>
matrix<type> matrix<type>::operator*(const matrix<type> &m2)
{
if(scalar)
	{
	matrix<type> temp(m2);
	temp=temp*mat[0][0];

	return temp;
	}
else if(m2.scalar)
	{
	matrix<type> temp(*this);
	temp=temp*m2.mat[0][0];

	return temp;
	}

if(nColumns!=m2.nRows)
	throw invalid_operator();

size_t i,j;
matrix<type> temp(nRows,m2.nColumns);
type sum;

for(i=0; i<nRows; i++)
	for(j=0; j<m2.nColumns; j++)
		{
		sum=0;
		for(size_t h=0; h<nColumns; h++)
			sum+=mat[i][h]*m2[h][j];
		temp[i][j]=sum;
		}
return temp;
}

template <class type>
matrix<type> matrix<type>::operator*(const type &factor)
{
size_t i;
matrix<type> temp(*this);
for(i=0; i<nRows; i++)
	temp[i]*=factor;
	
return temp;
}

template <class type>
matrix<type> matrix<type>::identity(size_t dimension)
{
matrix <type> temp(dimension,dimension);
size_t i;
for(i=0; i<dimension; ++i)
	temp[i][i]=type(true);

return temp;
}

template <class type>
matrix<type> GaussJordan(const matrix<type> &m)
{
size_t k;
size_t nVars=m.nColumns-1, limit;
matrix<type> temp(m);

if(temp.nRows>=temp.nColumns)
	limit=nVars;
else
	limit=temp.nRows;

for(k=0; k<limit; k++)
	{
	temp.arrangeRows(k);
	if(!(temp.mat[k][k]==type()))
		{
		temp.deleteOtherItemsInColumn(k);
		}
	}
if(m.pOstream)
	*m.pOstream<<endl;
return temp;
}

template <class type>
matrix<type> Gauss(const matrix<type> &m, size_t *pSwaps)
{
size_t k, swaps=0;
size_t nVars=m.nColumns-1, limit;
matrix<type> temp(m);

if(temp.nRows>=temp.nColumns)
	limit=nVars;
else
	limit=temp.nRows;

for(k=0; k<limit; k++)
	{
	if(temp.arrangeRows(k))
		++swaps;
	if(!(temp.mat[k][k]==type()))
		{
		temp.deleteDownItemsInColumn(k);
		}
	}

if(pSwaps)
	*pSwaps=swaps;

if(m.pOstream)
	*m.pOstream<<endl;
return temp;
}

template <class type>
matrix<type> GaussX(const matrix<type> &m)
{
return Gauss(m,0);
}

template <class type>
matrix<type> inverse(const matrix<type> &m)
{
matrix<type> det;

if(m.nRows!=m.nColumns)
	throw matrix<type>::invalid_operator();
det=determinant(m);
if(det[0][0]==type())
	throw matrix<type>::invalid_operator();

matrix<type> temp(m);

return temp.calculateInverse();
}

template <class type>
matrix<type> identity(const matrix<type> &m)
{
if(m.nRows==m.nColumns)
	{
	matrix <type> temp(m.nColumns,m.nRows);
	size_t i;
	for(i=0; i<m.nColumns; ++i)
		temp[i][i]=type(true);

	return temp;
	}
else
	throw matrix<type>::invalid_operator();
}

template <class type>
matrix<type> transposed(const matrix<type> &m)
{
matrix<type> temp(m.nColumns,m.nRows);
size_t i, j;

for(i=0; i<m.nColumns; i++)
	for(j=0; j<m.nRows; j++)
		temp[i][j]=m.mat[j][i];

return temp;
}

template <class type>
matrix<type> determinant(const matrix<type> &m)
{
if(m.nRows!=m.nColumns)
	throw matrix<type>::invalid_operator();

matrix<type> d(1,1);
//type d1;
size_t swaps=0;

if(m.nRows==1)
	d[0][0]=m.mat[0][0];
else if(m.nRows==2)
	d[0][0]=m.mat[0][0]*m.mat[1][1]-m.mat[0][1]*m.mat[1][0];
else
	{
	d=mul_diag(Gauss(m,&swaps));
	if(swaps%2)
		d[0][0]=-d[0][0];
	}
/*
matrix<type> temp(m);

bitarray C(temp.nColumns);
bitarray R(temp.nRows);
C=true;
R=true;
d1=temp.calculateDeterminant(C,R);
*/
return d;
}

template <class type>
matrix<type> plan(const matrix<type> &m)
{
if(m.nRows!=m.nColumns)
	{
	throw matrix<type>::invalid_operator();
	}

type sum(false);
for(size_t i=0; i<m.nRows; i++)
	sum+=m.mat[i][i];
return sum;
}

template <class type>
matrix<type> mul_diag(const matrix<type> &m)
{
if(m.nRows!=m.nColumns)
	{
	throw matrix<type>::invalid_operator();
	}

type product(true);
for(size_t i=0; i<m.nRows; i++)
	product*=m.mat[i][i];
return product;
}

template <class type>
matrix<type> sin(const matrix<type> &m)
{
type x;
binaryTree<function<type> > t=type::getFunctionsTree(x);

function<type> *pFunction;
pFunction=t.find(function<type>("sin",0));
if(!pFunction)
	{
	//throw incorrectFunction()
	}

return m.apply(pFunction->getFunction());
}

template <class type>
matrix<type> cos(const matrix<type> &m)
{
type x;
binaryTree<function<type> > t=type::getFunctionsTree(x);

function<type> *pFunction;
pFunction=t.find(function<type>("cos",0));
if(!pFunction)
	{
	//throw incorrectFunction()
	}

return m.apply(pFunction->getFunction());
}

template <class type>
matrix<type> tan(const matrix<type> &m)
{
type x;
binaryTree<function<type> > t=type::getFunctionsTree(x);

function<type> *pFunction;
pFunction=t.find(function<type>("tan",0));
if(!pFunction)
	{
	//throw incorrectFunction()
	}

return m.apply(pFunction->getFunction());
}

template <class type>
ostream &operator<<(ostream &os, const matrix<type> &m)
{
size_t i,j;
size_t ancho=os.width();
if(!ancho)
	{
	os.width(12);
	ancho=os.width();
	}
os<<'\n';

os<<'[';
os<<'\n';
for(i=0; i<m.nRows; i++)
	{
	for(j=0; j<m.nColumns; j++)
		{
		os.width(ancho);
		os<<m.mat[i][j]<<'\t';
		}
	os<<'\n';
	}
os<<']';
os<<'\n';

return os;
}

template <class type>
istream &operator>>(istream &is, matrix<type> &m)
{
size_t i,j;
deque<deque<type> > d;
deque<type>::size_type k=0, cols=0;
char matrixLine[4096], groupLine[4096];
char c;
int n;
type t;

if(is>>c && c=='[')
	{
	is.getline(matrixLine,4096,']');
	istrstream s(matrixLine,strlen(matrixLine));
	while(s)
		{
		if(s.getline(groupLine,4096))
			{
			istrstream g(groupLine,strlen(groupLine));
			d.push_back(deque<type>());
			g.eatwhite();
			while(g>>t)
				{
				d[k].push_back(t);
				g.eatwhite();
				}
			++k;
			}
		}
	for(int l=0; l<d.size(); ++l)
		if(d[l].size()>cols)
			cols=d[l].size();
	if(cols)
		{
		m=matrix<type>(d.size(),cols);
		for(i=0; i<d.size(); ++i)
			for(j=0; j<d[i].size(); ++j)
				{
				m.mat[i][j]=d[i][j];
				}
		}
	}
else
	{
	is.putback(c);

	if(is>>t)
		m=matrix<type>(t,m.pOstream);
	else
		is.clear(ios::failbit);
	}

return is;
}

}

#endif