| 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, const 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 |