#include <Origin.h>

void Orthogonalize_Active_Matrix()
{
	// Assumes the Matrix to be orthogonalized is active in the project.
	// Missing value in the Matrix is not checked in this function, 
	// so input all the data and set the matrix dimensions (select matrix > set dimensions in the menu)
	// exactly the real number of input data dimensions;
	
	//I have tested this using the data as following, it works well. 
	//Note that you should set the dimensions correctly, e.g. here it is 4 and 4.
	//
	//1		-1		0		0
	//1		1		-1		0
	//0		1		1		-1
	//0		0		1		1 
	
//the results is:
	//0.70711		-0.70711	0			0
	//0.57735		0.57735	 	-0.57735	0
	//0.31623		0.31623		0.63246		-0.63246
	//0.2582		0.2582		0.5164		0.7746

	//Get the Matrix data from the active matrix.
	MatrixLayer mlA;                           
	mlA = (MatrixLayer) Project.ActiveLayer(); 
	Matrix matA(mlA); 
	
	int	nRows = matA.GetNumRows();
	int	nCols = matA.GetNumCols();
	
	//create a new matrix page to store the orthogonalized matirx.
	MatrixPage mpB;
	mpB.Create("origin");
	if( !mpB.IsValid() )
	{
		printf("Create new matrix page failed!");
		return;
	}
	
	MatrixLayer mlB = mpB.Layers(0);
	Matrix matB(mlB);
	// Set the size as same as the original matrix
	matB.SetSize(nRows, nCols);
	
	//performing the Gram_Schmidt orthogonalization.
	Gram_Schmidt_Orthogonalization(matA, matB);
}

void Gram_Schmidt_Orthogonalization(const Matrix& mInput, Matrix& mOutput)
{
	int nRows = mInput.GetNumRows();
	
	for(int ii=0; ii < nRows; ii++)
	{
		vector vAlpha, vBeta;
		mInput.GetRow(vAlpha, ii);
		vBeta = vAlpha;
		vector vBetajj;
		double dAB, dBB; // dAB = vAlpha'*vBetajj; dBB = ||vBetajj||^2 = vBetajj'*vBetajj;
		
		//beta_i = alpha_i - sum_(j=0...i-1){( (alpha_i' * beta_j)/||beta_j||^2 ) * beta_j }; 
		for(int jj=0; jj<ii; jj++)
		{
			mOutput.GetRow(vBetajj, jj);
			dAB = vectors_inner_product(vAlpha, vBetajj);
			dBB = vectors_inner_product(vBetajj, vBetajj);
			vBeta = vBeta - (dAB / dBB) * vBetajj;
		}
		
		// normalize the vector
		double dBii;
		// dBii = ||vBeta_ii|| = sqrt( vBeta_ii*vBeta_ii');
		dBii = sqrt(vectors_inner_product(vBeta, vBeta));
		vBeta = vBeta / dBii;
		mOutput.SetRow(vBeta, ii);
	}
}

double vectors_inner_product(const vector vA, const vector vB)
{
	// this function calculates the inner product of two vectors, the result is a double value; 
	// if A=(a1, a2, ... an); B=(b1,b2,...bn); then  Sum=sum_(i=1...n){ai*bi}
	double dSum = 0;
	int nA = vA.GetSize();
	int nB = vB.GetSize();
	if(nA != nB) 
	{
		printf("the two vectors shoule have same size");
		return NANUM;		
	}
	for (int ii = 0; ii < nA; ii++)
	{
		dSum = dSum + vA[ii] * vB[ii];
	}
	return dSum;
}