Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
4349 | Serge | 1 | /* |
2 | * principal component analysis (PCA) |
||
3 | * Copyright (c) 2004 Michael Niedermayer |
||
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(sizeof(*pca->z) * n); |
||
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 |
||
65 | pca->mean[i] += v[i]; |
||
66 | for(j=i; 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 |
||
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 |
||
96 | for(j=i+1; j |
||
97 | sum += fabs(pca->covariance[j + i*n]); |
||
98 | |||
99 | if(sum == 0){ |
||
100 | for(i=0; i |
||
101 | double maxvalue= -1; |
||
102 | for(j=i; 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 |
||
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 |
||
120 | for(j=i+1; 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 |
||
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 |
||
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 |
||
171 | #include |
||
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 |
||
198 | v[j] -= v[pos]; |
||
199 | }*/ |
||
200 | // sum += av_lfg_get(&prng) % 10; |
||
201 | /* for(j=0; 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 |
||
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 |
||
219 | printf("%f ", pca->covariance[i + j*LEN]); |
||
220 | } |
||
221 | printf("\n"); |
||
222 | } |
||
223 | |||
224 | for(i=0; i |
||
225 | double v[LEN]; |
||
226 | double error=0; |
||
227 | memset(v, 0, sizeof(v)); |
||
228 | for(j=0; j |
||
229 | for(k=0; 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 |
||
240 | for(j=0; 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=pos)>8;>9000000;>>32)){ |