#ifndef __Matrix_CPP
#define __Matrix_CPP

#include "matrix.h"
#include <cstdlib>
#include<cstdio>
#include "Eigen/Dense"
#include<time.h>
#include <typeinfo>
#include <iostream>

#define MAX_RANDOM 200
const double eps = 1e-12;
// Parameter Constructor
using namespace std;


template<typename T>
MyMatrix<T>::MyMatrix(unsigned _rows, unsigned _cols, const T& _initial) {
  mat.resize(_rows);
  for (unsigned i=0; i<mat.size(); i++) {
    mat[i].resize(_cols, _initial);
  }

  rows = _rows;
  cols = _cols;
}

// Copy Constructor
template<typename T>
MyMatrix<T>::MyMatrix(const MyMatrix<T>& rhs) {
  mat = rhs.mat;
  rows = rhs.get_rows();
  cols = rhs.get_cols();
}

// (Virtual) Destructor
template<typename T>
MyMatrix<T>::~MyMatrix() {}

// Assignment Operator
template<typename T>
MyMatrix<T>& MyMatrix<T>::operator=(const MyMatrix<T>& rhs) {
  if (&rhs == this)
    return *this;

  unsigned new_rows = rhs.get_rows();
  unsigned new_cols = rhs.get_cols();

  mat.resize(new_rows);
  for (unsigned i=0; i<mat.size(); i++) {
    mat[i].resize(new_cols);
  }

  for (unsigned i=0; i<new_rows; i++) {
    for (unsigned j=0; j<new_cols; j++) {
      mat[i][j] = rhs(i, j);
    }
  }
  rows = new_rows;
  cols = new_cols;

  return *this;
}

// Addition of two matrices
template<typename T>
MyMatrix<T> MyMatrix<T>::operator+(const MyMatrix<T>& rhs) {
  Fraction f(1,1);
  MyMatrix<T> result(rows, cols, f);

  for (unsigned i=0; i<rows; i++) {
    for (unsigned j=0; j<cols; j++) {
      result(i,j) = this->mat[i][j] + rhs(i,j);
    }
  }

  return result;
}




// Left multiplication of this MyMatrix and another
template<typename T>
MyMatrix<T> MyMatrix<T>::operator*(const MyMatrix<T>& rhs) {
  unsigned rows = rhs.get_rows();
  unsigned cols = rhs.get_cols();
  Fraction f(0,1);
  MyMatrix<T> result(rows, cols, f);

  for (unsigned i=0; i<rows; i++) {
    for (unsigned j=0; j<cols; j++) {
      for (unsigned k=0; k<rows; k++) {
        result(i,j) += this->mat[i][k] * rhs(k,j);
      }
    }
  }

  return result;
}

// Multiply a MyMatrix with a vector
template<typename T>
std::vector<T> MyMatrix<T>::operator*(const std::vector<T>& rhs) {
  std::vector<T> result(rhs.size());

  for (unsigned i=0; i<rows; i++) {
    for (unsigned j=0; j<cols; j++) {
      result[i] += this->mat[i][j] * rhs[j];
    }
  }

  return result;
}



// Access the individual elements
template<typename T>
T& MyMatrix<T>::operator()(const unsigned& row, const unsigned& col) {
  return this->mat[row][col];
}

// Access the individual elements (const)
template<typename T>
const T& MyMatrix<T>::operator()(const unsigned& row, const unsigned& col) const {
  return this->mat[row][col];
}

// Get the number of rows of the MyMatrix
template<typename T>
unsigned MyMatrix<T>::get_rows() const {
  return this->rows;
}

// Get the number of columns of the MyMatrix
template<typename T>
unsigned MyMatrix<T>::get_cols() const {
  return this->cols;
}

template<typename T>
MyMatrix<T> MyMatrix<T>::swap_rows(unsigned row1, unsigned row2) {
  for (int k=0; k<mat.size(); k++) {
      double temp = this->mat[row1][k];
      this->mat[row1][k] = this->mat[row2][k];
      this->mat[row2][k] = temp;
  }

  return *this;
}

template<typename T>
void MyMatrix<T>::print() {
  for (int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
       cout << this->mat[i][j] << " ";
    }
    printf("\n");
  }
}

template<typename T>
Eigen::MatrixXd MyMatrix<T>::create_double_eigen() {
  Eigen::MatrixXd matE(rows, cols);

  for(int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
      matE(i,j) = (double) this->mat[i][j];
    }
  }

  return matE;
}

template<typename T>
Eigen::MatrixXf MyMatrix<T>::create_float_eigen() {
  Eigen::MatrixXf matE(rows, cols);

  for(int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
      matE(i,j) = this->mat[i][j];
    }
  }

  return matE;
}

template<typename T>
void MyMatrix<T>::gwp(std::vector<T> vec) {
  int i,j,k;
  double m,s,t;
  double X[rows], AB[rows][cols+1];
  Fraction f(1,1);

  for(i = 0; i < rows; i++) {
     AB[i][rows] = (double) vec[i] ;
  }

  for(int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
      AB[i][j] = (double) this->mat[i][j];
    }
  }

  for(i=0; i<rows-1; i++) {
    for(j=i+1; j<rows; j++) {
      if(fabs(AB[i][i]) < eps)
        return;
      m=-AB[j][i]/AB[i][i];
      for(k=i+1; k<=rows; k++)
        AB[j][k]+= m*AB[i][k];
    }
  }

  for(i=rows-1; i>=0; i--) {
    s = AB[i][rows];
    for(j=rows-1; j>=i+1; j--)
      s-=AB[i][j]*X[j];
    if(fabs(AB[i][i]) < eps)
      return;
    X[i]=s/AB[i][i];
  }

  for(i = 0; i < rows; i++)
    printf("%2.4f\n",X[i]);
}

template<typename T>
void MyMatrix<T>::gpp(const std::vector<T>& vec) {
  int i,j,k;
  double X[rows], AB[rows][cols+1];

  for(i = 0; i < rows; i++) {
    AB[i][rows] = (double) vec[i];
  }

  for(int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
      AB[i][j] = (double) this->mat[i][j];
    }
  }

  for(i=0; i<rows; i++)                    //Pivotisation
    for(k=i+1; k<rows; k++)
      if(abs(AB[i][i]) < abs(AB[k][i])) {
        for(j=0; j<=rows; j++){
          double temp=AB[i][j];
          AB[i][j]=AB[k][j];
          AB[k][j]=temp;
        }
      }

  for(i=0; i<rows-1; i++)            //loop to perform the gauss elimination
      for (k=i+1;k<rows;k++) {
            double t=AB[k][i]/AB[i][i];
            for (j=0;j<=rows;j++)
                AB[k][j]=AB[k][j]-t*AB[i][j];    //make the elements below the pivot elements equal to zero or elimnate the variables
          }

  for (i=rows-1;i>=0;i--)                //back-substitution
  {                        //x is an array whose values correspond to the values of x,y,z..
      X[i]=AB[i][rows];                //make the variable to be calculated equal to the rhs of the last equation
      for (j=i+1;j<rows;j++)
          if (j!=i)            //then subtract all the lhs values except the coefficient of the variable whose value                                   is being calculated
              X[i]=X[i]-AB[i][j]*X[j];
      X[i]=X[i]/AB[i][i];            //now finally divide the rhs by the coefficient of the variable to be calculated
  }

  for(i = 0; i < rows; i++)
    printf("%2.4f\n",X[i]);
}

template<typename T>
void MyMatrix<T>::gcp(const std::vector<T>& vec) {
  int n = rows, MAT1 = rows;
  int lambda[n], omega[n], p[rows], q[rows];
  int i,j,k;
  int pi,pj,tmp;
  double max;
  double ftmp,temp;
  double x[rows], AB[rows][cols+1];
  pi=0,pj=0,max=0.0;

  for(i = 0; i < rows; i++) {
    AB[i][rows] = vec[i];
  }

  for(int i=0; i<rows; i++) {
    for (int j=0; j<cols; j++) {
      AB[i][j] = this->mat[i][j];
    }
  }
  puts("Before pivot");
  for(i=0;i<rows;i++) {
    for(j=0;j<cols+1;j++) {
      printf("%2.4f ",AB[i][j]);
    }
    puts("");
  }

  for (k=0;k<n;k++){
      for (i=k;i<n;i++) {
          for (j=k;j<n;j++) {
              temp = fabs(AB[i][j]);
              if (temp>max){
                  max = temp;
                  pi=i;
                  pj=j;
              }
          }
      }

      tmp=p[k];
      p[k]=p[pi];
      p[pi]=tmp;
      for(i=0; i<n; i++) {
        ftmp = AB[k][i];
        AB[k][i] = AB[pi][i];
        AB[pi][i] = ftmp;
      }

      //Swap Col
      tmp=q[k];
      q[k]=q[pj];
      q[pj]=tmp;
      for (i=0;i<n;i++){
          ftmp=AB[i][k];
          AB[i][k]=AB[i][pj];
          AB[i][pj]=ftmp;
      }
      //END PIVOT

      //check pivot size and decompose
      if (AB[k][k]!=0){
          for (i=k+1;i<n;i++){
              //Column normalisation
              AB[i][k] = AB[i][k] / AB[k][k];
          }
          for (j=k+1;j<n;j++){
              for(i=k+1; i<n; i++){
                AB[i][j] = AB[i][j] - AB[i][k] * AB[k][j];
              }
          }
      }
      //END DECOMPOSE
  }
  puts("After pivot");
  for(i=0;i<rows;i++) {
    for(j=0;j<cols+1;j++) {
      printf("%2.4f ",AB[i][j]);
    }
    puts("");
  }
  /*
  int ii=0,ip;
  double xtmp[MAT1];
  double c[3];
  //Swap rows (x=Px)
  for (i=0; i<MAT1; i++){
      xtmp[i]=x[p[i]]; //value that should be here
  }
  //Lx=x
  for (i=0;i<MAT1;i++){
      ftmp=xtmp[i];
      for (j=ii-1;j<i;j++)
          ftmp-=AB[i][j]*xtmp[j];
      xtmp[i]=ftmp;
  }

  //xtmp[MAT1-1]/=a[MAT1-1][MAT1-1];
  for (i=MAT1-2;i>=0;i--){
      ftmp=xtmp[i];
      for (j=i+1;j<MAT1;j++){
          ftmp-=AB[i][j]*xtmp[j];
      }
      xtmp[i]=(ftmp)/AB[i][i];
  }

  for (i=0;i<MAT1;i++){
      x[q[i]]=xtmp[i];
  }

  for(i = 0; i < rows; i++)
    printf("%2.4f\n",x[i]); */
}
#endif