annotate lls.c @ 992:a13125b5be3a libavutil

bswap: change ME to NE in macro names Other parts of FFmpeg use NE (native endian) rather than ME (machine). This makes it consistent.
author mru
date Sat, 10 Jul 2010 22:09:01 +0000
parents 0795a743bda1
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
1 /*
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
2 * linear least squares model
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
3 *
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
4 * Copyright (c) 2006 Michael Niedermayer <michaelni@gmx.at>
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
5 *
116
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
6 * This file is part of FFmpeg.
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
7 *
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
8 * FFmpeg is free software; you can redistribute it and/or
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
9 * modify it under the terms of the GNU Lesser General Public
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
10 * License as published by the Free Software Foundation; either
116
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
11 * version 2.1 of the License, or (at your option) any later version.
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
12 *
116
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
13 * FFmpeg is distributed in the hope that it will be useful,
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
16 * Lesser General Public License for more details.
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
17 *
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
18 * You should have received a copy of the GNU Lesser General Public
116
d76a36742464 Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents: 90
diff changeset
19 * License along with FFmpeg; if not, write to the Free Software
358
f13e5473611e license header consistency cosmetics
diego
parents: 116
diff changeset
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
21 */
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
22
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
23 /**
899
0795a743bda1 Remove explicit filename from Doxygen @file commands.
diego
parents: 763
diff changeset
24 * @file
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
25 * linear least squares model
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
26 */
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
27
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
28 #include <math.h>
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
29 #include <string.h>
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
30
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
31 #include "lls.h"
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
32
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
33 void av_init_lls(LLSModel *m, int indep_count){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
34 memset(m, 0, sizeof(LLSModel));
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
35
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
36 m->indep_count= indep_count;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
37 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
38
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
39 void av_update_lls(LLSModel *m, double *var, double decay){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
40 int i,j;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
41
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
42 for(i=0; i<=m->indep_count; i++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
43 for(j=i; j<=m->indep_count; j++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
44 m->covariance[i][j] *= decay;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
45 m->covariance[i][j] += var[i]*var[j];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
46 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
47 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
48 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
49
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
50 void av_solve_lls(LLSModel *m, double threshold, int min_order){
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
51 int i,j,k;
441
6335e6e63d91 Fix the following using void* casts, proper casts are less readable and
michael
parents: 424
diff changeset
52 double (*factor)[MAX_VARS+1]= (void*)&m->covariance[1][0];
6335e6e63d91 Fix the following using void* casts, proper casts are less readable and
michael
parents: 424
diff changeset
53 double (*covar )[MAX_VARS+1]= (void*)&m->covariance[1][1];
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
54 double *covar_y = m->covariance[0];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
55 int count= m->indep_count;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
56
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
57 for(i=0; i<count; i++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
58 for(j=i; j<count; j++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
59 double sum= covar[i][j];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
60
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
61 for(k=i-1; k>=0; k--)
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
62 sum -= factor[i][k]*factor[j][k];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
63
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
64 if(i==j){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
65 if(sum < threshold)
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
66 sum= 1.0;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
67 factor[i][i]= sqrt(sum);
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
68 }else
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
69 factor[j][i]= sum / factor[i][i];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
70 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
71 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
72 for(i=0; i<count; i++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
73 double sum= covar_y[i+1];
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
74 for(k=i-1; k>=0; k--)
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
75 sum -= factor[i][k]*m->coeff[0][k];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
76 m->coeff[0][i]= sum / factor[i][i];
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
77 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
78
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
79 for(j=count-1; j>=min_order; j--){
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
80 for(i=j; i>=0; i--){
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
81 double sum= m->coeff[0][i];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
82 for(k=i+1; k<=j; k++)
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
83 sum -= factor[k][i]*m->coeff[j][k];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
84 m->coeff[j][i]= sum / factor[i][i];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
85 }
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
86
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
87 m->variance[j]= covar_y[0];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
88 for(i=0; i<=j; i++){
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
89 double sum= m->coeff[j][i]*covar[i][i] - 2*covar_y[i+1];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
90 for(k=0; k<i; k++)
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
91 sum += 2*m->coeff[j][k]*covar[k][i];
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
92 m->variance[j] += m->coeff[j][i]*sum;
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
93 }
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
94 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
95 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
96
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
97 double av_evaluate_lls(LLSModel *m, double *param, int order){
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
98 int i;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
99 double out= 0;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
100
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
101 for(i=0; i<=order; i++)
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
102 out+= param[i]*m->coeff[order][i];
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
103
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
104 return out;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
105 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
106
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
107 #ifdef TEST
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
108
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
109 #include <stdlib.h>
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
110 #include <stdio.h>
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
111
404
f9a4c04ebb0e main() --> main(void)
diego
parents: 358
diff changeset
112 int main(void){
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
113 LLSModel m;
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
114 int i, order;
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
115
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
116 av_init_lls(&m, 3);
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
117
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
118 for(i=0; i<100; i++){
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
119 double var[4];
424
1cdbf12cb116 Remove unused variable variance.
diego
parents: 404
diff changeset
120 double eval;
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
121 var[0] = (rand() / (double)RAND_MAX - 0.5)*2;
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
122 var[1] = var[0] + rand() / (double)RAND_MAX - 0.5;
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
123 var[2] = var[1] + rand() / (double)RAND_MAX - 0.5;
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
124 var[3] = var[2] + rand() / (double)RAND_MAX - 0.5;
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
125 av_update_lls(&m, var, 0.99);
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
126 av_solve_lls(&m, 0.001, 0);
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
127 for(order=0; order<3; order++){
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
128 eval= av_evaluate_lls(&m, var+1, order);
700
e110508543ac Align test program output columns.
diego
parents: 642
diff changeset
129 printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n",
79
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
130 var[0], order, eval, sqrt(m.variance[order] / (i+1)),
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
131 m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]);
adbb5540fa47 calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders
michael
parents: 78
diff changeset
132 }
76
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
133 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
134 return 0;
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
135 }
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
136
8c75234388b5 linear least squares solver using cholesky factorization
michael
parents:
diff changeset
137 #endif