#include <iostream>
#include <vector>
#include <stdlib.h>
#include <string>
#include <math.h> //fabs


/* TO BE DONE
 - need to check operations perform correctly
 */
 
/*create a matrix input string of form [1,2,3][4,5,6][7,8,9] - or allow commas between brackets, [],[],[] (they are ignored) or add B vector - [1,2,3|1][4,5,6|2][7,8,9|3]
take input in as string and convert it to a vector<T> A and vector<T> B inside of the mtx class

if the format works (implicitly) then the matrix should read it in - otherwise it should try a different form
*/

class line;

////////////////////////////////////////////////////////////////////////////////
//CLASS MTX DEFINITION
template <class T>
class mtx		
{
	private:
	friend class line;
	std::vector<T> A;	//matrix
	std::vector<T> B;	//solution vector - size = number of rows for matrix - not necessary for some matrix operations
	std::size_t width;	
	std::size_t height;	
	std::size_t x, y;
	//PRIVATE FUNCTIONS
	bool isMatrixSquare();	//handles an error
	void swapRows (int first_in, int second_in);
	class line
	{
		mtx * p_mtx;
		public:
		line(mtx *p_mtx_in);
		T& operator [] (std::size_t y_in);
		const T& operator [] (std::size_t y_in) const;
	};
	public:
	//CONSTRUCTORS
	mtx ();
	mtx (std::vector<T>& Ain, std::vector<T>& Bin);
	//INDEX ACCESS OPERATOR OVERLOADS
	line& operator [] (std::size_t index);
	T& operator () (std::size_t x_in, std::size_t y_in);
	const T& operator () (std::size_t x_in, std::size_t y_in) const;
	//MATRIX PUBLIC FUNCTIONS
	mtx operator = (mtx& in);	//ASSIGNMENT OPERATOR
	T& at (int i, int j);
	
	friend void swap(mtx<T>& first_in, mtx<T>& second_in);
	//UNARY FUNCTIONS
	mtx& operator += (const mtx& rhs);
	mtx& operator -= (const mtx& rhs);
	mtx& operator *= (const mtx& rhs);
	//BINARY FUNCTIONS
	mtx operator + (const mtx& rhs);
	mtx operator - (const mtx& rhs);
	mtx operator * (const mtx& rhs);
	mtx operator * (const T& rhs);	//scalar multiplication
	
	//void Apush_back (double);
	//double Aat (int row, int col);													//reads element at row,col where 0,0 is the first element in the matrix
	//void Bpush_back (double);
	//double Bat (int x);
	//bool 
	//T determinant();	//calculates the determinant of square matrix
	mtx inverse();	//returns the inverse matrix
	void rref ();																//Reduce Row Echelon Form - basic solution for square matrix
	void print();
	
	/*at some point more functions will be added but this should be enough for the beam calc
	void inverse (void);
	void eigen (void);
	mtx operator* (mtx const rhs&);	//A = B*C
	mtx& operator*= (mtx const rhs&);
	
	*/
	private:
	line m_line;
};

////////////////////////////////////////////////////////////////////////////////
//MATRIX CLASS DECLARATION
//PRIVATE FUNCTIONS*************************************************************
template <class T>
bool mtx<T>::isMatrixSquare()
{
	if (width == height)
	{
		return true;
	}
	std::cerr << "ERROR at bool mtx<T>::isMatrixSquare() - matrix is not square." << std::endl;
	return false;
}
template <class T>
void mtx<T>::swapRows (int row1_in, int row2_in)
{
	if (row1_in >= height || row2_in >= height)
	{
		std::cerr <<"Error - row parameters not within matrix and vector dimension";
		return;
	}
	T temp;
	for (int i; i < width; i++)
	{
		temp = A[i + row1_in * width];
		A[i + row1_in * width] = A[i + row2_in * width];
		A[i + row2_in * width] = temp;
	}
	temp = B.at(row1_in);	//swaps values in Bvector
	B.at(row1_in) = B.at(row2_in);
	B.at(row2_in) = temp;
	return;
}
template <class T>
T& mtx<T>::line::operator [] (std::size_t y_in)
{
	if (y_in >= p_mtx->height)
	{
		std::cerr << "row value is out of bounds and will be replaced with 0";
		p_mtx->y = 0;
	}
	else
	{
		p_mtx->y = y_in;
	}
	return p_mtx->A[p_mtx->x+p_mtx->y*p_mtx->width];
}
template <class T>
const T& mtx<T>::line::operator [] (std::size_t y_in) const
{
	if (y_in >= p_mtx->height)
	{
		std::cerr << "row value is out of bounds and will be replaced with 0";
		p_mtx->y = 0;
	}
	else
	{
		p_mtx->y = y_in;
	}
	return p_mtx->A[p_mtx->x+p_mtx->y*p_mtx->width];
}
//MATRIX CLASS CONSTRUCTORS*****************************************************
template <class T>
mtx<T>::mtx ():A(), B(), width(0), height(0), m_line(this)
{
	
}
template <class T>
mtx<T>::mtx (std::vector<T>& Ain, std::vector<T>& Bin):A(Ain), B(Bin), m_line(this)
{
	std::cout << "constructor started" <<std::endl;
	if(Ain.size() != 0 && Bin.size() != 0 && Ain.size() % Bin.size() == 0)
	{
		height = Bin.size();
		width = Ain.size()/height;
	}
	else
	{
		std::cerr << "Error - Matrix A and Vector B are of incompatible sizes - dummy matrix created";
		A.clear();
		B.clear();
	}
}
//MATRIX CLASS ACCESS OPERATOR OVERLOADS****************************************
template <class T>
typename mtx<T>::line& mtx<T>::operator [] (std::size_t x_in)
{
	if (x_in >= width)
	{
		std::cerr << "row coordinate is out of bounds - setting to 0";
		x = 0;
	}
	else
	{
		x = x_in;
	}
	return m_line;
}
template <class T>
T& mtx<T>::operator () (std::size_t x_in, std::size_t y_in)
{
	return A[x_in + y_in * width];
}
template <class T>
const T& mtx<T>::operator () (std::size_t x_in, std::size_t y_in) const
{
	return A[x_in + y_in * width];
}
//MATRIX CLASS - UNARY FUNCTION OVERLOADS
template <class T>
mtx<T>& mtx<T>::operator += (const mtx& rhs)
{
	if (height != rhs.height || width != rhs.width)
	{
		std::cerr << "ERROR at mtx& mtx<T>::operator += (const mtx& rhs) - matrix sizes do not match.  No addition performed." << std::endl;
		return *this;
	}
	for (int i = 0; i<= A.size(); i++)
	{
		A[i] += rhs.A[i];
	}
	return *this;
}
template <class T>
mtx<T>& mtx<T>::operator -= (const mtx& rhs)
{
	if (height != rhs.height || width != rhs.width)
	{
		std::cerr << "ERROR at mtx& mtx<T>::operator -= (const mtx& rhs) - matrix sizes do not match.  No subtraction performed." << std::endl;
		return *this;
	}
	for (int i = 0; i<= A.size(); i++)
	{
		A[i] -= rhs.A[i];
	}
	return *this;
}
template <class T>
mtx<T>& mtx<T>::operator *= (const mtx& rhs)
{
		if (width != rhs.height)
	{
		std::cout << "ERROR at mtx& mtx<T>::operator *= (const mtx& rhs) - matrix sizes are not compatible.  No multiplication performed." << std::endl;
		return *this;
	}
	//create new matrix then swap with this
	std::vector<T> Atemp;
	std::vector<T> Btemp(0,height);
	for (int j = 0; j < height; j++)//use j for height in this and width in rhs
	{
		for (int i = 0; i < rhs.width; i++)
		{
			T mtxValue = 0;
			for (int n = 0; n < width; n++)
			{
				mtxValue += A[n+j*width] * rhs(i,n); 
			}
			Atemp.push_back (mtxValue);
		}
	}
	mtx<T> temp(Atemp, Btemp);
	//swap(*this, temp);
	return temp;
}





template <class T>
T& mtx<T>::at (int x_in, int y_in)
{
	return this->A.at(x_in + y_in * width);
}



template <class T>
void mtx<T>::rref ()
{
	//WHAT HAPPENS WHEN WE HAVE MORE COLUMNS THAN ROWS?
	//MORE ROWS THAN COLUMNS?
	//WHEN A COLUMN DOESNT HAVE A 1 ON THE DIAGONAL
	if(A.size() == 0 || B.size() == 0 || A.size() % B.size() != 0)
	{
		std::cerr << "ERROR at [void mtx<T>::rref ()] - Matrix A and Vector B are of incompatible sizes - RREF not performed";
		return;
	}	
	T largest, coefficient;
	int largestRow, count;
	count = 0;
	for (int i = 0; i< width; i++)	//iterate through each column
	{
		largest = 0;
		for (int j = count; j < height; j++)	//find largest in each row
		{
			if (fabs(A[i + j * width]) > largest)
			{
				largest = fabs(A[i + j * width]);
				largestRow = j;
			}
		}
		if (largestRow != count)	//not on the diagonal then swap it onto the diagonal
		{
			this->swapRows(largestRow, count);
		}
		if (i != height && largest != 0)	//if largest is 0 then skip
		{
			for (int k = 0; k < height; k++)	//iterate through rows and 0 out column
			{
				if (k != count)	//skip over diagonal
				{
					coefficient = -A[i + k * width]/A[i + i * width];
					std::cout << "coefficient = " << coefficient << std::endl;
					A[i + k * width] = 0;	//first coefficient is 0 - no need to calc
					for (int m = i + 1; m < width; m++)	//iterate across width
					{
						A[m + k * width] += A[m + i * width] * coefficient; 
					}
					B[k] += B[i] * coefficient;	
				}
			}
			this->print();
		}
		if (A[i + i * width] != 0)	//set diagonal value to 1
		{
			coefficient = 1.0 / A[i + i * width];
			A[i + i * width] = 1.0;
			for (int n = i + 1; n < width; n++)
			{
				A[n + i * width] *= coefficient;
			}
			B[i] *= coefficient;
			count++;
			if (count == width || count == height)
			{
				return;
			}
		}
	}
}

template <class T>
mtx<T>::line::line (mtx * p_mtx_in):p_mtx(p_mtx_in)
{
	
}

template <class T>
void swap (mtx<T>& first_in, mtx<T>& second_in)
{
	using std::swap;
	swap(first_in.A, second_in.A);
	swap(first_in.B, second_in.B);
	swap(first_in.width, second_in.width);
	swap(first_in.height, second_in.height);
	swap(first_in.x, second_in.x);
	swap(first_in.y, second_in.y);
	swap(first_in.m_line, second_in.m_line);
	return;
}

template <class T>
void mtx<T>::print ()
{
	for (int j=0; j<height; j++)
	{
		for (int i = 0; i<width; i++)
		{
			std::cout << A[i + j * width] << " ";
		}
		std::cout << "| " << B[j] << std::endl;
	}
}

int main() {
	// your code goes here
	std::vector<double> atemp, btemp;
	
	for (int i = 0; i < 12; i++)
	{
		atemp.push_back(double(pow(i,3)/4.12785));
	}
	
	for (int i = 0; i < 3; i++)
	{
		btemp.push_back(double(4.0-i));
	}
	mtx<double> M (atemp, btemp);
	M.print();
	M *= M;
	M.rref();
	M.print();
	double pi = 3.14159;
	M [1][2] = pi;
	M(0,1) = pi;
	M.print();
	M.rref();
	M.print();
	
	return 0;
}