Mercurial > libavutil.hg
annotate pca.c @ 554:29d015080b3d libavutil
Do not mix declarations and statements (by ramiro).
author | michael |
---|---|
date | Sun, 17 Aug 2008 15:35:55 +0000 |
parents | 288f054e06e8 |
children | 2fb78abcf725 |
rev | line source |
---|---|
550 | 1 /* |
2 * Principal component analysis | |
3 * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at> | |
4 * | |
5 * This library is free software; you can redistribute it and/or | |
6 * modify it under the terms of the GNU Lesser General Public | |
7 * License as published by the Free Software Foundation; either | |
8 * version 2 of the License, or (at your option) any later version. | |
9 * | |
10 * This library is distributed in the hope that it will be useful, | |
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 * Lesser General Public License for more details. | |
14 * | |
15 * You should have received a copy of the GNU Lesser General Public | |
16 * License along with this library; if not, write to the Free Software | |
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
18 * | |
19 */ | |
20 | |
21 /** | |
22 * @file pca.c | |
23 * Principal component analysis | |
24 */ | |
25 | |
551 | 26 #include "common.h" |
550 | 27 #include "pca.h" |
28 | |
29 int ff_pca_init(PCA *pca, int n){ | |
30 if(n<=0) | |
31 return -1; | |
32 | |
33 pca->n= n; | |
34 pca->count=0; | |
35 pca->covariance= av_mallocz(sizeof(double)*n*n); | |
36 pca->mean= av_mallocz(sizeof(double)*n); | |
37 | |
38 return 0; | |
39 } | |
40 | |
41 void ff_pca_free(PCA *pca){ | |
42 av_freep(&pca->covariance); | |
43 av_freep(&pca->mean); | |
44 } | |
45 | |
46 void ff_pca_add(PCA *pca, double *v){ | |
47 int i, j; | |
48 const int n= pca->n; | |
49 | |
50 for(i=0; i<n; i++){ | |
51 pca->mean[i] += v[i]; | |
52 for(j=i; j<n; j++) | |
53 pca->covariance[j + i*n] += v[i]*v[j]; | |
54 } | |
55 pca->count++; | |
56 } | |
57 | |
58 int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){ | |
59 int i, j, k, pass; | |
60 const int n= pca->n; | |
61 double z[n]; | |
62 | |
63 memset(eigenvector, 0, sizeof(double)*n*n); | |
64 | |
65 for(j=0; j<n; j++){ | |
66 pca->mean[j] /= pca->count; | |
67 eigenvector[j + j*n] = 1.0; | |
68 for(i=0; i<=j; i++){ | |
69 pca->covariance[j + i*n] /= pca->count; | |
70 pca->covariance[j + i*n] -= pca->mean[i] * pca->mean[j]; | |
71 pca->covariance[i + j*n] = pca->covariance[j + i*n]; | |
72 } | |
73 eigenvalue[j]= pca->covariance[j + j*n]; | |
74 z[j]= 0; | |
75 } | |
76 | |
77 for(pass=0; pass < 50; pass++){ | |
78 double sum=0; | |
79 | |
80 for(i=0; i<n; i++) | |
81 for(j=i+1; j<n; j++) | |
82 sum += fabs(pca->covariance[j + i*n]); | |
83 | |
84 if(sum == 0){ | |
85 for(i=0; i<n; i++){ | |
86 double maxvalue= -1; | |
87 for(j=i; j<n; j++){ | |
88 if(eigenvalue[j] > maxvalue){ | |
89 maxvalue= eigenvalue[j]; | |
90 k= j; | |
91 } | |
92 } | |
93 eigenvalue[k]= eigenvalue[i]; | |
94 eigenvalue[i]= maxvalue; | |
95 for(j=0; j<n; j++){ | |
96 double tmp= eigenvector[k + j*n]; | |
97 eigenvector[k + j*n]= eigenvector[i + j*n]; | |
98 eigenvector[i + j*n]= tmp; | |
99 } | |
100 } | |
101 return pass; | |
102 } | |
103 | |
104 for(i=0; i<n; i++){ | |
105 for(j=i+1; j<n; j++){ | |
106 double covar= pca->covariance[j + i*n]; | |
107 double t,c,s,tau,theta, h; | |
108 | |
109 if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3 | |
110 continue; | |
111 if(fabs(covar) == 0.0) //FIXME shouldnt be needed | |
112 continue; | |
113 if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){ | |
114 pca->covariance[j + i*n]=0.0; | |
115 continue; | |
116 } | |
117 | |
118 h= (eigenvalue[j]+z[j]) - (eigenvalue[i]+z[i]); | |
119 theta=0.5*h/covar; | |
120 t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); | |
121 if(theta < 0.0) t = -t; | |
122 | |
123 c=1.0/sqrt(1+t*t); | |
124 s=t*c; | |
125 tau=s/(1.0+c); | |
126 z[i] -= t*covar; | |
127 z[j] += t*covar; | |
128 | |
554
29d015080b3d
Do not mix declarations and statements (by ramiro).
michael
parents:
553
diff
changeset
|
129 #define ROTATE(a,i,j,k,l) {\ |
550 | 130 double g=a[j + i*n];\ |
131 double h=a[l + k*n];\ | |
132 a[j + i*n]=g-s*(h+g*tau);\ | |
554
29d015080b3d
Do not mix declarations and statements (by ramiro).
michael
parents:
553
diff
changeset
|
133 a[l + k*n]=h+s*(g-h*tau); } |
550 | 134 for(k=0; k<n; k++) { |
135 if(k!=i && k!=j){ | |
136 ROTATE(pca->covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j)) | |
137 } | |
138 ROTATE(eigenvector,k,i,k,j) | |
139 } | |
140 pca->covariance[j + i*n]=0.0; | |
141 } | |
142 } | |
143 for (i=0; i<n; i++) { | |
144 eigenvalue[i] += z[i]; | |
145 z[i]=0.0; | |
146 } | |
147 } | |
148 | |
149 return -1; | |
150 } | |
151 | |
553 | 152 #ifdef TEST |
550 | 153 |
154 #undef printf | |
552 | 155 #undef random |
550 | 156 #include <stdio.h> |
157 #include <stdlib.h> | |
158 | |
159 int main(){ | |
160 PCA pca; | |
161 int i, j, k; | |
162 #define LEN 8 | |
163 double eigenvector[LEN*LEN]; | |
164 double eigenvalue[LEN]; | |
165 | |
166 ff_pca_init(&pca, LEN); | |
167 | |
168 for(i=0; i<9000000; i++){ | |
169 double v[2*LEN+100]; | |
170 double sum=0; | |
171 int pos= random()%LEN; | |
172 int v2= (random()%101) - 50; | |
173 v[0]= (random()%101) - 50; | |
174 for(j=1; j<8; j++){ | |
175 if(j<=pos) v[j]= v[0]; | |
176 else v[j]= v2; | |
177 sum += v[j]; | |
178 } | |
179 /* for(j=0; j<LEN; j++){ | |
180 v[j] -= v[pos]; | |
181 }*/ | |
182 // sum += random()%10; | |
183 /* for(j=0; j<LEN; j++){ | |
184 v[j] -= sum/LEN; | |
185 }*/ | |
186 // lbt1(v+100,v+100,LEN); | |
187 ff_pca_add(&pca, v); | |
188 } | |
189 | |
190 | |
191 ff_pca(&pca, eigenvector, eigenvalue); | |
192 for(i=0; i<LEN; i++){ | |
193 pca.count= 1; | |
194 pca.mean[i]= 0; | |
195 | |
196 // (0.5^|x|)^2 = 0.5^2|x| = 0.25^|x| | |
197 | |
198 | |
199 // pca.covariance[i + i*LEN]= pow(0.5, fabs | |
200 for(j=i; j<LEN; j++){ | |
201 printf("%f ", pca.covariance[i + j*LEN]); | |
202 } | |
203 printf("\n"); | |
204 } | |
205 | |
206 #if 1 | |
207 for(i=0; i<LEN; i++){ | |
208 double v[LEN]; | |
209 double error=0; | |
210 memset(v, 0, sizeof(v)); | |
211 for(j=0; j<LEN; j++){ | |
212 for(k=0; k<LEN; k++){ | |
213 v[j] += pca.covariance[FFMIN(k,j) + FFMAX(k,j)*LEN] * eigenvector[i + k*LEN]; | |
214 } | |
215 v[j] /= eigenvalue[i]; | |
216 error += fabs(v[j] - eigenvector[i + j*LEN]); | |
217 } | |
218 printf("%f ", error); | |
219 } | |
220 printf("\n"); | |
221 #endif | |
222 for(i=0; i<LEN; i++){ | |
223 for(j=0; j<LEN; j++){ | |
224 printf("%9.6f ", eigenvector[i + j*LEN]); | |
225 } | |
226 printf(" %9.1f %f\n", eigenvalue[i], eigenvalue[i]/eigenvalue[0]); | |
227 } | |
228 | |
229 return 0; | |
230 } | |
231 #endif |