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