mymatrix.cpp

PRINCIPAL COMPONENT ANALYSIS,using C++ language

#include "MyMatrix.h"
#include <cassert>

MyMatrix::MyMatrix(double matrix[][MAXSIZE]/* =NULL */,
				   int row/* =MAXSIZE */,int col/* =MAXSIZE */){
	this->col=col;
	this->row=row;
	if(matrix!=NULL){
		for(int i=0;i<row;i++)
			for(int j=0;j<col;j++)
				this->matrix[i][j]=matrix[i][j];
	}
	else{
		for(int i=0;i<row;i++){
			for(int j=0;j<col;j++){
				if(i!=j)	this->matrix[i][j]=0;
				else		this->matrix[i][j]=1;
			}
		}
	}
}

MyMatrix::MyMatrix(MyMatrix &A){
	col=A.col;
	row=A.row;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			matrix[i][j]=A.matrix[i][j];
}

void MyMatrix::operator =(MyMatrix B){
	col=B.col;
	row=B.row;
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			matrix[i][j]=B.matrix[i][j];
}

void MyMatrix::Display(){
	cout<<setiosflags(ios::fixed|ios::showpoint);
	for(int i=0;i<row;i++){
		for(int j=0;j<col;j++)
			cout<<setw(10)<<setprecision(5)<<matrix[i][j];
		cout<<endl;
	}
	cout<<endl;
}

double MyMatrix::GetMaxAbsNondialogElement(int &r,int &c){
	int i,j;
	double max=matrix[0][1];
	r=0;
	c=1;
	for(i=0;i<row;i++){
		for(j=0;j<i;j++){
			if(i==j)	continue;
			if(fabs(matrix[i][j])>fabs(max)){
				max=matrix[i][j];
				r=i;
				c=j;
			}
		}
	}

	return max;
}

double MyMatrix::GetMaxAbsElement(){
	double max=matrix[0][0];
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			if(fabs(matrix[i][j])>fabs(max))
				max=matrix[i][j];
	return max;
}

MyMatrix operator*(MyMatrix A,MyMatrix B){
	if(A.col!=B.row)
	{
		cerr<<"Error : invalid demension , matrix multiplication failed!"<<endl;
		return MyMatrix();
	}
	int i,j,k;
	double sum;
	MyMatrix product(NULL,A.row,B.col);

	for(i=0;i<A.row;i++){
		for(j=0;j<B.col;j++){
			sum=0;
			for(k=0;k<A.col;k++)
				sum+=A.matrix[i][k]*B.matrix[k][j];
			product.matrix[i][j]=sum;
		}
	}
	return product;
}

MyMatrix operator*(MyMatrix A,double coefficient){
	for(int i=0;i<A.row;i++)
		for(int j=0;j<A.col;j++)
			A.matrix[i][j]*=coefficient;
	return A;
}

MyMatrix MyMatrix::JacobiMethod(MyMatrix& eigenvalue){
	MyMatrix temp=(*this);
	eigenvalue=MyMatrix(NULL,this->row,this->col);

	double cosTheta,sinTheta,tan2Theta,sin2Theta,cos2Theta,max;
	int rpos,cpos;

	for(int i=0;i<TIMES;i++){
		MyMatrix P(NULL,row,col);
		max=temp.GetMaxAbsNondialogElement(rpos,cpos);
		if(fabs(temp.matrix[rpos][rpos]-temp.matrix[cpos][cpos])<10e-30){
			sinTheta=cosTheta=1.0/sqrt(2.0);
		}
		else{
			tan2Theta=2*temp.matrix[rpos][cpos]/(temp.matrix[rpos][rpos]-temp.matrix[cpos][cpos]);
			cos2Theta=1.0/sqrt(1+tan2Theta*tan2Theta);
			sin2Theta=tan2Theta*cos2Theta;
			cosTheta=sqrt((1+cos2Theta)/2);
			sinTheta=sin2Theta/(2*cosTheta);
		}
		P.matrix[rpos][rpos]=P.matrix[cpos][cpos]=cosTheta;
		P.matrix[rpos][cpos]=-sinTheta;
		P.matrix[cpos][rpos]=sinTheta;
		temp=P.TransposeMatrix()*temp*P;
		eigenvalue=eigenvalue*P;
	}
	return temp;
}

MyMatrix MyMatrix::TransposeMatrix(){
	MyMatrix temp(NULL,col,row);
	for(int i=0;i<row;i++)
		for(int j=0;j<col;j++)
			temp.matrix[j][i]=matrix[i][j];

	return temp;
}

double MyMatrix::GetElement(int row, int col){
	if(row>=this->row||col>=this->col){
		cerr<<"Error : illegal index , out of range"<<endl;
		return 0.0;
	}

	return matrix[row][col];
}

void MyMatrix::SwapColumn(int a, int b){
	if(a>=col||b>=col){
		cerr<<"Error : illegal index , out of range"<<endl;
		return;
	}

	for(int i=0;i<row;i++){
		double t=matrix[i][a];
		matrix[i][a]=matrix[i][b];
		matrix[i][b]=t;
	}
}

bool MyMatrix::SetElement(int row, int col, double value){
	if(row>=this->row||col>=this->col){
		cerr<<"Error : illegal index , out of range"<<endl;
		return false;
	}

	matrix[row][col]=value;
	return true;
}

int MyMatrix::GetRow(){
	return row;
}

int MyMatrix::GetCol(){
	return col;
}

double MyMatrix::GetMeanOfCol(int col){
	assert(col>=0&&col<this->col);

	double res=0;
	for(int i=0;i<row;i++){
		res+=matrix[i][col];
	}
	return res/row;
}

double MyMatrix::GetMeanOfRow(int row){
	assert(row>=0&&row<this->row);

	double res=0;
	for(int i=0;i<col;i++){
		res+=matrix[row][i];
	}

	return res/col;
}

MyMatrix MyMatrix::StandarlizeVectors(){
	MyMatrix res(NULL,row,col);
	double mean,Lii,t;
	for(int i=0;i<row;i++){
		Lii=0;
		mean=GetMeanOfRow(i);
		for(int k=0;k<col;k++)
			Lii+=(matrix[i][k]-mean)*(matrix[i][k]-mean);
		t=sqrt(Lii/(res.col-1));
		for(int j=0;j<col;j++){
			res.matrix[i][j]=(matrix[i][j]-mean)/t;
		}
	}

	return res;
}

double MyMatrix::GetRowStandardDeviration(int row){
	assert(row>=0&&row<this->row);
	double mean=GetMeanOfRow(row);

	double res=0;
	for(int i=0;i<col;i++){
		res+=(matrix[row][i]-mean)*(matrix[row][i]-mean);
	}
	res=sqrt(res/(col-1));

	return res;
}

double MyMatrix::GetRowCovirance(int r1, int r2){
	assert(r1>=0&&r1<this->row);
	assert(r2>=0&&r2<this->row);

	double res=0;
	double mean1=GetMeanOfRow(r1),mean2=GetMeanOfRow(r2);
	for(int i=0;i<col;i++){
		res+=(matrix[r1][i]-mean1)*(matrix[r2][i]-mean2);
	}
	res/=(col-1);

	return res;
}

MyMatrix MyMatrix::GetCorelationMatrix(){
	MyMatrix res(NULL,row,row);
	double *standardVariance=new double[row];
	int i,j;
	for(i=0;i<row;i++)
		standardVariance[i]=GetRowStandardDeviration(i);

	for(i=0;i<row;i++){
		for(j=0;j<i;j++){
			res.matrix[j][i]=res.matrix[i][j]=
				GetRowCovirance(i,j)/(standardVariance[i]*standardVariance[j]);
		}
	}

	delete []standardVariance;
	return res;
}