2 * Principal component analysis
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
37 PCA
*ff_pca_init(int n
){
42 pca
= av_mallocz(sizeof(PCA
));
45 pca
->covariance
= av_mallocz(sizeof(double)*n
*n
);
46 pca
->mean
= av_mallocz(sizeof(double)*n
);
51 void ff_pca_free(PCA
*pca
){
52 av_freep(&pca
->covariance
);
57 void ff_pca_add(PCA
*pca
, double *v
){
64 pca
->covariance
[j
+ i
*n
] += v
[i
]*v
[j
];
69 int ff_pca(PCA
*pca
, double *eigenvector
, double *eigenvalue
){
74 memset(eigenvector
, 0, sizeof(double)*n
*n
);
77 pca
->mean
[j
] /= pca
->count
;
78 eigenvector
[j
+ j
*n
] = 1.0;
80 pca
->covariance
[j
+ i
*n
] /= pca
->count
;
81 pca
->covariance
[j
+ i
*n
] -= pca
->mean
[i
] * pca
->mean
[j
];
82 pca
->covariance
[i
+ j
*n
] = pca
->covariance
[j
+ i
*n
];
84 eigenvalue
[j
]= pca
->covariance
[j
+ j
*n
];
88 for(pass
=0; pass
< 50; pass
++){
93 sum
+= fabs(pca
->covariance
[j
+ i
*n
]);
99 if(eigenvalue
[j
] > maxvalue
){
100 maxvalue
= eigenvalue
[j
];
104 eigenvalue
[k
]= eigenvalue
[i
];
105 eigenvalue
[i
]= maxvalue
;
107 double tmp
= eigenvector
[k
+ j
*n
];
108 eigenvector
[k
+ j
*n
]= eigenvector
[i
+ j
*n
];
109 eigenvector
[i
+ j
*n
]= tmp
;
116 for(j
=i
+1; j
<n
; j
++){
117 double covar
= pca
->covariance
[j
+ i
*n
];
118 double t
,c
,s
,tau
,theta
, h
;
120 if(pass
< 3 && fabs(covar
) < sum
/ (5*n
*n
)) //FIXME why pass < 3
122 if(fabs(covar
) == 0.0) //FIXME shouldnt be needed
124 if(pass
>=3 && fabs((eigenvalue
[j
]+z
[j
])/covar
) > (1LL<<32) && fabs((eigenvalue
[i
]+z
[i
])/covar
) > (1LL<<32)){
125 pca
->covariance
[j
+ i
*n
]=0.0;
129 h
= (eigenvalue
[j
]+z
[j
]) - (eigenvalue
[i
]+z
[i
]);
131 t
=1.0/(fabs(theta
)+sqrt(1.0+theta
*theta
));
132 if(theta
< 0.0) t
= -t
;
140 #define ROTATE(a,i,j,k,l) {\
141 double g=a[j + i*n];\
142 double h=a[l + k*n];\
143 a[j + i*n]=g-s*(h+g*tau);\
144 a[l + k*n]=h+s*(g-h*tau); }
147 ROTATE(pca
->covariance
,FFMIN(k
,i
),FFMAX(k
,i
),FFMIN(k
,j
),FFMAX(k
,j
))
149 ROTATE(eigenvector
,k
,i
,k
,j
)
151 pca
->covariance
[j
+ i
*n
]=0.0;
154 for (i
=0; i
<n
; i
++) {
155 eigenvalue
[i
] += z
[i
];
174 double eigenvector
[LEN
*LEN
];
175 double eigenvalue
[LEN
];
177 pca
= ff_pca_init(LEN
);
179 for(i
=0; i
<9000000; i
++){
182 int pos
= random()%LEN
;
183 int v2
= (random()%101) - 50;
184 v
[0]= (random()%101) - 50;
186 if(j
<=pos
) v
[j
]= v
[0];
190 /* for(j=0; j<LEN; j++){
193 // sum += random()%10;
194 /* for(j=0; j<LEN; j++){
197 // lbt1(v+100,v+100,LEN);
202 ff_pca(pca
, eigenvector
, eigenvalue
);
203 for(i
=0; i
<LEN
; i
++){
207 // (0.5^|x|)^2 = 0.5^2|x| = 0.25^|x|
210 // pca.covariance[i + i*LEN]= pow(0.5, fabs
211 for(j
=i
; j
<LEN
; j
++){
212 printf("%f ", pca
->covariance
[i
+ j
*LEN
]);
218 for(i
=0; i
<LEN
; i
++){
221 memset(v
, 0, sizeof(v
));
222 for(j
=0; j
<LEN
; j
++){
223 for(k
=0; k
<LEN
; k
++){
224 v
[j
] += pca
->covariance
[FFMIN(k
,j
) + FFMAX(k
,j
)*LEN
] * eigenvector
[i
+ k
*LEN
];
226 v
[j
] /= eigenvalue
[i
];
227 error
+= fabs(v
[j
] - eigenvector
[i
+ j
*LEN
]);
229 printf("%f ", error
);
233 for(i
=0; i
<LEN
; i
++){
234 for(j
=0; j
<LEN
; j
++){
235 printf("%9.6f ", eigenvector
[i
+ j
*LEN
]);
237 printf(" %9.1f %f\n", eigenvalue
[i
], eigenvalue
[i
]/eigenvalue
[0]);