pca in c

Expand Post »
could anyone help me with this code?
the algo is alroght when tested manually. guess there is some problem with the declaration part which gives huge values(garbage )
thank you
C Syntax (Toggle Plain Text)
#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>
 
#define m 10   // m samples 
#define n 5   // n variables
#define dim 5 //whatever var dim is ***dim of covariance matrix
#define pertar 0.99
 
double scale_esti[m][n],mul[m][n];
double lambda,max,mu,w,sintheta,costheta,totvar,pervar,pervar1,temp,tmp;
double a[dim][dim],b[dim][dim],v[dim][dim],c[dim][dim];
double ix,im[10]={0},siginv[dim][dim]={0},det=1;// do not disturb siginv
double tousqr[1][1],zedy[1][5],evectsort[n][n];
double tl[m][2],tls[m][2],transload[2][5],tlstl[m][n];
double y[1][n],ytrans[n][1],v1[n],var[n],std[n];
double cov[n][n],cov1[n][n],sigma[dim][dim],t[m][n],mat[m][n];// do not disturb sigma
double mean[n],mu11[n],sel_evect[n][5];// do not disturb sel_evect
double combi[n+1][n],zed[m][n],zedytrans[n][1],spe[m][n];// spe should be 1,1;
double mattrans[2][5];
int i,h,j,z,edim,k,p,q,myflag=0,count=0;
 
double data[10][5]={{92,99,1,8,15},{98,80,7,14,16},{4,81,88,20,22},{85,87,19,21,3},{86,93,25,2,9},{17,24,76,83,90},{23,5,82,89,91},{79,6,13,95,97},{10,	12,	94,	96,	78},{11,	18,	100,	77,	84}};
 
FILE *val,*vec,*tou,*in,*scores,*sperror;
// functions
void   print(double q[][n],int,int);
void   newline();
void   print1(double v[],int);
double calmean(double q[][n],int ,int);
double calstd(double q[][n],int,int);
double autoscale(int, int, double q[][n],double r[],double s[]);
double covariance(int row,int col,double w[][n] );//,double u[n]
double evd(double r[][n]);
double combine(int dimn, double u[][n],double v[][n]);
double sort(int dimn,double w[][n]);
double loading(int dimn,double w[][n],double t[][n]);// w is a[][],z is evectsort[][]
double fscores();
double multiply(int g,int h,int r,double w[][n],double l[][n]);//multiply two matrices
double matinverse(int r,double w[][n]);
double computesigma(int dimn,double w[][n]);
double sigmainverse(int);
double ftousqr();
double mattranspose(int p,int k,double w[][n]);
double estimate(double w[][2],double r[][dim]);
double fspe();
 
main()
{
	val=fopen("eigen_val3.dat","w");
	vec=fopen("eigen_vect3.dat","w");
	tou=fopen("tousqr3.xls","w");
	scores=fopen("scores3.xls","w");
	sperror=fopen("spe.xls","w");
	scores=fopen("scores3.xls","r");
 
	print(data,m,n);
getch();
	calmean(data,m,n);
	print1(mean,n);
	newline();
getch();
	calstd(data,m,n);
	print1(std,n);
	newline();
getch();
	autoscale(m,n,data,mean,std);
	print(mat,m,n);
	newline();
getch();
	calmean(mat,m,n);
	print1(mean,n);
	newline();
getch();
	covariance(10,5,mat);
	print(cov,n,n);
	newline();
getch();
	printf("evd of covariance\n");
	evd(cov);
getch();
	printf("eigen value matrix\n");
	print(a,n,n);
	newline();
getch();
	printf("eigen vector matrix\n");
	print(c,n,n);
	newline();
getch();
	printf("\ncombined matrix >>>>>>>>>>>\n");
	combine(n,a,c);
	print(combi,n+1,n);
	newline();
getch();
	printf("\ncombined matrix sorted>>>>>>>>>>>\n");
	sort(n,combi);
	newline();
	printf("*********************sort done\n");
getch();
	loading(n,a,evectsort);
	newline();
	print(sel_evect,dim,edim);
	newline();
	printf("*******************selection done\n");
getch();
	fscores();
	printf("*******************scores done\n");
	print(mul,m,edim);
	newline();
getch();
 
 
	/*for(i=0;i<m;i++)
		for(j=0;j<edim;j++)
			mul[i][j]=t[i][j];
	printf("*******************ttttttttttt\n");
	print(t,m,edim);
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			printf("%lf\t",t[i][j]);	
		}
		printf("\n");
 
	}//keep this saved keep this saved*/
getch();
	printf("\nsigma[][]>>>>>>>>>>>>>>\n");
	computesigma(edim,a);
	print(sigma,edim,edim);
	printf("*******************sigma done\n");
getch();
	sigmainverse(edim);
	print(siginv,edim,edim);
	printf("*********************sigma inverse done\n");
getch();
	ftousqr();
	printf("*********************tou sqr done\n");
getch();
	printf("*********************\n");
	print(sel_evect,n,edim);
	mattranspose(n,edim,sel_evect);
	print(mattrans,edim,n);
	printf("*************************transload done \n");
getch();
/*printf("*******************ttttttttttt\n");
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			fscanf(scores,"%lf\t",&t[i][j]);
			printf("%lf\t",&t[i][j]);
		}
		printf("\n");
		fscanf(scores,"\n");
	}
	*/
	//print(t,m,edim);
getch();
	estimate(mul,transload); // w is t[][],q is transload[][]
	print(scale_esti,m,n);
	printf("esti!@$@%@#%#@%#@%#@\n");
getch();
	fspe();
}
/*********************external functions********************************/
// print matrix
void print(double q[][n],int row,int col)
{
	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
			printf("%lf\t",q[i][j]);
		printf("\n");
 
	}			
}
 
// print 1D array
void print1(double v[],int dimn)
{
	for(i=0;i<dimn;i++)
		printf("%lf\t",v[i]);
}
// mean 
double calmean(double q[][n],int row,int col)
{
	for(i=0;i<col;i++)
	{
		mu11[i]=0.0;
		for(j=0;j<row;j++)
		{
			mu11[i]=mu11[i]+q[j][i];
		}
		mean[i]=mu11[i]/row;
 
	}
	for(i=0;i<col;i++)	
		return(mean[i]);
}
// std
double calstd(double q[][n],int row,int col)
{
	for(j=0;j<col;j++)
	{
		v1[j]=0.0;
		for(i=0;i<row;i++)
		{
			v1[j]=v1[j]+pow((q[i][j]-mean[j]),2);
		}
		var[j]=v1[j]/(row-1);
		std[j]=sqrt(var[j]);
	}
	for(i=0;i<col;i++)
		return(std[i]);	
}
// autoscaling
double autoscale(int row,int col,double q[][n],double r[],double s[])// r is mean, s is std
{  
	for(j=0;j<col;j++)
		for(i=0;i<row;i++)
			mat[i][j]=(q[i][j]-r[j])/s[j];
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
			return(mat[i][j]);
}
 
// covariance
double covariance(int row,int col,double w[][n] )//,double u[]-u[i]-u[j]
{
	for(i=0;i<col;i++)
	{
		for(j=0;j<col;j++)
		{
			cov1[i][j]=0.0;
			for(k=0;k<row;k++) //sum over number of samples
			{
			cov1[i][j]=cov1[i][j]+w[k][i]*w[k][j];
			}
			cov[i][j]=cov1[i][j]/(row-1);
		}
	}
	for(i=0;i<col;i++)
		for(j=0;j<col;j++)
			return(cov[i][j]);		
}
 
// evd 
double evd(double r[][n])
// evd of covariance matrix
{
	//copy covariance matrix to a[][] matrix
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			a[i][j]=cov[i][j];
			//printf("%f\t",a[i][j]);
		}//printf("\dim");
	}
 
	//printing a[][] matrix
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			printf("%lf\t",a[i][j]);	
		}
		printf("\n");
 
	}
	printf("\n");
	////getch();
 
 
	//copy to v 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			v[i][j]=a[i][j];
			//printf("%f\t",a[i][j]);
		}//printf("\dim");
	}
 
myflag=1;
while(myflag==1)
{
	myflag=0;
	count++;
	max=0.0;
	for(i=0;i<dim;i++)
	{
		for(j=i+1;j<dim;j++)
		{
			if((a[i][j]>0)&&(a[i][j]>max))
			{
				max=a[i][j];
				p=i;
				q=j;
			}
			else if((a[i][j]<0) && (a[i][j]<-1*max))
			{
				max=-1*a[i][j];
				p=i;
				q=j;
			}
		}//j loop
	} //i loop
	//printf("p=%d\tq=%d\dim",p,q);
 
	lambda=2*a[p][q];
	mu=a[q][q]-a[p][p];
	w=sqrt((lambda*lambda)+(mu*mu));
 
	if(mu==0)
		sintheta=1/sqrt(2);
	if((mu>0)&&(lambda>0))
		sintheta=sqrt((w-mu)/(2*w));
	if((mu<0)&&(lambda>0))
		sintheta=-sqrt((w+mu)/(2*w));
	if((mu<0)&&(lambda<0))
		sintheta=sqrt((w+mu)/(2*w));
	if((mu>0)&&(lambda<0))
		sintheta=-sqrt((w-mu)/(2*w));
 
	costheta=sqrt(1-(sintheta*sintheta));
	i=0;
	k=0;
	for(i=0;i<dim;i++)
	{
		if((i!=p)&&(i!=q))
		{
			for(k=0;k<dim;k++)
			{
				if((k!=p)&&(k!=q))
				{
					b[i][k]=a[i][k];
				}
				else
				{
					b[i][p]=(a[i][p]*costheta)-(a[i][q]*sintheta);
					b[i][q]=(a[i][p]*sintheta)+(a[i][q]*costheta);
				}
			}
		}
		else
		{
			for(k=0;k<dim;k++)
			{
				if((k!=p)&&(k!=q))
				{
					b[p][k]=(a[p][k]*costheta)-(a[q][k]*sintheta);
					b[q][k]=(a[p][k]*sintheta)+(a[q][k]*costheta);
				}
			}
		}//else
	}//for
 
	b[p][p]=(a[p][p]*costheta*costheta)+(a[q][q]*sintheta*sintheta)-(2*a[p][q]*sintheta*costheta);
	b[q][q]=(a[p][p]*sintheta*sintheta)+(a[q][q]*costheta*costheta)+(2*a[p][q]*sintheta*costheta);
	b[p][q]=0;
	b[q][p]=0;
 
	//creating C matrix
 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			if(j==p)
				c[i][p]=v[i][p]*costheta-v[i][q]*sintheta;
			else if(j==q)
				c[i][q]=v[i][p]*sintheta+v[i][q]*costheta;
			else
				c[i][j]=v[i][j];
		}
	}
 
	//copy c to v for next iteration
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		v[i][j]=c[i][j];
	}
	//copy b to a for next iteration
 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			a[i][j]=b[i][j];
	}
	printf("iter=%d\n",count);
	printf("p=%d\tq=%d\n",p,q);
	printf("maximum element of the upper triangular matrix=%f\n",a[p][q]);
	printf("sintheta=%f\tcostheta=%f\n\n",sintheta,costheta);
 
	printf("A matrix after rotation\n");
 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			printf("%f\t",a[i][j]);
		printf("\n");
	}
	printf("C matrix after rotation\n");
 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			printf("%f\t",c[i][j]);
		printf("\n");
	}
	printf("\n");
 
//checking for looping for next iteration
 
	for(i=0;i<dim;i++)
	{
		for(j=i+1;j<dim;j++)
		{
			if(fabs(b[i][j])!=0.0)
			{
				myflag=1;
				break;
			}
			else if(fabs(b[i][j])==0.0)
				myflag=0;
		}//j loop
		if(myflag==1)
			break;
	}//i loop
	if(myflag==0)
		break;
 
}//while loop
for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			fprintf(val,"%lf\t",a[i][j]);
		fprintf(val,"\n");
	}
for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			fprintf(vec,"%lf\t",c[i][j]);
		fprintf(vec,"\n");
	}
 
 
// eval matrix
 
for(i=0;i<dim;i++)
	for(j=0;j<dim;j++)
		return(a[i][j]);
// e vector matrix
 
for(i=0;i<dim;i++)
	for(j=0;j<dim;j++)
		return(c[i][j]);
}	
 
// combining e val e vect matrices
double combine(int dimn, double u[][n],double v[][n])
{	
	for(i=0;i<dimn+1;i++)
	{
		for(j=0;j<dimn;j++)
		{
			combi[i][j]=v[i][j];
			combi[dimn][j]=u[j][j];
		}	
	}
	printf("\neigen value,eigen vector matrices combined>>>>>>>>>>>\n");
	for(i=0;i<dimn+1;i++)
		for(j=0;j<dimn;j++)
			return(combi[i][j]);
}
// sort()
double sort(int dimn,double w[][n])
{
	for (j=0; j<dimn; j++)
	{
		for (i=0; i<dimn+1; i++)
		{  
			if (combi[dimn][j+1] > combi[dimn][j]) 
			{   
			temp = combi[i][j];
			combi[i][j] = combi[i][j+1];
			combi[i][j+1]= temp;
			}
		}
	}
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			evectsort[i][j]=combi[i][j];
	print(evectsort,n,n);
	printf(">>>>>>\n");
// sorting eigen value matrix,a[][], in descending order of diagonal elements
	for (i=0; i<dimn-1; i++)
	{
		for (j=0; j<dimn-1-i; j++)
		if (a[j+1][j+1] > a[j][j]) 
		{  
		tmp = a[j][j];         
		a[j][j] = a[j+1][j+1];
		a[j+1][j+1] = tmp;
		}
	}
	print(a,n,n);
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(evectsort[i][j]);
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(a[i][j]);
}
// loading
double loading(int dimn,double w[][n],double t[][n])// w is a[][],t is evectsort[][]
{
// calculating totvar
	totvar=0.0;
	for(i=0;i<dimn;i++)
	{
		totvar=totvar+w[i][i];//*a[i][i];
	}
	printf("\ntotvar=%lf\n",totvar);
getch();
 
// selecting loading vectors
	k=0;
	while(pervar<=pertar && k<dim) 
	{		
		pervar1=0.0;
		for(i=0;i<k+1;i++)
		{	
			pervar1=pervar1+w[i][i];//*a[i][i];
			pervar=pervar1/totvar;
		}
		printf("\nk=%d\tpervar=%lf\tpertar=%lf\n",k,pervar,pertar);
		k=k+1;
	} 
 
	edim=k-1;
	printf("edim=%d",edim);
// eigen vectors(from evectsort[][]) from col 0 to col k-1 are chosen as the loading vectors******dim by edim
	for(i=0;i<dimn;i++)
	{
		for(j=0;j<edim;j++)
		{
			sel_evect[i][j]=t[i][j];
		}
	}
	for(i=0;i<dimn;i++)
		for(j=0;j<edim;j++)
			return(sel_evect[i][j]);
}
// fscores
double fscores()
{	
	multiply(m,n,edim,mat,sel_evect);
	print(mul,m,edim);
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			fprintf(scores,"%lf\t",mul[i][j]);
		}
		fprintf(scores,"\n");
	}
	printf("*******************scores to file done\n");
	for(i=0;i<m;i++)
		for(j=0;j<edim;j++)
			return(mul[i][j]);
}
// multiply
double multiply(int g,int h,int r,double w[][n],double l[][n])//multiply two matrices
{
	for(i = 0; i < g; i++) 
		for( j = 0; j < r; j++)
            for( k = 0; k < h; k++) 
                  mul[i][j] +=  w[i][k] * l[k][j];
}
double computesigma(int dimn,double w[][n])// dimn is the edim,w is a[][]
// compute sigma[][]**edim by edim=eigen value matrix corresponding to the selected eigen vectors
{	
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			sigma[i][j]=a[i][j];
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(sigma[i][j]);
}
double sigmainverse(int r) // r is dimension of edim
{	
	matinverse(edim,sigma);
	for(i=0;i<r;i++)
		for(j=0;j<r;j++)
			return(siginv[i][j]);
}
 
double matinverse(int r,double w[][n]) // r is dimension of sigma,w is sigma[][]
{
	for(i=0;i<r;i++)
		siginv[i][i]=1;
 
	for(k=0;k<r-1;k++)
	{
		ix=w[k][k];
		for(i=0;i<r;i++)
		{
			w[k][i]=w[k][i]/ix;
			siginv[k][i]=siginv[k][i]/ix;
		}
		for(i=k+1;i<r;i++)
			im[i-k-1]=w[i][k]/w[k][k];
		for(i=k+1;i<r;i++)
			for(j=0;j<r;j++)
			{
				w[i][j]-=w[k][j]*im[i-k-1];
				siginv[i][j]-=siginv[k][j]*im[i-k-1];
			}
		det*=ix;
	}
	det*=w[k][k];
	if(det==0)
		printf("The inverse doesn't exist as det is zero\n");
	else
	{
		for(j=0;j<r;j++)
			siginv[k][j]/=w[k][k];
		w[k][k]=1;
		for(k=edim-1;k>0;k--)
		{
			for(i=k;i<r;i++)
				im[i-k]=w[k-1][i]/w[i][i];
			for(i=k;i<r;i++)
				for(j=0;j<r;j++)
				{
					w[k-1][j]-=w[i][j]*im[i-k];
					siginv[k-1][j]-=siginv[i][j]*im[i-k];
				}
		}
	for(i=0;i<r;i++)
		for(j=0;j<r;j++)
			return(siginv[i][j]);
	}
}
 
double ftousqr()
// compute tousqr***tousqr=
//(transpose of x)*(loading vectors)*(sigmainv)*(transpose of loading vectors)*(x)
{
	for(h=0;h<m;h++)
	{
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<n;j++)
		  {
			y[i][j]=data[i][j];
 
		  }
 
		}
	// multiply y[1][n]*sel_evect*******1 by edim
		for(i = h; i < h+1; i++) 
           for( j = 0; j < edim; j++)
             for( k = 0; k < n; k++) 
                   tl[i][j] +=  y[i][k] * sel_evect[k][j];
	// multiply tl[][]*siginv[][]
		for(i = h; i < h+1; i++) 
           for( j = 0; j < edim; j++)
             for( k = 0; k < edim; k++) 
                   tls[i][j] +=  tl[i][k] * siginv[k][j];
		// compute transpose of sel_evect
		for(i=0;i<edim;i++)
			for(j=0;j<dim;j++)
				transload[i][j]=sel_evect[j][i];
 
		// multiply tls[][]*transload
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < edim; k++) 
					 tlstl[i][j] += tls[i][k] * transload[k][j];
	// transpose of y[][]****** n by 1
		for(i=0;i<n;i++)
			for(j=h;j<h+1;j++)
				ytrans[i][j]=y[j][i];
	// multiply tlstl[][]*data[][]
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < 1; k++) 
					 tousqr[i][j] += tlstl[i][k] * ytrans[k][j];
 
		printf("tousqr[][]>>>>>>>>>>>>>>>>>\n");
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			printf("%lf\n",tousqr[i][j]);
			fprintf(tou,"%lf\n",tousqr[i][j]);
		  }
		  printf("\n");
		  fprintf(tou,"\n");
		}
		printf("\n");
 
  }	//h loop
}
 
double mattranspose(int p,int k,double w[][n])//specify dimensions of w[][] 
{
	for(i = 0; i < k; i++) 
       for( j = 0; j < p; j++)
		   mattrans[i][j]=w[j][i];
	for(i = 0; i < k; i++) 
       for( j = 0; j < p; j++)
		   return(mattrans[i][j]);
}
double estimate(double w[][2],double r[][dim]) // w is t[][],q is transload[][]
{// scores*transload*******************(mby edim) (edim by dim)
	for(i = 0; i < m; i++) 
       for( j = 0; j < dim; j++)
         for( k = 0; k < edim; k++) 
            scale_esti[i][j] +=  mul[i][k] * transload[k][j];
 
	for(i=0;i<m;i++)
		for(j=0;j<dim;j++)
			return(scale_esti[i][j]);	
}
 
double fspe()
{// compute Scaled-scale_esti
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			zed[i][j]=mat[i][j]-scale_esti[i][j];
// compute spe	
	for(h=0;h<m;h++)
	{
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<n;j++)
		  {
			zedy[i][j]=zed[i][j];
			//printf("%lf\t",zedy[i][j]);
		  }
		  printf("\n");
		}
		printf("\n");
// transpose of y[][]****** n by 1
		for(i=0;i<n;i++)
			for(j=h;j<h+1;j++)
				zedytrans[i][j]=zedy[j][i];
 
// multiply tlstl[][]*data[][]
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < 1; k++) 
					 spe[i][j] += zedy[i][k] * zedytrans[k][j];
 
		printf("spe[][]>>>>>>>>>>>>>>>>>\n");
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			printf("%lf\n",spe[i][j]);
		  }
		  printf("\n");
		}
		printf("\n");
 
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			fprintf(sperror,"%lf\n",spe[i][j]);
		  }
		  printf(tou,"\n");
		}
  }	//h loop	
}
void newline()
{
	printf("\n\n");
}
#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>

#define m 10   // m samples 
#define n 5   // n variables
#define dim 5 //whatever var dim is ***dim of covariance matrix
#define pertar 0.99

double scale_esti[m][n],mul[m][n];
double lambda,max,mu,w,sintheta,costheta,totvar,pervar,pervar1,temp,tmp;
double a[dim][dim],b[dim][dim],v[dim][dim],c[dim][dim];
double ix,im[10]={0},siginv[dim][dim]={0},det=1;// do not disturb siginv
double tousqr[1][1],zedy[1][5],evectsort[n][n];
double tl[m][2],tls[m][2],transload[2][5],tlstl[m][n];
double y[1][n],ytrans[n][1],v1[n],var[n],std[n];
double cov[n][n],cov1[n][n],sigma[dim][dim],t[m][n],mat[m][n];// do not disturb sigma
double mean[n],mu11[n],sel_evect[n][5];// do not disturb sel_evect
double combi[n+1][n],zed[m][n],zedytrans[n][1],spe[m][n];// spe should be 1,1;
double mattrans[2][5];
int i,h,j,z,edim,k,p,q,myflag=0,count=0;

double data[10][5]={{92,99,1,8,15},{98,80,7,14,16},{4,81,88,20,22},{85,87,19,21,3},{86,93,25,2,9},{17,24,76,83,90},{23,5,82,89,91},{79,6,13,95,97},{10,	12,	94,	96,	78},{11,	18,	100,	77,	84}};

FILE *val,*vec,*tou,*in,*scores,*sperror;
// functions
void   print(double q[][n],int,int);
void   newline();
void   print1(double v[],int);
double calmean(double q[][n],int ,int);
double calstd(double q[][n],int,int);
double autoscale(int, int, double q[][n],double r[],double s[]);
double covariance(int row,int col,double w[][n] );//,double u[n]
double evd(double r[][n]);
double combine(int dimn, double u[][n],double v[][n]);
double sort(int dimn,double w[][n]);
double loading(int dimn,double w[][n],double t[][n]);// w is a[][],z is evectsort[][]
double fscores();
double multiply(int g,int h,int r,double w[][n],double l[][n]);//multiply two matrices
double matinverse(int r,double w[][n]);
double computesigma(int dimn,double w[][n]);
double sigmainverse(int);
double ftousqr();
double mattranspose(int p,int k,double w[][n]);
double estimate(double w[][2],double r[][dim]);
double fspe();

main()
{
	val=fopen("eigen_val3.dat","w");
	vec=fopen("eigen_vect3.dat","w");
	tou=fopen("tousqr3.xls","w");
	scores=fopen("scores3.xls","w");
	sperror=fopen("spe.xls","w");
	scores=fopen("scores3.xls","r");
	
	print(data,m,n);
getch();
	calmean(data,m,n);
	print1(mean,n);
	newline();
getch();
	calstd(data,m,n);
	print1(std,n);
	newline();
getch();
	autoscale(m,n,data,mean,std);
	print(mat,m,n);
	newline();
getch();
	calmean(mat,m,n);
	print1(mean,n);
	newline();
getch();
	covariance(10,5,mat);
	print(cov,n,n);
	newline();
getch();
	printf("evd of covariance\n");
	evd(cov);
getch();
	printf("eigen value matrix\n");
	print(a,n,n);
	newline();
getch();
	printf("eigen vector matrix\n");
	print(c,n,n);
	newline();
getch();
	printf("\ncombined matrix >>>>>>>>>>>\n");
	combine(n,a,c);
	print(combi,n+1,n);
	newline();
getch();
	printf("\ncombined matrix sorted>>>>>>>>>>>\n");
	sort(n,combi);
	newline();
	printf("*********************sort done\n");
getch();
	loading(n,a,evectsort);
	newline();
	print(sel_evect,dim,edim);
	newline();
	printf("*******************selection done\n");
getch();
	fscores();
	printf("*******************scores done\n");
	print(mul,m,edim);
	newline();
getch();
	
	
	/*for(i=0;i<m;i++)
		for(j=0;j<edim;j++)
			mul[i][j]=t[i][j];
	printf("*******************ttttttttttt\n");
	print(t,m,edim);
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			printf("%lf\t",t[i][j]);	
		}
		printf("\n");
		
	}//keep this saved keep this saved*/
getch();
	printf("\nsigma[][]>>>>>>>>>>>>>>\n");
	computesigma(edim,a);
	print(sigma,edim,edim);
	printf("*******************sigma done\n");
getch();
	sigmainverse(edim);
	print(siginv,edim,edim);
	printf("*********************sigma inverse done\n");
getch();
	ftousqr();
	printf("*********************tou sqr done\n");
getch();
	printf("*********************\n");
	print(sel_evect,n,edim);
	mattranspose(n,edim,sel_evect);
	print(mattrans,edim,n);
	printf("*************************transload done \n");
getch();
/*printf("*******************ttttttttttt\n");
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			fscanf(scores,"%lf\t",&t[i][j]);
			printf("%lf\t",&t[i][j]);
		}
		printf("\n");
		fscanf(scores,"\n");
	}
	*/
	//print(t,m,edim);
getch();
	estimate(mul,transload); // w is t[][],q is transload[][]
	print(scale_esti,m,n);
	printf("esti!@$@%@#%#@%#@%#@\n");
getch();
	fspe();
}
/*********************external functions********************************/
// print matrix
void print(double q[][n],int row,int col)
{
	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
			printf("%lf\t",q[i][j]);
		printf("\n");
		
	}			
}

// print 1D array
void print1(double v[],int dimn)
{
	for(i=0;i<dimn;i++)
		printf("%lf\t",v[i]);
}
// mean 
double calmean(double q[][n],int row,int col)
{
	for(i=0;i<col;i++)
	{
		mu11[i]=0.0;
		for(j=0;j<row;j++)
		{
			mu11[i]=mu11[i]+q[j][i];
		}
		mean[i]=mu11[i]/row;
				
	}
	for(i=0;i<col;i++)	
		return(mean[i]);
}
// std
double calstd(double q[][n],int row,int col)
{
	for(j=0;j<col;j++)
	{
		v1[j]=0.0;
		for(i=0;i<row;i++)
		{
			v1[j]=v1[j]+pow((q[i][j]-mean[j]),2);
		}
		var[j]=v1[j]/(row-1);
		std[j]=sqrt(var[j]);
	}
	for(i=0;i<col;i++)
		return(std[i]);	
}
// autoscaling
double autoscale(int row,int col,double q[][n],double r[],double s[])// r is mean, s is std
{  
	for(j=0;j<col;j++)
		for(i=0;i<row;i++)
			mat[i][j]=(q[i][j]-r[j])/s[j];
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
			return(mat[i][j]);
}

// covariance
double covariance(int row,int col,double w[][n] )//,double u[]-u[i]-u[j]
{
	for(i=0;i<col;i++)
	{
		for(j=0;j<col;j++)
		{
			cov1[i][j]=0.0;
			for(k=0;k<row;k++) //sum over number of samples
			{
			cov1[i][j]=cov1[i][j]+w[k][i]*w[k][j];
			}
			cov[i][j]=cov1[i][j]/(row-1);
		}
	}
	for(i=0;i<col;i++)
		for(j=0;j<col;j++)
			return(cov[i][j]);		
}

// evd 
double evd(double r[][n])
// evd of covariance matrix
{
	//copy covariance matrix to a[][] matrix
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			a[i][j]=cov[i][j];
			//printf("%f\t",a[i][j]);
		}//printf("\dim");
	}

	//printing a[][] matrix
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			printf("%lf\t",a[i][j]);	
		}
		printf("\n");
		
	}
	printf("\n");
	////getch();


	//copy to v 
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			v[i][j]=a[i][j];
			//printf("%f\t",a[i][j]);
		}//printf("\dim");
	}

myflag=1;
while(myflag==1)
{
	myflag=0;
	count++;
	max=0.0;
	for(i=0;i<dim;i++)
	{
		for(j=i+1;j<dim;j++)
		{
			if((a[i][j]>0)&&(a[i][j]>max))
			{
				max=a[i][j];
				p=i;
				q=j;
			}
			else if((a[i][j]<0) && (a[i][j]<-1*max))
			{
				max=-1*a[i][j];
				p=i;
				q=j;
			}
		}//j loop
	} //i loop
	//printf("p=%d\tq=%d\dim",p,q);

	lambda=2*a[p][q];
	mu=a[q][q]-a[p][p];
	w=sqrt((lambda*lambda)+(mu*mu));

	if(mu==0)
		sintheta=1/sqrt(2);
	if((mu>0)&&(lambda>0))
		sintheta=sqrt((w-mu)/(2*w));
	if((mu<0)&&(lambda>0))
		sintheta=-sqrt((w+mu)/(2*w));
	if((mu<0)&&(lambda<0))
		sintheta=sqrt((w+mu)/(2*w));
	if((mu>0)&&(lambda<0))
		sintheta=-sqrt((w-mu)/(2*w));

	costheta=sqrt(1-(sintheta*sintheta));
	i=0;
	k=0;
	for(i=0;i<dim;i++)
	{
		if((i!=p)&&(i!=q))
		{
			for(k=0;k<dim;k++)
			{
				if((k!=p)&&(k!=q))
				{
					b[i][k]=a[i][k];
				}
				else
				{
					b[i][p]=(a[i][p]*costheta)-(a[i][q]*sintheta);
					b[i][q]=(a[i][p]*sintheta)+(a[i][q]*costheta);
				}
			}
		}
		else
		{
			for(k=0;k<dim;k++)
			{
				if((k!=p)&&(k!=q))
				{
					b[p][k]=(a[p][k]*costheta)-(a[q][k]*sintheta);
					b[q][k]=(a[p][k]*sintheta)+(a[q][k]*costheta);
				}
			}
		}//else
	}//for

	b[p][p]=(a[p][p]*costheta*costheta)+(a[q][q]*sintheta*sintheta)-(2*a[p][q]*sintheta*costheta);
	b[q][q]=(a[p][p]*sintheta*sintheta)+(a[q][q]*costheta*costheta)+(2*a[p][q]*sintheta*costheta);
	b[p][q]=0;
	b[q][p]=0;

	//creating C matrix

	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		{
			if(j==p)
				c[i][p]=v[i][p]*costheta-v[i][q]*sintheta;
			else if(j==q)
				c[i][q]=v[i][p]*sintheta+v[i][q]*costheta;
			else
				c[i][j]=v[i][j];
		}
	}

	//copy c to v for next iteration
	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
		v[i][j]=c[i][j];
	}
	//copy b to a for next iteration

	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			a[i][j]=b[i][j];
	}
	printf("iter=%d\n",count);
	printf("p=%d\tq=%d\n",p,q);
	printf("maximum element of the upper triangular matrix=%f\n",a[p][q]);
	printf("sintheta=%f\tcostheta=%f\n\n",sintheta,costheta);
	
	printf("A matrix after rotation\n");

	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			printf("%f\t",a[i][j]);
		printf("\n");
	}
	printf("C matrix after rotation\n");

	for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			printf("%f\t",c[i][j]);
		printf("\n");
	}
	printf("\n");

//checking for looping for next iteration

	for(i=0;i<dim;i++)
	{
		for(j=i+1;j<dim;j++)
		{
			if(fabs(b[i][j])!=0.0)
			{
				myflag=1;
				break;
			}
			else if(fabs(b[i][j])==0.0)
				myflag=0;
		}//j loop
		if(myflag==1)
			break;
	}//i loop
	if(myflag==0)
		break;

}//while loop
for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			fprintf(val,"%lf\t",a[i][j]);
		fprintf(val,"\n");
	}
for(i=0;i<dim;i++)
	{
		for(j=0;j<dim;j++)
			fprintf(vec,"%lf\t",c[i][j]);
		fprintf(vec,"\n");
	}


// eval matrix

for(i=0;i<dim;i++)
	for(j=0;j<dim;j++)
		return(a[i][j]);
// e vector matrix

for(i=0;i<dim;i++)
	for(j=0;j<dim;j++)
		return(c[i][j]);
}	

// combining e val e vect matrices
double combine(int dimn, double u[][n],double v[][n])
{	
	for(i=0;i<dimn+1;i++)
	{
		for(j=0;j<dimn;j++)
		{
			combi[i][j]=v[i][j];
			combi[dimn][j]=u[j][j];
		}	
	}
	printf("\neigen value,eigen vector matrices combined>>>>>>>>>>>\n");
	for(i=0;i<dimn+1;i++)
		for(j=0;j<dimn;j++)
			return(combi[i][j]);
}
// sort()
double sort(int dimn,double w[][n])
{
	for (j=0; j<dimn; j++)
	{
		for (i=0; i<dimn+1; i++)
		{  
			if (combi[dimn][j+1] > combi[dimn][j]) 
			{   
			temp = combi[i][j];
			combi[i][j] = combi[i][j+1];
			combi[i][j+1]= temp;
			}
		}
	}
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			evectsort[i][j]=combi[i][j];
	print(evectsort,n,n);
	printf(">>>>>>\n");
// sorting eigen value matrix,a[][], in descending order of diagonal elements
	for (i=0; i<dimn-1; i++)
	{
		for (j=0; j<dimn-1-i; j++)
		if (a[j+1][j+1] > a[j][j]) 
		{  
		tmp = a[j][j];         
		a[j][j] = a[j+1][j+1];
		a[j+1][j+1] = tmp;
		}
	}
	print(a,n,n);
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(evectsort[i][j]);
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(a[i][j]);
}
// loading
double loading(int dimn,double w[][n],double t[][n])// w is a[][],t is evectsort[][]
{
// calculating totvar
	totvar=0.0;
	for(i=0;i<dimn;i++)
	{
		totvar=totvar+w[i][i];//*a[i][i];
	}
	printf("\ntotvar=%lf\n",totvar);
getch();

// selecting loading vectors
	k=0;
	while(pervar<=pertar && k<dim) 
	{		
		pervar1=0.0;
		for(i=0;i<k+1;i++)
		{	
			pervar1=pervar1+w[i][i];//*a[i][i];
			pervar=pervar1/totvar;
		}
		printf("\nk=%d\tpervar=%lf\tpertar=%lf\n",k,pervar,pertar);
		k=k+1;
	} 

	edim=k-1;
	printf("edim=%d",edim);
// eigen vectors(from evectsort[][]) from col 0 to col k-1 are chosen as the loading vectors******dim by edim
	for(i=0;i<dimn;i++)
	{
		for(j=0;j<edim;j++)
		{
			sel_evect[i][j]=t[i][j];
		}
	}
	for(i=0;i<dimn;i++)
		for(j=0;j<edim;j++)
			return(sel_evect[i][j]);
}
// fscores
double fscores()
{	
	multiply(m,n,edim,mat,sel_evect);
	print(mul,m,edim);
	for(i=0;i<m;i++)
	{
		for(j=0;j<edim;j++)
		{
			fprintf(scores,"%lf\t",mul[i][j]);
		}
		fprintf(scores,"\n");
	}
	printf("*******************scores to file done\n");
	for(i=0;i<m;i++)
		for(j=0;j<edim;j++)
			return(mul[i][j]);
}
// multiply
double multiply(int g,int h,int r,double w[][n],double l[][n])//multiply two matrices
{
	for(i = 0; i < g; i++) 
		for( j = 0; j < r; j++)
            for( k = 0; k < h; k++) 
                  mul[i][j] +=  w[i][k] * l[k][j];
}
double computesigma(int dimn,double w[][n])// dimn is the edim,w is a[][]
// compute sigma[][]**edim by edim=eigen value matrix corresponding to the selected eigen vectors
{	
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			sigma[i][j]=a[i][j];
	for(i=0;i<dimn;i++)
		for(j=0;j<dimn;j++)
			return(sigma[i][j]);
}
double sigmainverse(int r) // r is dimension of edim
{	
	matinverse(edim,sigma);
	for(i=0;i<r;i++)
		for(j=0;j<r;j++)
			return(siginv[i][j]);
}

double matinverse(int r,double w[][n]) // r is dimension of sigma,w is sigma[][]
{
	for(i=0;i<r;i++)
		siginv[i][i]=1;

	for(k=0;k<r-1;k++)
	{
		ix=w[k][k];
		for(i=0;i<r;i++)
		{
			w[k][i]=w[k][i]/ix;
			siginv[k][i]=siginv[k][i]/ix;
		}
		for(i=k+1;i<r;i++)
			im[i-k-1]=w[i][k]/w[k][k];
		for(i=k+1;i<r;i++)
			for(j=0;j<r;j++)
			{
				w[i][j]-=w[k][j]*im[i-k-1];
				siginv[i][j]-=siginv[k][j]*im[i-k-1];
			}
		det*=ix;
	}
	det*=w[k][k];
	if(det==0)
		printf("The inverse doesn't exist as det is zero\n");
	else
	{
		for(j=0;j<r;j++)
			siginv[k][j]/=w[k][k];
		w[k][k]=1;
		for(k=edim-1;k>0;k--)
		{
			for(i=k;i<r;i++)
				im[i-k]=w[k-1][i]/w[i][i];
			for(i=k;i<r;i++)
				for(j=0;j<r;j++)
				{
					w[k-1][j]-=w[i][j]*im[i-k];
					siginv[k-1][j]-=siginv[i][j]*im[i-k];
				}
		}
	for(i=0;i<r;i++)
		for(j=0;j<r;j++)
			return(siginv[i][j]);
	}
}

double ftousqr()
// compute tousqr***tousqr=
//(transpose of x)*(loading vectors)*(sigmainv)*(transpose of loading vectors)*(x)
{
	for(h=0;h<m;h++)
	{
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<n;j++)
		  {
			y[i][j]=data[i][j];
			
		  }
		  
		}
	// multiply y[1][n]*sel_evect*******1 by edim
		for(i = h; i < h+1; i++) 
           for( j = 0; j < edim; j++)
             for( k = 0; k < n; k++) 
                   tl[i][j] +=  y[i][k] * sel_evect[k][j];
	// multiply tl[][]*siginv[][]
		for(i = h; i < h+1; i++) 
           for( j = 0; j < edim; j++)
             for( k = 0; k < edim; k++) 
                   tls[i][j] +=  tl[i][k] * siginv[k][j];
		// compute transpose of sel_evect
		for(i=0;i<edim;i++)
			for(j=0;j<dim;j++)
				transload[i][j]=sel_evect[j][i];
		
		// multiply tls[][]*transload
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < edim; k++) 
					 tlstl[i][j] += tls[i][k] * transload[k][j];
	// transpose of y[][]****** n by 1
		for(i=0;i<n;i++)
			for(j=h;j<h+1;j++)
				ytrans[i][j]=y[j][i];
	// multiply tlstl[][]*data[][]
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < 1; k++) 
					 tousqr[i][j] += tlstl[i][k] * ytrans[k][j];
		
		printf("tousqr[][]>>>>>>>>>>>>>>>>>\n");
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			printf("%lf\n",tousqr[i][j]);
			fprintf(tou,"%lf\n",tousqr[i][j]);
		  }
		  printf("\n");
		  fprintf(tou,"\n");
		}
		printf("\n");

  }	//h loop
}

double mattranspose(int p,int k,double w[][n])//specify dimensions of w[][] 
{
	for(i = 0; i < k; i++) 
       for( j = 0; j < p; j++)
		   mattrans[i][j]=w[j][i];
	for(i = 0; i < k; i++) 
       for( j = 0; j < p; j++)
		   return(mattrans[i][j]);
}
double estimate(double w[][2],double r[][dim]) // w is t[][],q is transload[][]
{// scores*transload*******************(mby edim) (edim by dim)
	for(i = 0; i < m; i++) 
       for( j = 0; j < dim; j++)
         for( k = 0; k < edim; k++) 
            scale_esti[i][j] +=  mul[i][k] * transload[k][j];
	
	for(i=0;i<m;i++)
		for(j=0;j<dim;j++)
			return(scale_esti[i][j]);	
}

double fspe()
{// compute Scaled-scale_esti
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			zed[i][j]=mat[i][j]-scale_esti[i][j];
// compute spe	
	for(h=0;h<m;h++)
	{
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<n;j++)
		  {
			zedy[i][j]=zed[i][j];
			//printf("%lf\t",zedy[i][j]);
		  }
		  printf("\n");
		}
		printf("\n");
// transpose of y[][]****** n by 1
		for(i=0;i<n;i++)
			for(j=h;j<h+1;j++)
				zedytrans[i][j]=zedy[j][i];
		
// multiply tlstl[][]*data[][]
		for(i = h; i < h+1; i++) 
			for( j = 0; j < dim; j++)
				for( k = 0; k < 1; k++) 
					 spe[i][j] += zedy[i][k] * zedytrans[k][j];
		
		printf("spe[][]>>>>>>>>>>>>>>>>>\n");
		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			printf("%lf\n",spe[i][j]);
		  }
		  printf("\n");
		}
		printf("\n");

		for(i=h;i<h+1;i++)
		{
		  for(j=0;j<1;j++)
		  {
			fprintf(sperror,"%lf\n",spe[i][j]);
		  }
		  printf(tou,"\n");
		}
  }	//h loop	
}
void newline()
{
	printf("\n\n");
}
Similar Threads
I've been highjacked by Ownbox.com in Viruses, Spyware and other Nasties
COM Memory Layout, Vtable, and Function Pointer Issues And Questions in C++
can someone review in Viruses, Spyware and other Nasties
Ghost 2003 & USB restore in Windows Software
wireless router not transmitting in Networking
sv_1124
Reputation Points: 10
Solved Threads: 0
Newbie Poster
Offline
1 posts
since Sep 2009
Permalink
Sep 20th, 2009
Re: pca in c

Are you kidding me? Why are you posting this as a code snippet? You should have posted this in the C forum as a "normal" thread". And quite frankly unless you are willing to put some effort in and pinpoint the minimal code that demonstrates your problem, I'd guess you're going to have a hard time finding people here prepared to help you when you dump the better part of 800 lines of code and virtually say "please look and fix for me".
yellowSnow
Reputation Points: 651
Solved Threads: 35
Posting Whiz in Training
Offline
201 posts
since Jul 2009
Permalink
Sep 20th, 2009
Re: pca in c

In addition (personally one of my favourite links):
http://cboard.cprogramming.com/c-pro...t-process.html
tux4life
Reputation Points: 2125
Solved Threads: 243
Postaholic
Offline
2,105 posts
since Feb 2009
This thread is more than three months old

No one has posted to this discussion for at least three months. Please let old threads die and do not reply to them unless you feel you have something new and valuable to contribute that absolutely must be added to make the discussion complete. Otherwise, please start a new thread in this forum instead.
Previous Thread in C Forum Timeline: Can anyone help me..
Next Thread in C Forum Timeline: can anyone help me with writing output
DaniWeb Message

Cancel Changes
User Name
Password