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
|
|
129 #define ROTATE(a,i,j,k,l)\
|
|
130 double g=a[j + i*n];\
|
|
131 double h=a[l + k*n];\
|
|
132 a[j + i*n]=g-s*(h+g*tau);\
|
|
133 a[l + k*n]=h+s*(g-h*tau);
|
|
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
|