Principal component analysis
[ffmpeg.git] / libavutil / pca.c
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
26 #include <math.h>
27 #include "avcodec.h"
28 #include "pca.h"
29
30 int ff_pca_init(PCA *pca, int n){
31     if(n<=0)
32         return -1;
33
34     pca->n= n;
35     pca->count=0;
36     pca->covariance= av_mallocz(sizeof(double)*n*n);
37     pca->mean= av_mallocz(sizeof(double)*n);
38
39     return 0;
40 }
41
42 void ff_pca_free(PCA *pca){
43     av_freep(&pca->covariance);
44     av_freep(&pca->mean);
45 }
46
47 void ff_pca_add(PCA *pca, double *v){
48     int i, j;
49     const int n= pca->n;
50
51     for(i=0; i<n; i++){
52         pca->mean[i] += v[i];
53         for(j=i; j<n; j++)
54             pca->covariance[j + i*n] += v[i]*v[j];
55     }
56     pca->count++;
57 }
58
59 int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){
60     int i, j, k, pass;
61     const int n= pca->n;
62     double z[n];
63
64     memset(eigenvector, 0, sizeof(double)*n*n);
65
66     for(j=0; j<n; j++){
67         pca->mean[j] /= pca->count;
68         eigenvector[j + j*n] = 1.0;
69         for(i=0; i<=j; i++){
70             pca->covariance[j + i*n] /= pca->count;
71             pca->covariance[j + i*n] -= pca->mean[i] * pca->mean[j];
72             pca->covariance[i + j*n] = pca->covariance[j + i*n];
73         }
74         eigenvalue[j]= pca->covariance[j + j*n];
75         z[j]= 0;
76     }
77
78     for(pass=0; pass < 50; pass++){
79         double sum=0;
80
81         for(i=0; i<n; i++)
82             for(j=i+1; j<n; j++)
83                 sum += fabs(pca->covariance[j + i*n]);
84
85         if(sum == 0){
86             for(i=0; i<n; i++){
87                 double maxvalue= -1;
88                 for(j=i; j<n; j++){
89                     if(eigenvalue[j] > maxvalue){
90                         maxvalue= eigenvalue[j];
91                         k= j;
92                     }
93                 }
94                 eigenvalue[k]= eigenvalue[i];
95                 eigenvalue[i]= maxvalue;
96                 for(j=0; j<n; j++){
97                     double tmp= eigenvector[k + j*n];
98                     eigenvector[k + j*n]= eigenvector[i + j*n];
99                     eigenvector[i + j*n]= tmp;
100                 }
101             }
102             return pass;
103         }
104
105         for(i=0; i<n; i++){
106             for(j=i+1; j<n; j++){
107                 double covar= pca->covariance[j + i*n];
108                 double t,c,s,tau,theta, h;
109
110                 if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3
111                     continue;
112                 if(fabs(covar) == 0.0) //FIXME shouldnt be needed
113                     continue;
114                 if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){
115                     pca->covariance[j + i*n]=0.0;
116                     continue;
117                 }
118
119                 h= (eigenvalue[j]+z[j]) - (eigenvalue[i]+z[i]);
120                 theta=0.5*h/covar;
121                 t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
122                 if(theta < 0.0) t = -t;
123
124                 c=1.0/sqrt(1+t*t);
125                 s=t*c;
126                 tau=s/(1.0+c);
127                 z[i] -= t*covar;
128                 z[j] += t*covar;
129
130 #define ROTATE(a,i,j,k,l)\
131     double g=a[j + i*n];\
132     double h=a[l + k*n];\
133     a[j + i*n]=g-s*(h+g*tau);\
134     a[l + k*n]=h+s*(g-h*tau);
135                 for(k=0; k<n; k++) {
136                     if(k!=i && k!=j){
137                         ROTATE(pca->covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j))
138                     }
139                     ROTATE(eigenvector,k,i,k,j)
140                 }
141                 pca->covariance[j + i*n]=0.0;
142             }
143         }
144         for (i=0; i<n; i++) {
145             eigenvalue[i] += z[i];
146             z[i]=0.0;
147         }
148     }
149
150     return -1;
151 }
152
153 #if 1
154
155 #undef printf
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