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
#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");
}
