f7ae350b48c73732886d6dcbe0c9ddd4808d3cdc
2 * principal component analysis (PCA)
3 * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>
5 * This file is part of FFmpeg.
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.
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.
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
24 * principal component analysis (PCA)
38 PCA
*ff_pca_init(int n
){
43 pca
= av_mallocz(sizeof(*pca
));
45 pca
->z
= av_malloc_array(n
, sizeof(*pca
->z
));
47 pca
->covariance
= av_calloc(n
*n
, sizeof(double));
48 pca
->mean
= av_calloc(n
, sizeof(double));
53 void ff_pca_free(PCA
*pca
){
54 av_freep(&pca
->covariance
);
60 void ff_pca_add(PCA
*pca
, const double *v
){
67 pca
->covariance
[j
+ i
*n
] += v
[i
]*v
[j
];
72 int ff_pca(PCA
*pca
, double *eigenvector
, double *eigenvalue
){
78 memset(eigenvector
, 0, sizeof(double)*n
*n
);
81 pca
->mean
[j
] /= pca
->count
;
82 eigenvector
[j
+ j
*n
] = 1.0;
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
];
88 eigenvalue
[j
]= pca
->covariance
[j
+ j
*n
];
92 for(pass
=0; pass
< 50; pass
++){
97 sum
+= fabs(pca
->covariance
[j
+ i
*n
]);
103 if(eigenvalue
[j
] > maxvalue
){
104 maxvalue
= eigenvalue
[j
];
108 eigenvalue
[k
]= eigenvalue
[i
];
109 eigenvalue
[i
]= maxvalue
;
111 double tmp
= eigenvector
[k
+ j
*n
];
112 eigenvector
[k
+ j
*n
]= eigenvector
[i
+ j
*n
];
113 eigenvector
[i
+ j
*n
]= tmp
;
120 for(j
=i
+1; j
<n
; j
++){
121 double covar
= pca
->covariance
[j
+ i
*n
];
122 double t
,c
,s
,tau
,theta
, h
;
124 if(pass
< 3 && fabs(covar
) < sum
/ (5*n
*n
)) //FIXME why pass < 3
126 if(fabs(covar
) == 0.0) //FIXME should not be needed
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;
133 h
= (eigenvalue
[j
]+z
[j
]) - (eigenvalue
[i
]+z
[i
]);
135 t
=1.0/(fabs(theta
)+sqrt(1.0+theta
*theta
));
136 if(theta
< 0.0) t
= -t
;
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); }
151 ROTATE(pca
->covariance
,FFMIN(k
,i
),FFMAX(k
,i
),FFMIN(k
,j
),FFMAX(k
,j
))
153 ROTATE(eigenvector
,k
,i
,k
,j
)
155 pca
->covariance
[j
+ i
*n
]=0.0;
158 for (i
=0; i
<n
; i
++) {
159 eigenvalue
[i
] += z
[i
];
178 double eigenvector
[LEN
*LEN
];
179 double eigenvalue
[LEN
];
182 av_lfg_init(&prng
, 1);
184 pca
= ff_pca_init(LEN
);
186 for(i
=0; i
<9000000; i
++){
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;
193 if(j
<=pos
) v
[j
]= v
[0];
197 /* for(j=0; j<LEN; j++){
200 // sum += av_lfg_get(&prng) % 10;
201 /* for(j=0; j<LEN; j++){
204 // lbt1(v+100,v+100,LEN);
209 ff_pca(pca
, eigenvector
, eigenvalue
);
210 for(i
=0; i
<LEN
; i
++){
214 // (0.5^|x|)^2 = 0.5^2|x| = 0.25^|x|
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
]);
224 for(i
=0; i
<LEN
; i
++){
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
];
232 v
[j
] /= eigenvalue
[i
];
233 error
+= fabs(v
[j
] - eigenvector
[i
+ j
*LEN
]);
235 printf("%f ", error
);
239 for(i
=0; i
<LEN
; i
++){
240 for(j
=0; j
<LEN
; j
++){
241 printf("%9.6f ", eigenvector
[i
+ j
*LEN
]);
243 printf(" %9.1f %f\n", eigenvalue
[i
], eigenvalue
[i
]/eigenvalue
[0]);