# include "Matrix.h"
#include "Global.h"

Matrix::Matrix():col(Global::g_dimNumber),row(Global::g_dimNumber){
	if(row==0) return;
	allocateMemory(row,col);

}
Matrix::Matrix(const int c, const int r):col(c),row(r){
	allocateMemory(row,col);

}
void Matrix::setDataRow(const double *d,const int c,const int r){
	if(r!=row) return;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=d[j];
}
void Matrix::setDataCol(const double *d, const int r,const int c){
	if(c!=col) return;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=d[i];

}
void Matrix::setData(const double * const * d){
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=d[i][j];
}
Matrix::~Matrix(){

    freeMemory();
}
void Matrix::operator *(const Matrix &m){

	if(col!=m.row) Throw(Logic_error("can not *, col of 1st matrix must be equal to row of 2nd matrix"));

	Matrix r(m.col,row);
	for(int i=0;i<row;i++){
		for(int j=0;j<m.col;j++){
			r.vec[i].data[j]=0;
			for(int k=0;k<col;k++)
				r.vec[i].data[j]+=vec[i].data[k]*m.vec[k].data[j];
		}
	}
	if(row!=m.row||col!=m.col){
       freeMemory();
	   row=r.row;col=r.col;
       allocateMemory(row,col);
	}
	*this=r;
}
void Matrix::operator *(const double x){
	for(int i=0;i<row;i++){
		for(int j=0;j<col;j++){
			vec[i].data[j]*=x;
		}
	}

}
Matrix&Matrix::operator =(const Matrix &m){
	if(row!=m.row||col!=m.col) return *this;
	if(this==&m) return *this;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=m.vec[i].data[j];
	return *this;
}
bool Matrix::Identity(){
	if(row!=col) return false;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]= (j==i);

	return true;
}
bool Matrix::isIdentity(){
    if(row!=col) return false;
    for(int i=0;i<row;i++){
        for(int j=0;j<col;j++) {
        if(vec[i].data[j]!=(i==j)){
                return false;
            }
        }
    }
    return true;
}
void Matrix::Diagonal(const double CondiNum ){
	if(row!=col) Throw(Logic_error("can not be diagonal, matrix must be squre"));

	double min,max;
	double *r=new double[row];
	for(int i=0;i<row;i++)	r[i]=gRandFloat(1,row,false);

	min=max=r[0];
	for(int i=0;i<row;i++)	{
		if(min>r[i])min=r[i];
		if(max<r[i]) max=r[i];
	}
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(j!=i) vec[i].data[j]=0.;
			else vec[i].data[j]=pow(CondiNum,(r[i]-min)/(max-min));
	delete []r;
}
void Matrix::InitialToZero(){
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=0.;
}
void Matrix::Set_Rotation(const int &r, const int &c,const double &angle){
	Identity();
	vec[r].data[r]=cos(angle);
	vec[r].data[c]=-sin(angle);
	vec[c].data[r]=sin(angle);
	vec[c].data[c]=cos(angle);
}


void Matrix::GenerateRotationMatrix(const double CondiNum,bool rMode){

	Matrix ortholeft,orthoright,diagonal;
	ortholeft.Randomize(rMode);
	ortholeft.Orthonormalization();
	orthoright.Randomize(rMode);
	orthoright.Orthonormalization();
	diagonal.Diagonal(CondiNum);
	ortholeft*diagonal;
	ortholeft*orthoright;
	*this=ortholeft;

}
void Matrix::Randomize(bool rMode){
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			vec[i].data[j]=gRandFloat(-1,1,rMode);
}
void Matrix::Orthonormalization(){
	if(row!=col) Throw(Logic_error("can not orthonormalization, matrix must be squre"));

	MyVector t1(col),t2(col);
	//Gram schmidt process
	for(int i=1;i<row;i++){
		t1=vec[i];
		t2=vec[i];
		for(int j=0;j<i;j++){
			t1.ProjectionToV(vec[j]);
			t2-t1;
			t1=vec[i];
		}
		vec[i]=t2;
	}
	// Normalization each vector
	for(int i=0;i<row;i++) vec[i].Normalization();
}
void Matrix::Transpose(){
	double t;
	for(int i=0;i<row;i++)
		for(int j=0;j<i;j++){
			t=vec[i].data[j];
			vec[i].data[j]=vec[j].data[i];
			vec[j].data[i]=t;
		}

}

void Matrix::inverse(){
	double **a;

	a=new double*[row];
	for(int i=0;i<row;i++ ) a[i]=new double[2*col];

	for(int i=0;i<row;i++){
        for(int j=0;j<col;j++) {
			a[i][j]=vec[i].data[j];
        }
    }

	for(int i=0;i<row;i++){
       for(int j=col;j<col*2;j++)
          a[i][j]=(j-i==col)?(1.0):(0.0);
    }

    for(int i=0;i<row;i++)
    {
       if(a[i][i]!=1.0)
       {
          double tmp=a[i][i];
          a[i][i]=1.0;
          for(int j=i;j<col*2;j++)
             a[i][j]/=tmp;
       }
       for(int k=0;k<row;k++)
       {
          if(i!=k)
          {
             double tmp=a[k][i];
             for(int j=0;j<col*2;j++)
                a[k][j]-=(tmp*a[i][j]);
          }
          else continue;
       }
    }

    for(int i=0;i<row;i++){
       for(int j=col;j<col*2;j++)
         vec[i].data[j-col]=a[i][j];
    }

	for(int i=0;i<row;i++) delete [] a[i];
	delete [] a;
}

void Matrix::freeMemory(){
    delete [] vec;
}
void Matrix::allocateMemory(const int r, const int c){
    vec=new MyVector[r];
	for(int i=0;i<r;i++)
		vec[i].allocateMemory(c);
}
double Matrix::getData(const int r, const int c){
    if(r<row&&c<col) return vec[r].data[c];
	else Throw(Out_of_range("out of range error in Matrix getData()"));
}
