ffmpeg / libavutil / pca.c @ 4869f47e
History  View  Annotate  Download (6.22 KB)
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 021111307 USA

18 
*

19 
*/

20  
21 
/**

22 
* @file pca.c

23 
* Principal component analysis

24 
*/

25  
26 
#include "common.h" 
27 
#include "pca.h" 
28  
29 
typedef struct PCA{ 
30 
int count;

31 
int n;

32 
double *covariance;

33 
double *mean;

34 
}PCA; 
35  
36 
PCA *ff_pca_init(int n){

37 
PCA *pca; 
38 
if(n<=0) 
39 
return NULL; 
40  
41 
pca= av_mallocz(sizeof(PCA));

42 
pca>n= n; 
43 
pca>count=0;

44 
pca>covariance= av_mallocz(sizeof(double)*n*n); 
45 
pca>mean= av_mallocz(sizeof(double)*n); 
46  
47 
return pca;

48 
} 
49  
50 
void ff_pca_free(PCA *pca){

51 
av_freep(&pca>covariance); 
52 
av_freep(&pca>mean); 
53 
av_free(pca); 
54 
} 
55  
56 
void ff_pca_add(PCA *pca, double *v){ 
57 
int i, j;

58 
const int n= pca>n; 
59  
60 
for(i=0; i<n; i++){ 
61 
pca>mean[i] += v[i]; 
62 
for(j=i; j<n; j++)

63 
pca>covariance[j + i*n] += v[i]*v[j]; 
64 
} 
65 
pca>count++; 
66 
} 
67  
68 
int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){ 
69 
int i, j, k, pass;

70 
const int n= pca>n; 
71 
double z[n];

72  
73 
memset(eigenvector, 0, sizeof(double)*n*n); 
74  
75 
for(j=0; j<n; j++){ 
76 
pca>mean[j] /= pca>count; 
77 
eigenvector[j + j*n] = 1.0; 
78 
for(i=0; i<=j; i++){ 
79 
pca>covariance[j + i*n] /= pca>count; 
80 
pca>covariance[j + i*n] = pca>mean[i] * pca>mean[j]; 
81 
pca>covariance[i + j*n] = pca>covariance[j + i*n]; 
82 
} 
83 
eigenvalue[j]= pca>covariance[j + j*n]; 
84 
z[j]= 0;

85 
} 
86  
87 
for(pass=0; pass < 50; pass++){ 
88 
double sum=0; 
89  
90 
for(i=0; i<n; i++) 
91 
for(j=i+1; j<n; j++) 
92 
sum += fabs(pca>covariance[j + i*n]); 
93  
94 
if(sum == 0){ 
95 
for(i=0; i<n; i++){ 
96 
double maxvalue= 1; 
97 
for(j=i; j<n; j++){

98 
if(eigenvalue[j] > maxvalue){

99 
maxvalue= eigenvalue[j]; 
100 
k= j; 
101 
} 
102 
} 
103 
eigenvalue[k]= eigenvalue[i]; 
104 
eigenvalue[i]= maxvalue; 
105 
for(j=0; j<n; j++){ 
106 
double tmp= eigenvector[k + j*n];

107 
eigenvector[k + j*n]= eigenvector[i + j*n]; 
108 
eigenvector[i + j*n]= tmp; 
109 
} 
110 
} 
111 
return pass;

112 
} 
113  
114 
for(i=0; i<n; i++){ 
115 
for(j=i+1; j<n; j++){ 
116 
double covar= pca>covariance[j + i*n];

117 
double t,c,s,tau,theta, h;

118  
119 
if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3 
120 
continue;

121 
if(fabs(covar) == 0.0) //FIXME shouldnt be needed 
122 
continue;

123 
if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){ 
124 
pca>covariance[j + i*n]=0.0; 
125 
continue;

126 
} 
127  
128 
h= (eigenvalue[j]+z[j])  (eigenvalue[i]+z[i]); 
129 
theta=0.5*h/covar; 
130 
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); 
131 
if(theta < 0.0) t = t; 
132  
133 
c=1.0/sqrt(1+t*t); 
134 
s=t*c; 
135 
tau=s/(1.0+c); 
136 
z[i] = t*covar; 
137 
z[j] += t*covar; 
138  
139 
#define ROTATE(a,i,j,k,l) {\

140 
double g=a[j + i*n];\

141 
double h=a[l + k*n];\

142 
a[j + i*n]=gs*(h+g*tau);\ 
143 
a[l + k*n]=h+s*(gh*tau); } 
144 
for(k=0; k<n; k++) { 
145 
if(k!=i && k!=j){

146 
ROTATE(pca>covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j)) 
147 
} 
148 
ROTATE(eigenvector,k,i,k,j) 
149 
} 
150 
pca>covariance[j + i*n]=0.0; 
151 
} 
152 
} 
153 
for (i=0; i<n; i++) { 
154 
eigenvalue[i] += z[i]; 
155 
z[i]=0.0; 
156 
} 
157 
} 
158  
159 
return 1; 
160 
} 
161  
162 
#ifdef TEST

163  
164 
#undef printf

165 
#undef random

166 
#include <stdio.h> 
167 
#include <stdlib.h> 
168  
169 
int main(){

170 
PCA *pca; 
171 
int i, j, k;

172 
#define LEN 8 
173 
double eigenvector[LEN*LEN];

174 
double eigenvalue[LEN];

175  
176 
pca= ff_pca_init(LEN); 
177  
178 
for(i=0; i<9000000; i++){ 
179 
double v[2*LEN+100]; 
180 
double sum=0; 
181 
int pos= random()%LEN;

182 
int v2= (random()%101)  50; 
183 
v[0]= (random()%101)  50; 
184 
for(j=1; j<8; j++){ 
185 
if(j<=pos) v[j]= v[0]; 
186 
else v[j]= v2;

187 
sum += v[j]; 
188 
} 
189 
/* for(j=0; j<LEN; j++){

190 
v[j] = v[pos];

191 
}*/

192 
// sum += random()%10;

193 
/* for(j=0; j<LEN; j++){

194 
v[j] = sum/LEN;

195 
}*/

196 
// lbt1(v+100,v+100,LEN);

197 
ff_pca_add(pca, v); 
198 
} 
199  
200  
201 
ff_pca(pca, eigenvector, eigenvalue); 
202 
for(i=0; i<LEN; i++){ 
203 
pca>count= 1;

204 
pca>mean[i]= 0;

205  
206 
// (0.5^x)^2 = 0.5^2x = 0.25^x

207  
208  
209 
// pca.covariance[i + i*LEN]= pow(0.5, fabs

210 
for(j=i; j<LEN; j++){

211 
printf("%f ", pca>covariance[i + j*LEN]);

212 
} 
213 
printf("\n");

214 
} 
215  
216 
#if 1 
217 
for(i=0; i<LEN; i++){ 
218 
double v[LEN];

219 
double error=0; 
220 
memset(v, 0, sizeof(v)); 
221 
for(j=0; j<LEN; j++){ 
222 
for(k=0; k<LEN; k++){ 
223 
v[j] += pca>covariance[FFMIN(k,j) + FFMAX(k,j)*LEN] * eigenvector[i + k*LEN]; 
224 
} 
225 
v[j] /= eigenvalue[i]; 
226 
error += fabs(v[j]  eigenvector[i + j*LEN]); 
227 
} 
228 
printf("%f ", error);

229 
} 
230 
printf("\n");

231 
#endif

232 
for(i=0; i<LEN; i++){ 
233 
for(j=0; j<LEN; j++){ 
234 
printf("%9.6f ", eigenvector[i + j*LEN]);

235 
} 
236 
printf(" %9.1f %f\n", eigenvalue[i], eigenvalue[i]/eigenvalue[0]); 
237 
} 
238  
239 
return 0; 
240 
} 
241 
#endif
