//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



// Template Implementation


template< class T >
Matrix<T>::Matrix()
	: ndims(2), length(0)
{
	data = 0;
	xDim = 0;
	yDim = 0;
	columns = 0;
	clean = 0;

}


template< class T >
Matrix<T>::Matrix(const int rows, const int cols)
{
	if(rows < 1 || cols < 1)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
	}

	ndims = 2;
	length = rows*cols;
	xDim = cols;
	yDim = rows;
	data = new (std::nothrow) T[length];
	Utility::CheckPointer(data);

	columns = new (std::nothrow) Vector<T>[cols];
	Utility::CheckPointer(columns);
	for(int i=0; i<cols; i++)
	{
		columns[i].Set(&(data[i*rows]), rows);
	}
	
	clean = new (std::nothrow) Cleaner<T>(data, columns);
	Utility::CheckPointer(clean);
	
}

template< class T >
Matrix<T>::Matrix(const int rows, const int cols, T init)
{
	if(rows < 1 || cols < 1)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
	}

	ndims = 2;
	length = rows*cols;
	xDim = cols;
	yDim = rows;
	data = new (std::nothrow) T[length];
	Utility::CheckPointer(data);

	columns = new (std::nothrow) Vector<T>[cols];
	Utility::CheckPointer(columns);
	for(int i=0; i<cols; i++)
	{
		columns[i].Set(&(data[i*rows]), rows);
	}
	
	clean = new (std::nothrow) Cleaner<T>(data, columns);
	Utility::CheckPointer(clean);
	
	for(int i=0;i<length;i++)
	{
		data[i] = init;
	}

}




template< class T >
Matrix<T>::Matrix(string str)
{
	if(str.substr(0,1) == "[" && str.substr(str.size()-1,1) == "]")
	{
		str = str.substr(1,str.size()-2);
		vector<string> row_strs = Utility::Split(str, ";");
		if(row_strs.size() < 1 )
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
		}
		int rsize = (int)row_strs.size();
		
		vector<string> dummy = Utility::Split(row_strs[0]);
		if(dummy.size() < 1 )
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
		}
		int csize = (int)dummy.size();
		ndims = 2;
		length = rsize*csize;
		xDim = csize;
		yDim = rsize;
		data = new (std::nothrow) T[length];
		Utility::CheckPointer(data);
		columns = new (std::nothrow) Vector<T>[csize];
		Utility::CheckPointer(columns);
		for(int i=0; i<csize; i++)
		{
			columns[i].Set(&(data[i*rsize]), rsize);
		}
		clean = new (std::nothrow) Cleaner<T>(data, columns);
		Utility::CheckPointer(clean);
		
		for(int i=0; i<rsize; i++)
		{
			vector<string> elem_strs = Utility::Split(row_strs[i]);
			if(elem_strs.size() != xDim )
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("String parse error. Matrix dimensions are not consistent!");
			}
			for(int j=0; j<csize; j++)
			{
				columns[j].data[i] = (T)Utility::ToDouble(elem_strs[j]);
			}
		}
	}
	else
	{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Incorrect initialization. String cannot be parsed!");
	}
}





template< class T >
Matrix<T>::Matrix(Matrix<T> &m) 
{
	ndims = m.ndims;
	length = m.length;
	xDim = m.xDim;
	yDim = m.yDim;

	data = m.data;
	columns = m.columns;

	clean = new (std::nothrow) Cleaner<T>(data, columns);
	Utility::CheckPointer(clean);

}



template< class T >
Matrix<T>::~Matrix()
{
	
	if(clean != 0)
	{
		delete clean;
	}
}


template< class T >
void Matrix<T>::Clean()
{
	if(clean != 0)
	{
		delete clean;
		data = 0;
		columns = 0;
		clean = 0;
	}
}


template< class T >
void Matrix<T>::Set(const int rows, const int cols, T *_data)
{
	if(rows < 1 || cols < 1)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
	}

	ndims = 2;
	length = rows*cols;
	xDim = cols;
	yDim = rows;
	data = _data;

	columns = new (std::nothrow) Vector<T>[cols];
	Utility::CheckPointer(columns);
	for(int i=0; i<cols; i++)
	{
		columns[i].Set(&(data[i*rows]), rows);
	}
	
	delete clean;
	clean = new (std::nothrow) Cleaner<T>(data, columns);
	Utility::CheckPointer(clean);
	
}

template< class T >
inline const T* Matrix<T>::DataPtr()
{
	return data;
}

template< class T >
inline T* Matrix<T>::Data()
{
	return data;
}


template< class T >
Matrix<T> Matrix<T>::Clone() const
{
	Matrix<T> temp(Rows(), Columns());
	
	memcpy(temp.data, data, sizeof(T)*temp.length);

	return temp;
}


template< class T >
void Matrix<T>::SwitchColumns(int i, int j)
{
	T *temp_i = columns[i].data;
	columns[i].data = columns[j].data;
	columns[j].data = temp_i;
}


template< class T >
void Matrix<T>::SyncData2Columns()
{
	T* temp_data = new (std::nothrow) T[length];
	Utility::CheckPointer(temp_data);
	
	for(int i=0; i<xDim; i++)
	{
		memcpy(&(temp_data[i*yDim]), columns[i].data, sizeof(T)*yDim);
	}
	data = temp_data;
	
	columns = new (std::nothrow) Vector<T>[xDim];
	Utility::CheckPointer(columns);
	
	for(int i=0; i<xDim; i++)
	{
		columns[i].data = &(data[i*yDim]);
	}
	
	delete clean;
	clean = new (std::nothrow) Cleaner<T>(data, columns);
	Utility::CheckPointer(clean);
}


template< class T >
Matrix<T> Matrix<T>::Slice(const int row1, const int row2, const int col1, const int col2)
{
	if(row1<0 || row1>=Rows() || row2<0 || row2>=Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Index outside bounds!");
	}
	if(row1 > row2)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
	}


	if(col1<0 || col1>=Columns() || col2<0 || col2>=Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Index outside bounds!");
	}
	if(col1 > col2)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
	}

	
	Matrix<T> temp(row2-row1+1, col2-col1+1);
	
	for(int j=col1; j<=col2; j++)
	{
		memcpy(temp[ j-col1 ].data, &(columns[j].data[row1]), sizeof(T)*temp.Rows());
	}

	return temp;
}


//template< class T >
//SubMatrix<T> Matrix<T>::Slice(int row1, int row2, int col1, int col2)
//{
//	if(row1<0 || row1>=YDim() || row2<0 || row2>=YDim())
//	{
//		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//		Utility::RunTimeError("Index outside bounds!");
//	}
//	if(row1 > row2)
//	{
//		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//		Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
//	}
//
//
//	if(col1<0 || col1>=XDim() || col2<0 || col2>=XDim())
//	{
//		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//		Utility::RunTimeError("Index outside bounds!");
//	}
//	if(col1 > col2)
//	{
//		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//		Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
//	}
//
//	
//	SubMatrix<T> temp;
//	temp.xDim = col2-col1+1;
//	temp.yDim = row2-row1+1;
//	temp.columns = new (std::nothrow) Vector<T>[temp.xDim];
//	Utility::CheckPointer(temp.columns);
//	temp.clean = new (std::nothrow) Cleaner<T>(temp.columns);
//	Utility::CheckPointer(temp.clean);
//
//	for(int j=col1; j<=col2; j++)
//	{
//		temp.columns[j-col1].Set(&(columns[j].data[row1]),  row2-row1+1);
//	}
//
//	return temp;
//}



template< class T >
Matrix<T> Matrix<T>::Slice(string str1, string str2)
{
	int row1, row2, col1, col2;
	vector<string> bounds1 = Utility::Split(str1, ":");
	vector<string> bounds2 = Utility::Split(str2, ":");
	
	if(str1 == ":")
	{
		row1 = 0;
		row2 = yDim-1;
	}
	else if(bounds1.size() == 2)
	{
		row1 = Utility::ToInt(bounds1[0]);
		row2 = Utility::ToInt(bounds1[1]);
	}
	else
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Incorrect slice argument. String cannot be parsed!");
	}

	if(str2 == ":")
	{
		col1 = 0;
		col2 = xDim-1;
	}
	else if(bounds2.size() == 2)
	{
		col1 = Utility::ToInt(bounds2[0]);
		col2 = Utility::ToInt(bounds2[1]);
	}
	else
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Incorrect slice argument. String cannot be parsed!");
	}
	
	return this->Slice(row1, row2, col1, col2);

}



template< class T >
Vector<T> Matrix<T>::Slice(string str)
{
	Vector<T> temp = (Vector<T>)*this;
	return temp.Slice(str);
}





template< class T >
inline const int Matrix<T>::Columns() const
{
	return xDim;
}

template< class T >
inline const int Matrix<T>::Rows() const
{
	return yDim;
}


template< class T >
inline const int Matrix<T>::XDim() const
{
	return xDim;
}

template< class T >
inline const int Matrix<T>::YDim() const
{
	return yDim;
}




template< class T >
inline const int Matrix<T>::Length() const
{
	return length;
}

template< class T >
inline const int Matrix<T>::Numel() const
{
	return length;
}

template< class T >
inline const int Matrix<T>::NDims() const
{
	return ndims;
}

template< class T >
inline void Matrix<T>::Init(const T init)
{
	for(int i=0;i<length;i++)
	{
		data[i] = init;
	}
}

//template< class T >
//inline Matrix<T>&  Matrix<T>::Rand()
//{
//	if(!RandomGen::Initialized())
//	{
//		RandomGen::Initialize();
//	}
//	for(int i=0;i<length;i++)
//	{
//		data[i] = (T)rand();
//	}
//	return *this;
//}

template< class T >
inline Matrix<T>&  Matrix<T>::Rand(const double max)
{
	if(!RandomGen::Initialized())
	{
		RandomGen::Initialize();
	}
	for(int i=0;i<length;i++)
	{
		data[i] = (T)(rand()/(double)RAND_MAX*max);
	}
	return *this;
}


//template< class T >
//Matrix<T> Matrix<T>::Rand(const int rows, const int cols)
//{
//	Matrix<T> m(rows, cols);
//	if(!RandomGen::Initialized())
//	{
//		RandomGen::Initialize();
//	}
//	for(int i=0;i<rows*cols;i++)
//	{
//		m.data[i] = (T)rand();
//	}
//	return m;
//}

template< class T >
Matrix<T> Matrix<T>::Rand(const int rows, const int cols, const double max)
{
	Matrix<T> m(rows, cols);
	if(!RandomGen::Initialized())
	{
		RandomGen::Initialize();
	}
	for(int i=0;i<rows*cols;i++)
	{
		m.data[i] = (T)(rand()/(double)RAND_MAX*max);
	}
	return m;
}



template< class T >
void Matrix<T>::ReadFromMatrix(const Matrix<T>& m)
{
	if(xDim < m.xDim || yDim < m.yDim)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("You can't read from a larger Matrix to a smaller Matrix!");
	}

	for(int i=0;i<m.xDim;i++)
	{
		for(int j=0;j<m.yDim;j++)
		{
			columns[i].data[j] = m.columns[i].data[j];
		}
	}

}

template< class T >
void Matrix<T>::ReadFromMatrix(const Matrix<T>& m, const int rowStart, const int colStart)
{
	if(xDim < m.xDim+colStart || yDim < m.yDim+rowStart)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("You can't read from a larger Matrix (from the copy start point) to a smaller Matrix!");
	}

	for(int i=0;i<m.xDim;i++)
	{
		for(int j=0;j<m.yDim;j++)
		{
			columns[i+colStart].data[j+rowStart] = m.columns[i].data[j];
		}
	}

}









template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Matrix<T>& m1, Matrix<T>& m2)
{
	Matrix<T> temp;
	if(dimension == 1)
	{
		if(m1.Columns() != m2.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix sizes are not compatible for concatenation!");
		}
		
		temp = Matrix<T>(m1.Rows()+m2.Rows(), m1.Columns());
		temp.ReadFromMatrix(m1);
		temp.ReadFromMatrix(m2, m1.Rows(), 0);
		
	}
	else if(dimension == 2)
	{
		if(m1.Rows() != m2.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix sizes are not compatible for concatenation!");
		}

		temp = Matrix<T>(m1.Rows(), m1.Columns()+m2.Columns());
		temp.ReadFromMatrix(m1);
		temp.ReadFromMatrix(m2, 0, m1.Columns());
	}
	else
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
	}
	
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Matrix<T>& m1, Vector<T>& m2)
{
	return Matrix<T>::Cat(dimension, m1, (Matrix<T>)m2);
}

template< class T >