C RSS C RSS

pca in c

Reply

•
•
Join Date: Sep 2009
Posts: 1
Reputation: sv_1124 is an unknown quantity at this point 
Solved Threads: 0
sv_1124 sv_1124 is offline Offline
Newbie Poster

pca in c

 
0
  #1
Sep 20th, 2009
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) 
  1. #include<iostream.h>
  2. #include<stdio.h>
  3. #include<stdlib.h>
  4. #include<math.h>
  5. #include<conio.h>
  6.  
  7. #define m 10   // m samples 
  8. #define n 5   // n variables
  9. #define dim 5 //whatever var dim is ***dim of covariance matrix
  10. #define pertar 0.99
  11.  
  12. double scale_esti[m][n],mul[m][n];
  13. double lambda,max,mu,w,sintheta,costheta,totvar,pervar,pervar1,temp,tmp;
  14. double a[dim][dim],b[dim][dim],v[dim][dim],c[dim][dim];
  15. double ix,im[10]={0},siginv[dim][dim]={0},det=1;// do not disturb siginv
  16. double tousqr[1][1],zedy[1][5],evectsort[n][n];
  17. double tl[m][2],tls[m][2],transload[2][5],tlstl[m][n];
  18. double y[1][n],ytrans[n][1],v1[n],var[n],std[n];
  19. double cov[n][n],cov1[n][n],sigma[dim][dim],t[m][n],mat[m][n];// do not disturb sigma
  20. double mean[n],mu11[n],sel_evect[n][5];// do not disturb sel_evect
  21. double combi[n+1][n],zed[m][n],zedytrans[n][1],spe[m][n];// spe should be 1,1;
  22. double mattrans[2][5];
  23. int i,h,j,z,edim,k,p,q,myflag=0,count=0;
  24.  
  25. 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}};
  26.  
  27. FILE *val,*vec,*tou,*in,*scores,*sperror;
  28. // functions
  29. void   print(double q[][n],int,int);
  30. void   newline();
  31. void   print1(double v[],int);
  32. double calmean(double q[][n],int ,int);
  33. double calstd(double q[][n],int,int);
  34. double autoscale(int, int, double q[][n],double r[],double s[]);
  35. double covariance(int row,int col,double w[][n] );//,double u[n]
  36. double evd(double r[][n]);
  37. double combine(int dimn, double u[][n],double v[][n]);
  38. double sort(int dimn,double w[][n]);
  39. double loading(int dimn,double w[][n],double t[][n]);// w is a[][],z is evectsort[][]
  40. double fscores();
  41. double multiply(int g,int h,int r,double w[][n],double l[][n]);//multiply two matrices
  42. double matinverse(int r,double w[][n]);
  43. double computesigma(int dimn,double w[][n]);
  44. double sigmainverse(int);
  45. double ftousqr();
  46. double mattranspose(int p,int k,double w[][n]);
  47. double estimate(double w[][2],double r[][dim]);
  48. double fspe();
  49.  
  50. main()
  51. {
  52. 	val=fopen("eigen_val3.dat","w");
  53. 	vec=fopen("eigen_vect3.dat","w");
  54. 	tou=fopen("tousqr3.xls","w");
  55. 	scores=fopen("scores3.xls","w");
  56. 	sperror=fopen("spe.xls","w");
  57. 	scores=fopen("scores3.xls","r");
  58.  
  59. 	print(data,m,n);
  60. getch();
  61. 	calmean(data,m,n);
  62. 	print1(mean,n);
  63. 	newline();
  64. getch();
  65. 	calstd(data,m,n);
  66. 	print1(std,n);
  67. 	newline();
  68. getch();
  69. 	autoscale(m,n,data,mean,std);
  70. 	print(mat,m,n);
  71. 	newline();
  72. getch();
  73. 	calmean(mat,m,n);
  74. 	print1(mean,n);
  75. 	newline();
  76. getch();
  77. 	covariance(10,5,mat);
  78. 	print(cov,n,n);
  79. 	newline();
  80. getch();
  81. 	printf("evd of covariance\n");
  82. 	evd(cov);
  83. getch();
  84. 	printf("eigen value matrix\n");
  85. 	print(a,n,n);
  86. 	newline();
  87. getch();
  88. 	printf("eigen vector matrix\n");
  89. 	print(c,n,n);
  90. 	newline();
  91. getch();
  92. 	printf("\ncombined matrix >>>>>>>>>>>\n");
  93. 	combine(n,a,c);
  94. 	print(combi,n+1,n);
  95. 	newline();
  96. getch();
  97. 	printf("\ncombined matrix sorted>>>>>>>>>>>\n");
  98. 	sort(n,combi);
  99. 	newline();
  100. 	printf("*********************sort done\n");
  101. getch();
  102. 	loading(n,a,evectsort);
  103. 	newline();
  104. 	print(sel_evect,dim,edim);
  105. 	newline();
  106. 	printf("*******************selection done\n");
  107. getch();
  108. 	fscores();
  109. 	printf("*******************scores done\n");
  110. 	print(mul,m,edim);
  111. 	newline();
  112. getch();
  113.  
  114.  
  115. 	/*for(i=0;i<m;i++)
  116. 		for(j=0;j<edim;j++)
  117. 			mul[i][j]=t[i][j];
  118. 	printf("*******************ttttttttttt\n");
  119. 	print(t,m,edim);
  120. 	for(i=0;i<m;i++)
  121. 	{
  122. 		for(j=0;j<edim;j++)
  123. 		{
  124. 			printf("%lf\t",t[i][j]);	
  125. 		}
  126. 		printf("\n");
  127.  
  128. 	}//keep this saved keep this saved*/
  129. getch();
  130. 	printf("\nsigma[][]>>>>>>>>>>>>>>\n");
  131. 	computesigma(edim,a);
  132. 	print(sigma,edim,edim);
  133. 	printf("*******************sigma done\n");
  134. getch();
  135. 	sigmainverse(edim);
  136. 	print(siginv,edim,edim);
  137. 	printf("*********************sigma inverse done\n");
  138. getch();
  139. 	ftousqr();
  140. 	printf("*********************tou sqr done\n");
  141. getch();
  142. 	printf("*********************\n");
  143. 	print(sel_evect,n,edim);
  144. 	mattranspose(n,edim,sel_evect);
  145. 	print(mattrans,edim,n);
  146. 	printf("*************************transload done \n");
  147. getch();
  148. /*printf("*******************ttttttttttt\n");
  149. 	for(i=0;i<m;i++)
  150. 	{
  151. 		for(j=0;j<edim;j++)
  152. 		{
  153. 			fscanf(scores,"%lf\t",&t[i][j]);
  154. 			printf("%lf\t",&t[i][j]);
  155. 		}
  156. 		printf("\n");
  157. 		fscanf(scores,"\n");
  158. 	}
  159. 	*/
  160. 	//print(t,m,edim);
  161. getch();
  162. 	estimate(mul,transload); // w is t[][],q is transload[][]
  163. 	print(scale_esti,m,n);
  164. 	printf("esti!@$@%@#%#@%#@%#@\n");
  165. getch();
  166. 	fspe();
  167. }
  168. /*********************external functions********************************/
  169. // print matrix
  170. void print(double q[][n],int row,int col)
  171. {
  172. 	for(i=0;i<row;i++)
  173. 	{
  174. 		for(j=0;j<col;j++)
  175. 			printf("%lf\t",q[i][j]);
  176. 		printf("\n");
  177.  
  178. 	}			
  179. }
  180.  
  181. // print 1D array
  182. void print1(double v[],int dimn)
  183. {
  184. 	for(i=0;i<dimn;i++)
  185. 		printf("%lf\t",v[i]);
  186. }
  187. // mean 
  188. double calmean(double q[][n],int row,int col)
  189. {
  190. 	for(i=0;i<col;i++)
  191. 	{
  192. 		mu11[i]=0.0;
  193. 		for(j=0;j<row;j++)
  194. 		{
  195. 			mu11[i]=mu11[i]+q[j][i];
  196. 		}
  197. 		mean[i]=mu11[i]/row;
  198.  
  199. 	}
  200. 	for(i=0;i<col;i++)	
  201. 		return(mean[i]);
  202. }
  203. // std
  204. double calstd(double q[][n],int row,int col)
  205. {
  206. 	for(j=0;j<col;j++)
  207. 	{
  208. 		v1[j]=0.0;
  209. 		for(i=0;i<row;i++)
  210. 		{
  211. 			v1[j]=v1[j]+pow((q[i][j]-mean[j]),2);
  212. 		}
  213. 		var[j]=v1[j]/(row-1);
  214. 		std[j]=sqrt(var[j]);
  215. 	}
  216. 	for(i=0;i<col;i++)
  217. 		return(std[i]);	
  218. }
  219. // autoscaling
  220. double autoscale(int row,int col,double q[][n],double r[],double s[])// r is mean, s is std
  221. {  
  222. 	for(j=0;j<col;j++)
  223. 		for(i=0;i<row;i++)
  224. 			mat[i][j]=(q[i][j]-r[j])/s[j];
  225. 	for(i=0;i<row;i++)
  226. 		for(j=0;j<col;j++)
  227. 			return(mat[i][j]);
  228. }
  229.  
  230. // covariance
  231. double covariance(int row,int col,double w[][n] )//,double u[]-u[i]-u[j]
  232. {
  233. 	for(i=0;i<col;i++)
  234. 	{
  235. 		for(j=0;j<col;j++)
  236. 		{
  237. 			cov1[i][j]=0.0;
  238. 			for(k=0;k<row;k++) //sum over number of samples
  239. 			{
  240. 			cov1[i][j]=cov1[i][j]+w[k][i]*w[k][j];
  241. 			}
  242. 			cov[i][j]=cov1[i][j]/(row-1);
  243. 		}
  244. 	}
  245. 	for(i=0;i<col;i++)
  246. 		for(j=0;j<col;j++)
  247. 			return(cov[i][j]);		
  248. }
  249.  
  250. // evd 
  251. double evd(double r[][n])
  252. // evd of covariance matrix
  253. {
  254. 	//copy covariance matrix to a[][] matrix
  255. 	for(i=0;i<dim;i++)
  256. 	{
  257. 		for(j=0;j<dim;j++)
  258. 		{
  259. 			a[i][j]=cov[i][j];
  260. 			//printf("%f\t",a[i][j]);
  261. 		}//printf("\dim");
  262. 	}
  263.  
  264. 	//printing a[][] matrix
  265. 	for(i=0;i<n;i++)
  266. 	{
  267. 		for(j=0;j<n;j++)
  268. 		{
  269. 			printf("%lf\t",a[i][j]);	
  270. 		}
  271. 		printf("\n");
  272.  
  273. 	}
  274. 	printf("\n");
  275. 	////getch();
  276.  
  277.  
  278. 	//copy to v 
  279. 	for(i=0;i<dim;i++)
  280. 	{
  281. 		for(j=0;j<dim;j++)
  282. 		{
  283. 			v[i][j]=a[i][j];
  284. 			//printf("%f\t",a[i][j]);
  285. 		}//printf("\dim");
  286. 	}
  287.  
  288. myflag=1;
  289. while(myflag==1)
  290. {
  291. 	myflag=0;
  292. 	count++;
  293. 	max=0.0;
  294. 	for(i=0;i<dim;i++)
  295. 	{
  296. 		for(j=i+1;j<dim;j++)
  297. 		{
  298. 			if((a[i][j]>0)&&(a[i][j]>max))
  299. 			{
  300. 				max=a[i][j];
  301. 				p=i;
  302. 				q=j;
  303. 			}
  304. 			else if((a[i][j]<0) && (a[i][j]<-1*max))
  305. 			{
  306. 				max=-1*a[i][j];
  307. 				p=i;
  308. 				q=j;
  309. 			}
  310. 		}//j loop
  311. 	} //i loop
  312. 	//printf("p=%d\tq=%d\dim",p,q);
  313.  
  314. 	lambda=2*a[p][q];
  315. 	mu=a[q][q]-a[p][p];
  316. 	w=sqrt((lambda*lambda)+(mu*mu));
  317.  
  318. 	if(mu==0)
  319. 		sintheta=1/sqrt(2);
  320. 	if((mu>0)&&(lambda>0))
  321. 		sintheta=sqrt((w-mu)/(2*w));
  322. 	if((mu<0)&&(lambda>0))
  323. 		sintheta=-sqrt((w+mu)/(2*w));
  324. 	if((mu<0)&&(lambda<0))
  325. 		sintheta=sqrt((w+mu)/(2*w));
  326. 	if((mu>0)&&(lambda<0))
  327. 		sintheta=-sqrt((w-mu)/(2*w));
  328.  
  329. 	costheta=sqrt(1-(sintheta*sintheta));
  330. 	i=0;
  331. 	k=0;
  332. 	for(i=0;i<dim;i++)
  333. 	{
  334. 		if((i!=p)&&(i!=q))
  335. 		{
  336. 			for(k=0;k<dim;k++)
  337. 			{
  338. 				if((k!=p)&&(k!=q))
  339. 				{
  340. 					b[i][k]=a[i][k];
  341. 				}
  342. 				else
  343. 				{
  344. 					b[i][p]=(a[i][p]*costheta)-(a[i][q]*sintheta);
  345. 					b[i][q]=(a[i][p]*sintheta)+(a[i][q]*costheta);
  346. 				}
  347. 			}
  348. 		}
  349. 		else
  350. 		{
  351. 			for(k=0;k<dim;k++)
  352. 			{
  353. 				if((k!=p)&&(k!=q))
  354. 				{
  355. 					b[p][k]=(a[p][k]*costheta)-(a[q][k]*sintheta);
  356. 					b[q][k]=(a[p][k]*sintheta)+(a[q][k]*costheta);
  357. 				}
  358. 			}
  359. 		}//else
  360. 	}//for
  361.  
  362. 	b[p][p]=(a[p][p]*costheta*costheta)+(a[q][q]*sintheta*sintheta)-(2*a[p][q]*sintheta*costheta);
  363. 	b[q][q]=(a[p][p]*sintheta*sintheta)+(a[q][q]*costheta*costheta)+(2*a[p][q]*sintheta*costheta);
  364. 	b[p][q]=0;
  365. 	b[q][p]=0;
  366.  
  367. 	//creating C matrix
  368.  
  369. 	for(i=0;i<dim;i++)
  370. 	{
  371. 		for(j=0;j<dim;j++)
  372. 		{
  373. 			if(j==p)
  374. 				c[i][p]=v[i][p]*costheta-v[i][q]*sintheta;
  375. 			else if(j==q)
  376. 				c[i][q]=v[i][p]*sintheta+v[i][q]*costheta;
  377. 			else
  378. 				c[i][j]=v[i][j];
  379. 		}
  380. 	}
  381.  
  382. 	//copy c to v for next iteration
  383. 	for(i=0;i<dim;i++)
  384. 	{
  385. 		for(j=0;j<dim;j++)
  386. 		v[i][j]=c[i][j];
  387. 	}
  388. 	//copy b to a for next iteration
  389.  
  390. 	for(i=0;i<dim;i++)
  391. 	{
  392. 		for(j=0;j<dim;j++)
  393. 			a[i][j]=b[i][j];
  394. 	}
  395. 	printf("iter=%d\n",count);
  396. 	printf("p=%d\tq=%d\n",p,q);
  397. 	printf("maximum element of the upper triangular matrix=%f\n",a[p][q]);
  398. 	printf("sintheta=%f\tcostheta=%f\n\n",sintheta,costheta);
  399.  
  400. 	printf("A matrix after rotation\n");
  401.  
  402. 	for(i=0;i<dim;i++)
  403. 	{
  404. 		for(j=0;j<dim;j++)
  405. 			printf("%f\t",a[i][j]);
  406. 		printf("\n");
  407. 	}
  408. 	printf("C matrix after rotation\n");
  409.  
  410. 	for(i=0;i<dim;i++)
  411. 	{
  412. 		for(j=0;j<dim;j++)
  413. 			printf("%f\t",c[i][j]);
  414. 		printf("\n");
  415. 	}
  416. 	printf("\n");
  417.  
  418. //checking for looping for next iteration
  419.  
  420. 	for(i=0;i<dim;i++)
  421. 	{
  422. 		for(j=i+1;j<dim;j++)
  423. 		{
  424. 			if(fabs(b[i][j])!=0.0)
  425. 			{
  426. 				myflag=1;
  427. 				break;
  428. 			}
  429. 			else if(fabs(b[i][j])==0.0)
  430. 				myflag=0;
  431. 		}//j loop
  432. 		if(myflag==1)
  433. 			break;
  434. 	}//i loop
  435. 	if(myflag==0)
  436. 		break;
  437.  
  438. }//while loop
  439. for(i=0;i<dim;i++)
  440. 	{
  441. 		for(j=0;j<dim;j++)
  442. 			fprintf(val,"%lf\t",a[i][j]);
  443. 		fprintf(val,"\n");
  444. 	}
  445. for(i=0;i<dim;i++)
  446. 	{
  447. 		for(j=0;j<dim;j++)
  448. 			fprintf(vec,"%lf\t",c[i][j]);
  449. 		fprintf(vec,"\n");
  450. 	}
  451.  
  452.  
  453. // eval matrix
  454.  
  455. for(i=0;i<dim;i++)
  456. 	for(j=0;j<dim;j++)
  457. 		return(a[i][j]);
  458. // e vector matrix
  459.  
  460. for(i=0;i<dim;i++)
  461. 	for(j=0;j<dim;j++)
  462. 		return(c[i][j]);
  463. }	
  464.  
  465. // combining e val e vect matrices
  466. double combine(int dimn, double u[][n],double v[][n])
  467. {	
  468. 	for(i=0;i<dimn+1;i++)
  469. 	{
  470. 		for(j=0;j<dimn;j++)
  471. 		{
  472. 			combi[i][j]=v[i][j];
  473. 			combi[dimn][j]=u[j][j];
  474. 		}	
  475. 	}
  476. 	printf("\neigen value,eigen vector matrices combined>>>>>>>>>>>\n");
  477. 	for(i=0;i<dimn+1;i++)
  478. 		for(j=0;j<dimn;j++)
  479. 			return(combi[i][j]);
  480. }
  481. // sort()
  482. double sort(int dimn,double w[][n])
  483. {
  484. 	for (j=0; j<dimn; j++)
  485. 	{
  486. 		for (i=0; i<dimn+1; i++)
  487. 		{  
  488. 			if (combi[dimn][j+1] > combi[dimn][j]) 
  489. 			{   
  490. 			temp = combi[i][j];
  491. 			combi[i][j] = combi[i][j+1];
  492. 			combi[i][j+1]= temp;
  493. 			}
  494. 		}
  495. 	}
  496. 	for(i=0;i<dimn;i++)
  497. 		for(j=0;j<dimn;j++)
  498. 			evectsort[i][j]=combi[i][j];
  499. 	print(evectsort,n,n);
  500. 	printf(">>>>>>\n");
  501. // sorting eigen value matrix,a[][], in descending order of diagonal elements
  502. 	for (i=0; i<dimn-1; i++)
  503. 	{
  504. 		for (j=0; j<dimn-1-i; j++)
  505. 		if (a[j+1][j+1] > a[j][j]) 
  506. 		{  
  507. 		tmp = a[j][j];         
  508. 		a[j][j] = a[j+1][j+1];
  509. 		a[j+1][j+1] = tmp;
  510. 		}
  511. 	}
  512. 	print(a,n,n);
  513. 	for(i=0;i<dimn;i++)
  514. 		for(j=0;j<dimn;j++)
  515. 			return(evectsort[i][j]);
  516. 	for(i=0;i<dimn;i++)
  517. 		for(j=0;j<dimn;j++)
  518. 			return(a[i][j]);
  519. }
  520. // loading
  521. double loading(int dimn,double w[][n],double t[][n])// w is a[][],t is evectsort[][]
  522. {
  523. // calculating totvar
  524. 	totvar=0.0;
  525. 	for(i=0;i<dimn;i++)
  526. 	{
  527. 		totvar=totvar+w[i][i];//*a[i][i];
  528. 	}
  529. 	printf("\ntotvar=%lf\n",totvar);
  530. getch();
  531.  
  532. // selecting loading vectors
  533. 	k=0;
  534. 	while(pervar<=pertar && k<dim) 
  535. 	{		
  536. 		pervar1=0.0;
  537. 		for(i=0;i<k+1;i++)
  538. 		{	
  539. 			pervar1=pervar1+w[i][i];//*a[i][i];
  540. 			pervar=pervar1/totvar;
  541. 		}
  542. 		printf("\nk=%d\tpervar=%lf\tpertar=%lf\n",k,pervar,pertar);
  543. 		k=k+1;
  544. 	} 
  545.  
  546. 	edim=k-1;
  547. 	printf("edim=%d",edim);
  548. // eigen vectors(from evectsort[][]) from col 0 to col k-1 are chosen as the loading vectors******dim by edim
  549. 	for(i=0;i<dimn;i++)
  550. 	{
  551. 		for(j=0;j<edim;j++)
  552. 		{
  553. 			sel_evect[i][j]=t[i][j];
  554. 		}
  555. 	}
  556. 	for(i=0;i<dimn;i++)
  557. 		for(j=0;j<edim;j++)
  558. 			return(sel_evect[i][j]);
  559. }
  560. // fscores
  561. double fscores()
  562. {	
  563. 	multiply(m,n,edim,mat,sel_evect);
  564. 	print(mul,m,edim);
  565. 	for(i=0;i<m;i++)
  566. 	{
  567. 		for(j=0;j<edim;j++)
  568. 		{
  569. 			fprintf(scores,"%lf\t",mul[i][j]);
  570. 		}
  571. 		fprintf(scores,"\n");
  572. 	}
  573. 	printf("*******************scores to file done\n");
  574. 	for(i=0;i<m;i++)
  575. 		for(j=0;j<edim;j++)
  576. 			return(mul[i][j]);
  577. }
  578. // multiply
  579. double multiply(int g,int h,int r,double w[][n],double l[][n])//multiply two matrices
  580. {
  581. 	for(i = 0; i < g; i++) 
  582. 		for( j = 0; j < r; j++)
  583.             for( k = 0; k < h; k++) 
  584.                   mul[i][j] +=  w[i][k] * l[k][j];
  585. }
  586. double computesigma(int dimn,double w[][n])// dimn is the edim,w is a[][]
  587. // compute sigma[][]**edim by edim=eigen value matrix corresponding to the selected eigen vectors
  588. {	
  589. 	for(i=0;i<dimn;i++)
  590. 		for(j=0;j<dimn;j++)
  591. 			sigma[i][j]=a[i][j];
  592. 	for(i=0;i<dimn;i++)
  593. 		for(j=0;j<dimn;j++)
  594. 			return(sigma[i][j]);
  595. }
  596. double sigmainverse(int r) // r is dimension of edim
  597. {	
  598. 	matinverse(edim,sigma);
  599. 	for(i=0;i<r;i++)
  600. 		for(j=0;j<r;j++)
  601. 			return(siginv[i][j]);
  602. }
  603.  
  604. double matinverse(int r,double w[][n]) // r is dimension of sigma,w is sigma[][]
  605. {
  606. 	for(i=0;i<r;i++)
  607. 		siginv[i][i]=1;
  608.  
  609. 	for(k=0;k<r-1;k++)
  610. 	{
  611. 		ix=w[k][k];
  612. 		for(i=0;i<r;i++)
  613. 		{
  614. 			w[k][i]=w[k][i]/ix;
  615. 			siginv[k][i]=siginv[k][i]/ix;
  616. 		}
  617. 		for(i=k+1;i<r;i++)
  618. 			im[i-k-1]=w[i][k]/w[k][k];
  619. 		for(i=k+1;i<r;i++)
  620. 			for(j=0;j<r;j++)
  621. 			{
  622. 				w[i][j]-=w[k][j]*im[i-k-1];
  623. 				siginv[i][j]-=siginv[k][j]*im[i-k-1];
  624. 			}
  625. 		det*=ix;
  626. 	}
  627. 	det*=w[k][k];
  628. 	if(det==0)
  629. 		printf("The inverse doesn't exist as det is zero\n");
  630. 	else
  631. 	{
  632. 		for(j=0;j<r;j++)
  633. 			siginv[k][j]/=w[k][k];
  634. 		w[k][k]=1;
  635. 		for(k=edim-1;k>0;k--)
  636. 		{
  637. 			for(i=k;i<r;i++)
  638. 				im[i-k]=w[k-1][i]/w[i][i];
  639. 			for(i=k;i<r;i++)
  640. 				for(j=0;j<r;j++)
  641. 				{
  642. 					w[k-1][j]-=w[i][j]*im[i-k];
  643. 					siginv[k-1][j]-=siginv[i][j]*im[i-k];
  644. 				}
  645. 		}
  646. 	for(i=0;i<r;i++)
  647. 		for(j=0;j<r;j++)
  648. 			return(siginv[i][j]);
  649. 	}
  650. }
  651.  
  652. double ftousqr()
  653. // compute tousqr***tousqr=
  654. //(transpose of x)*(loading vectors)*(sigmainv)*(transpose of loading vectors)*(x)
  655. {
  656. 	for(h=0;h<m;h++)
  657. 	{
  658. 		for(i=h;i<h+1;i++)
  659. 		{
  660. 		  for(j=0;j<n;j++)
  661. 		  {
  662. 			y[i][j]=data[i][j];
  663.  
  664. 		  }
  665.  
  666. 		}
  667. 	// multiply y[1][n]*sel_evect*******1 by edim
  668. 		for(i = h; i < h+1; i++) 
  669.            for( j = 0; j < edim; j++)
  670.              for( k = 0; k < n; k++) 
  671.                    tl[i][j] +=  y[i][k] * sel_evect[k][j];
  672. 	// multiply tl[][]*siginv[][]
  673. 		for(i = h; i < h+1; i++) 
  674.            for( j = 0; j < edim; j++)
  675.              for( k = 0; k < edim; k++) 
  676.                    tls[i][j] +=  tl[i][k] * siginv[k][j];
  677. 		// compute transpose of sel_evect
  678. 		for(i=0;i<edim;i++)
  679. 			for(j=0;j<dim;j++)
  680. 				transload[i][j]=sel_evect[j][i];
  681.  
  682. 		// multiply tls[][]*transload
  683. 		for(i = h; i < h+1; i++) 
  684. 			for( j = 0; j < dim; j++)
  685. 				for( k = 0; k < edim; k++) 
  686. 					 tlstl[i][j] += tls[i][k] * transload[k][j];
  687. 	// transpose of y[][]****** n by 1
  688. 		for(i=0;i<n;i++)
  689. 			for(j=h;j<h+1;j++)
  690. 				ytrans[i][j]=y[j][i];
  691. 	// multiply tlstl[][]*data[][]
  692. 		for(i = h; i < h+1; i++) 
  693. 			for( j = 0; j < dim; j++)
  694. 				for( k = 0; k < 1; k++) 
  695. 					 tousqr[i][j] += tlstl[i][k] * ytrans[k][j];
  696.  
  697. 		printf("tousqr[][]>>>>>>>>>>>>>>>>>\n");
  698. 		for(i=h;i<h+1;i++)
  699. 		{
  700. 		  for(j=0;j<1;j++)
  701. 		  {
  702. 			printf("%lf\n",tousqr[i][j]);
  703. 			fprintf(tou,"%lf\n",tousqr[i][j]);
  704. 		  }
  705. 		  printf("\n");
  706. 		  fprintf(tou,"\n");
  707. 		}
  708. 		printf("\n");
  709.  
  710.   }	//h loop
  711. }
  712.  
  713. double mattranspose(int p,int k,double w[][n])//specify dimensions of w[][] 
  714. {
  715. 	for(i = 0; i < k; i++) 
  716.        for( j = 0; j < p; j++)
  717. 		   mattrans[i][j]=w[j][i];
  718. 	for(i = 0; i < k; i++) 
  719.        for( j = 0; j < p; j++)
  720. 		   return(mattrans[i][j]);
  721. }
  722. double estimate(double w[][2],double r[][dim]) // w is t[][],q is transload[][]
  723. {// scores*transload*******************(mby edim) (edim by dim)
  724. 	for(i = 0; i < m; i++) 
  725.        for( j = 0; j < dim; j++)
  726.          for( k = 0; k < edim; k++) 
  727.             scale_esti[i][j] +=  mul[i][k] * transload[k][j];
  728.  
  729. 	for(i=0;i<m;i++)
  730. 		for(j=0;j<dim;j++)
  731. 			return(scale_esti[i][j]);	
  732. }
  733.  
  734. double fspe()
  735. {// compute Scaled-scale_esti
  736. 	for(i=0;i<m;i++)
  737. 		for(j=0;j<n;j++)
  738. 			zed[i][j]=mat[i][j]-scale_esti[i][j];
  739. // compute spe	
  740. 	for(h=0;h<m;h++)
  741. 	{
  742. 		for(i=h;i<h+1;i++)
  743. 		{
  744. 		  for(j=0;j<n;j++)
  745. 		  {
  746. 			zedy[i][j]=zed[i][j];
  747. 			//printf("%lf\t",zedy[i][j]);
  748. 		  }
  749. 		  printf("\n");
  750. 		}
  751. 		printf("\n");
  752. // transpose of y[][]****** n by 1
  753. 		for(i=0;i<n;i++)
  754. 			for(j=h;j<h+1;j++)
  755. 				zedytrans[i][j]=zedy[j][i];
  756.  
  757. // multiply tlstl[][]*data[][]
  758. 		for(i = h; i < h+1; i++) 
  759. 			for( j = 0; j < dim; j++)
  760. 				for( k = 0; k < 1; k++) 
  761. 					 spe[i][j] += zedy[i][k] * zedytrans[k][j];
  762.  
  763. 		printf("spe[][]>>>>>>>>>>>>>>>>>\n");
  764. 		for(i=h;i<h+1;i++)
  765. 		{
  766. 		  for(j=0;j<1;j++)
  767. 		  {
  768. 			printf("%lf\n",spe[i][j]);
  769. 		  }
  770. 		  printf("\n");
  771. 		}
  772. 		printf("\n");
  773.  
  774. 		for(i=h;i<h+1;i++)
  775. 		{
  776. 		  for(j=0;j<1;j++)
  777. 		  {
  778. 			fprintf(sperror,"%lf\n",spe[i][j]);
  779. 		  }
  780. 		  printf(tou,"\n");
  781. 		}
  782.   }	//h loop	
  783. }
  784. void newline()
  785. {
  786. 	printf("\n\n");
  787. }
Reply With Quote Quick reply to this message  
•
•
Join Date: Jul 2009
Posts: 201
Reputation: yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold yellowSnow is a splendid one to behold 
Solved Threads: 35
yellowSnow's Avatar
yellowSnow yellowSnow is offline Offline
Posting Whiz in Training

Re: pca in c

 
3
  #2
Sep 20th, 2009
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".
Reply With Quote Quick reply to this message  
•
•
Join Date: Feb 2009
Posts: 1,968
Reputation: tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute tux4life has a reputation beyond repute 
Solved Threads: 214
tux4life's Avatar
tux4life tux4life is offline Offline
Posting Virtuoso

Re: pca in c

 
1
  #3
Sep 20th, 2009
In addition (personally one of my favourite links):
http://cboard.cprogramming.com/c-pro...t-process.html
"Never argue with idiots, they just drag you down to their level and then beat you with experience."
Reply With Quote Quick reply to this message  
Reply

Quick Reply to Thread
•
•
•
•
•
•
•
•
Message:
•
•
•
•



Similar Threads
Other Threads in the C Forum
Thread Tools Search this Thread



#include adobe ansi api array asterisks binarysearch changingto char character cm copyimagefile cprogramme creafecopyofanytypeoffileinc createcopyoffile csyntax database directory dynamic execv feet fgets file fork forloop frequency function getlasterror givemetehcodez global grade graphics gtkgcurlcompiling hacking hardware highest histogram i/o include incrementoperators infiniteloop input interest kernel keyboard kilometer license linked linkedlist linux linuxsegmentationfault list locate logical_drives looping loopinsideloop. lowest match matrix meter microsoft motherboard mqqueue mysql number odf opensource owf pattern pdf performance pointer posix probleminc process program programming radix recursion recv repetition research reversing scanf segmentationfault sequential shape socket socketprograming standard string systemcall threads turboc unix user voidmain() wab windows.h windowsapi
About Us | Contact Us | Advertise | DaniWeb | Acceptable Use Policy | RSS Feed

©2003 - 2009 DaniWeb® LLC