#ifndef _matrix_n_
#define _matrix_n_

#include "iostream"
#include <vector>
//#include <random>

#include "./utils/vector.h"


namespace utils{
	//#######################################
	//#            MATRIX TYPE              #
	//#######################################
	template <typename T>
	struct Matrix{
		utils::Vector<T> m;


		using vector_type = typename decltype(std::declval<T>().v)::value_type;

	   // Constructor to initialize matrix with rows × cols and a fill value
	    Matrix(uint64_t rows, uint64_t cols, typename T::value_type value = {}) {
	        fill(rows, cols, value);
	    }

	    Matrix() = default; // Default constructor

		T& operator[](uint64_t idx) { return m[idx]; } // Makes it able to do matr[1][1]
		const T& operator[](uint64_t idx) const { return m[idx]; } // Makes it able to do matr[1][1]


		void fill(uint64_t rows, uint64_t cols, const vector_type num=0){
			m.clear();
			for (uint64_t i = 0; i < rows; i++){
				T temp_vec;

				for (uint64_t j = 0; j < cols; j++){
					temp_vec.v.push_back(num);
				}
				m.push_back(temp_vec);
			}
		}

		void fill_RNG(const uint64_t rows, const uint64_t cols, const vector_type min = 0, const vector_type max = 1){
				m.clear();

				std::mt19937_64 rng{}; 
    			rng.seed( std::random_device{}());
    
				for (uint64_t i = 0; i < rows; i++){
					T temp_vec;
					for (uint64_t j = 0; j < cols; j++){
						temp_vec.v.push_back(std::uniform_real_distribution<>{min, max}(rng));
					}
					m.push_back(temp_vec);
				}
		}


		inline friend std::ostream& operator << (std::ostream& out, const Matrix& mat){
			out << "[";
			for (uint64_t i = 0; i < mat.m.size(); i++){
				out << "[";
				for (uint64_t j = 0; j < mat.m[i].v.size(); j++){
					if (j % mat.m[i].v.size() == mat.m[i].v.size() -1 && i == mat.m.size()-1){
						out << mat.m[i].v[j] << "]";
					}
					else if ((j % mat.m[i].v.size() == mat.m[i].v.size() -1)){
						out << mat.m[i].v[j] << "]," << std::endl;
					}
					else{
						out << mat.m[i].v[j] << ", ";
					}
				}
			}
			out << "]";
			return out;
		}
		void print() const{
			std::cout << *this << std::endl;
		}

		void inplace_transpose(){
			utils::Vector<T> temp_m = m;
			m.clear();

			uint64_t rows = temp_m.size();
			uint64_t cols = temp_m[0].v.size();

			for (uint64_t i = 0; i < cols; i++){
				T temp_vec;
				for (uint64_t j = 0; j < rows; j++){
					temp_vec.v.push_back(temp_m[j].v[i]);
				}
				m.push_back(temp_vec);
			}
	    }
 		Matrix<T> transpose()const{
 			Matrix<T> copy = *this;
 			copy.inplace_transpose();
 			return copy;
 		}


	   	void inplace_inverse(std::string method = "Gauss-Jordan"){
	   		//Matrix<T> temp_m = *this; // Copies the m into temp_m correctly (Before: utils::Vector<T> temp_m = m;)
	   		if (method == "Gauss-Jordan"){
	   			Matrix<T> temp_m(m.v.size(),m[0].v.size(),0);

	   			//std::cout << temp_m.m.v[0].size() << std::endl;
	   			//std::cout << m.v.size() << std::endl;

	   			uint64_t icol,irow,N=m.v.size(),M=temp_m.m.v[0].size();
	   			double big,dum,pivinv;
	   			Vi indxc(N,0),indxr(N,0),ipiv(N,0);

	   			//for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;}
   				for (uint64_t i = 0; i < N; i++){
   					big=0.0;
   					for (uint64_t j = 0; j < N; j++){
   						if (ipiv[j] != 1){
   							for (uint64_t k = 0; k < N; k++){
   								if (ipiv[k] == 0){
   									if (abs(m[j].v[k]) >= big){
   										big = abs(m[j].v[k]);
   										irow = j;
   										icol = k;
   									}
   								}
   							}
   						}
   					}
   					ipiv[icol]++;
   					if (irow != icol){
   						for (uint64_t l = 0; l < N; l++){ // SWAP 
   							double temp = m[irow].v[l];
   							m[irow].v[l] = m[icol].v[l];
   							m[icol].v[l] = temp;
   						}
   						for (uint64_t l = 0; l < M; l++){ // SWAP temp matrix
   							double temp = temp_m.m[irow].v[l];
   							temp_m.m[irow].v[l] = temp_m.m[icol].v[l];
   							temp_m.m[icol].v[l] = temp;
   						}
   					}

   					indxr[i] = irow;
   					indxc[i] = icol;
   					if (m[icol].v[icol] == 0.0){
   						throw std::runtime_error("utill:Matrix.Gauss-Jordan - Singular Matrix");
   					}
   					pivinv= 1.0/m[icol].v[icol];
   					m[icol].v[icol]=1.0;
   					for (uint64_t l = 0; l < N; l++){
   						m[icol].v[l] *= pivinv;
   					}
   					for (uint64_t l = 0; l < M; l++){
   						temp_m.m[icol].v[l] *= pivinv;
   					}
   					for (uint64_t ll = 0; ll < N; ll++){
   						if (ll != icol){
   							dum = m[ll].v[icol];
   							m[ll].v[icol] = 0;
   							for (uint64_t l = 0; l < N; l++){
   								m[ll].v[l] -= m[icol].v[l]*dum;
   							}
   							for (uint64_t l = 0; l < N; l++){
   								temp_m.m[ll].v[l] -= temp_m.m[icol].v[l]*dum;
   							}
   						}
   					}

   				}
   				//m = temp_m;
   				for (int64_t l = N-1; l >= 0; l--){
   					if (indxr[l] != indxc[l]){
   						for (uint64_t k = 0; k < N; k++){
   							double temp = m[k].v[indxr[l]];
							m[k].v[indxr[l]] = m[k].v[indxc[l]];
							m[k].v[indxc[l]] = temp;
   						}
   					}
   				}
	   		}
			else{
				throw std::runtime_error("utill:Matrix." + method + " - Not implemented yet \r \nImplemented: 'Gauss-Jordan',");
			}
	    }
		Matrix<T> inverse(std::string method = "Gauss-Jordan")const{
			Matrix<T> copy = *this;
			copy.inplace_inverse(method);
			return copy;
		}

	    utils::Vector<vector_type> vecmult(const utils::Vector<vector_type>& Vec)const{

	    	if (m[0].size() != Vec.size()){
	    		throw std::runtime_error("utill:Matrix.vecmult - Dimentions does not fit");
	    	}

		    // Create a temporary result vector
	    	utils::Vector<vector_type> copy(Vec.size(), 0);

		    
		    for (uint64_t i = 0; i < m.size(); ++i) {
		        for (uint64_t j = 0; j < m[0].size(); ++j) {
		            copy[i] += m[i][j] * Vec[j];
		        }
		    }
	    	return copy;
	    }


	    void inplace_matmult(const Matrix<T>& Mat){

	    	if (m.v[0].size() != Mat.m.v.size()){
	    		throw std::runtime_error("utill:Matrix.matmult - Dimentions does not fit");
	    	}

		    // Dimensions of the result
		    uint64_t rows = m.v.size();              // rows in *this
		    uint64_t cols = Mat.m[0].v.size();       // columns in Mat
		    uint64_t inner = m.v[0].size();          // shared dimension

		    // Create a temporary result matrix
	    	Matrix<T> temp_m(rows, cols, 0);

	    	// Perform matrix multiplication
	    	for (uint64_t i = 0; i < rows; i++){
	    		for (uint64_t j = 0; j < cols; j++){
	    			for (uint64_t k = 0; k < inner; k++){
	    				temp_m.m[i].v[j] += m[i].v[k] * Mat.m[k].v[j];
	    			}
	    		}
	    	}
	    	*this = temp_m;
	    }
	    Matrix<T> matmult(const Matrix<T>& Mat)const{
	    	Matrix<T> copy = *this;
	    	copy.inplace_matmult(Mat);
	    	return copy;
	    }


	};
	typedef Matrix<Vi> Mi;
	typedef Matrix<Vf> Mf;
	typedef Matrix<Vd> Md;
	
}


#endif // _numerics_n_