Mercurial > libavutil.hg
comparison lls.c @ 76:8c75234388b5 libavutil
linear least squares solver using cholesky factorization
author | michael |
---|---|
date | Fri, 14 Jul 2006 10:03:09 +0000 |
parents | |
children | 1eb5129d5030 |
comparison
equal
deleted
inserted
replaced
75:09ed77540442 | 76:8c75234388b5 |
---|---|
1 /* | |
2 * linear least squares model | |
3 * | |
4 * Copyright (c) 2006 Michael Niedermayer <michaelni@gmx.at> | |
5 * | |
6 * This library is free software; you can redistribute it and/or | |
7 * modify it under the terms of the GNU Lesser General Public | |
8 * License as published by the Free Software Foundation; either | |
9 * version 2 of the License, or (at your option) any later version. | |
10 * | |
11 * This library is distributed in the hope that it will be useful, | |
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 * Lesser General Public License for more details. | |
15 * | |
16 * You should have received a copy of the GNU Lesser General Public | |
17 * License along with this library; if not, write to the Free Software | |
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
19 */ | |
20 | |
21 /** | |
22 * @file lls.c | |
23 * linear least squares model | |
24 */ | |
25 | |
26 #include <math.h> | |
27 #include <string.h> | |
28 | |
29 #include "lls.h" | |
30 | |
31 #undef NDEBUG // allways check asserts, the speed effect is far too small to disable them | |
32 #include <assert.h> | |
33 | |
34 #ifdef TEST | |
35 #define av_log(a,b,...) printf(__VA_ARGS__) | |
36 #endif | |
37 | |
38 void av_init_lls(LLSModel *m, int indep_count){ | |
39 memset(m, 0, sizeof(LLSModel)); | |
40 | |
41 m->indep_count= indep_count; | |
42 } | |
43 | |
44 void av_update_lls(LLSModel *m, double *var, double decay){ | |
45 int i,j; | |
46 | |
47 for(i=0; i<=m->indep_count; i++){ | |
48 for(j=i; j<=m->indep_count; j++){ | |
49 m->covariance[i][j] *= decay; | |
50 m->covariance[i][j] += var[i]*var[j]; | |
51 } | |
52 } | |
53 } | |
54 | |
55 double av_solve_lls(LLSModel *m, double threshold){ | |
56 int i,j,k; | |
57 double (*factor)[MAX_VARS+1]= &m->covariance[1][0]; | |
58 double (*covar )[MAX_VARS+1]= &m->covariance[1][1]; | |
59 double *covar_y = m->covariance[0]; | |
60 double variance; | |
61 int count= m->indep_count; | |
62 | |
63 for(i=0; i<count; i++){ | |
64 for(j=i; j<count; j++){ | |
65 double sum= covar[i][j]; | |
66 | |
67 for(k=i-1; k>=0; k--) | |
68 sum -= factor[i][k]*factor[j][k]; | |
69 | |
70 if(i==j){ | |
71 if(sum < threshold) | |
72 sum= 1.0; | |
73 factor[i][i]= sqrt(sum); | |
74 }else | |
75 factor[j][i]= sum / factor[i][i]; | |
76 } | |
77 } | |
78 for(i=0; i<count; i++){ | |
79 double sum= covar_y[i+1]; | |
80 for(k=i-1; k>=0; k--) | |
81 sum -= factor[i][k]*m->coeff[k]; | |
82 m->coeff[i]= sum / factor[i][i]; | |
83 } | |
84 | |
85 for(i=count-1; i>=0; i--){ | |
86 double sum= m->coeff[i]; | |
87 for(k=i+1; k<count; k++) | |
88 sum -= factor[k][i]*m->coeff[k]; | |
89 m->coeff[i]= sum / factor[i][i]; | |
90 } | |
91 | |
92 variance= covar_y[0]; | |
93 for(i=0; i<count; i++){ | |
94 double sum= m->coeff[i]*covar[i][i] - 2*covar_y[i+1]; | |
95 for(j=0; j<i; j++) | |
96 sum += 2*m->coeff[j]*covar[j][i]; | |
97 variance += m->coeff[i]*sum; | |
98 } | |
99 return variance; | |
100 } | |
101 | |
102 double av_evaluate_lls(LLSModel *m, double *param){ | |
103 int i; | |
104 double out= 0; | |
105 | |
106 for(i=0; i<m->indep_count; i++) | |
107 out+= param[i]*m->coeff[i]; | |
108 | |
109 return out; | |
110 } | |
111 | |
112 #ifdef TEST | |
113 | |
114 #include <stdlib.h> | |
115 #include <stdio.h> | |
116 | |
117 int main(){ | |
118 LLSModel m; | |
119 int i; | |
120 | |
121 av_init_lls(&m, 3); | |
122 | |
123 for(i=0; i<100; i++){ | |
124 double var[4]; | |
125 double eval, variance; | |
126 var[1] = rand() / (double)RAND_MAX; | |
127 var[2] = rand() / (double)RAND_MAX; | |
128 var[3] = rand() / (double)RAND_MAX; | |
129 | |
130 var[2]= var[1] + var[3]; | |
131 | |
132 var[0] = var[1] + var[2] + var[3] + var[1]*var[2]/100; | |
133 | |
134 eval= av_evaluate_lls(&m, var+1); | |
135 av_update_lls(&m, var, 0.99); | |
136 variance= av_solve_lls(&m, 0.001); | |
137 av_log(NULL, AV_LOG_DEBUG, "real:%f pred:%f var:%f coeffs:%f %f %f\n", | |
138 var[0], eval, sqrt(variance / (i+1)), | |
139 m.coeff[0], m.coeff[1], m.coeff[2]); | |
140 } | |
141 return 0; | |
142 } | |
143 | |
144 #endif |