ffmpeg / libavutil / lls.c @ b2755007
History  View  Annotate  Download (3.98 KB)
1 
/*


2 
* linear least squares model

3 
*

4 
* Copyright (c) 2006 Michael Niedermayer <michaelni@gmx.at>

5 
*

6 
* This file is part of FFmpeg.

7 
*

8 
* FFmpeg is free software; you can redistribute it and/or

9 
* modify it under the terms of the GNU Lesser General Public

10 
* License as published by the Free Software Foundation; either

11 
* version 2.1 of the License, or (at your option) any later version.

12 
*

13 
* FFmpeg is distributed in the hope that it will be useful,

14 
* but WITHOUT ANY WARRANTY; without even the implied warranty of

15 
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

16 
* Lesser General Public License for more details.

17 
*

18 
* You should have received a copy of the GNU Lesser General Public

19 
* License along with FFmpeg; if not, write to the Free Software

20 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

21 
*/

22  
23 
/**

24 
* @file libavutil/lls.c

25 
* linear least squares model

26 
*/

27  
28 
#include <math.h> 
29 
#include <string.h> 
30  
31 
#include "lls.h" 
32  
33 
void av_init_lls(LLSModel *m, int indep_count){ 
34 
memset(m, 0, sizeof(LLSModel)); 
35  
36 
m>indep_count= indep_count; 
37 
} 
38  
39 
void av_update_lls(LLSModel *m, double *var, double decay){ 
40 
int i,j;

41  
42 
for(i=0; i<=m>indep_count; i++){ 
43 
for(j=i; j<=m>indep_count; j++){

44 
m>covariance[i][j] *= decay; 
45 
m>covariance[i][j] += var[i]*var[j]; 
46 
} 
47 
} 
48 
} 
49  
50 
void av_solve_lls(LLSModel *m, double threshold, int min_order){ 
51 
int i,j,k;

52 
double (*factor)[MAX_VARS+1]= (void*)&m>covariance[1][0]; 
53 
double (*covar )[MAX_VARS+1]= (void*)&m>covariance[1][1]; 
54 
double *covar_y = m>covariance[0]; 
55 
int count= m>indep_count;

56  
57 
for(i=0; i<count; i++){ 
58 
for(j=i; j<count; j++){

59 
double sum= covar[i][j];

60  
61 
for(k=i1; k>=0; k) 
62 
sum = factor[i][k]*factor[j][k]; 
63  
64 
if(i==j){

65 
if(sum < threshold)

66 
sum= 1.0; 
67 
factor[i][i]= sqrt(sum); 
68 
}else

69 
factor[j][i]= sum / factor[i][i]; 
70 
} 
71 
} 
72 
for(i=0; i<count; i++){ 
73 
double sum= covar_y[i+1]; 
74 
for(k=i1; k>=0; k) 
75 
sum = factor[i][k]*m>coeff[0][k];

76 
m>coeff[0][i]= sum / factor[i][i];

77 
} 
78  
79 
for(j=count1; j>=min_order; j){ 
80 
for(i=j; i>=0; i){ 
81 
double sum= m>coeff[0][i]; 
82 
for(k=i+1; k<=j; k++) 
83 
sum = factor[k][i]*m>coeff[j][k]; 
84 
m>coeff[j][i]= sum / factor[i][i]; 
85 
} 
86  
87 
m>variance[j]= covar_y[0];

88 
for(i=0; i<=j; i++){ 
89 
double sum= m>coeff[j][i]*covar[i][i]  2*covar_y[i+1]; 
90 
for(k=0; k<i; k++) 
91 
sum += 2*m>coeff[j][k]*covar[k][i];

92 
m>variance[j] += m>coeff[j][i]*sum; 
93 
} 
94 
} 
95 
} 
96  
97 
double av_evaluate_lls(LLSModel *m, double *param, int order){ 
98 
int i;

99 
double out= 0; 
100  
101 
for(i=0; i<=order; i++) 
102 
out+= param[i]*m>coeff[order][i]; 
103  
104 
return out;

105 
} 
106  
107 
#ifdef TEST

108  
109 
#include <stdlib.h> 
110 
#include <stdio.h> 
111  
112 
int main(void){ 
113 
LLSModel m; 
114 
int i, order;

115  
116 
av_init_lls(&m, 3);

117  
118 
for(i=0; i<100; i++){ 
119 
double var[4]; 
120 
double eval;

121 
#if 0

122 
var[1] = rand() / (double)RAND_MAX;

123 
var[2] = rand() / (double)RAND_MAX;

124 
var[3] = rand() / (double)RAND_MAX;

125 

126 
var[2]= var[1] + var[3]/2;

127 

128 
var[0] = var[1] + var[2] + var[3] + var[1]*var[2]/100;

129 
#else

130 
var[0] = (rand() / (double)RAND_MAX  0.5)*2; 
131 
var[1] = var[0] + rand() / (double)RAND_MAX  0.5; 
132 
var[2] = var[1] + rand() / (double)RAND_MAX  0.5; 
133 
var[3] = var[2] + rand() / (double)RAND_MAX  0.5; 
134 
#endif

135 
av_update_lls(&m, var, 0.99); 
136 
av_solve_lls(&m, 0.001, 0); 
137 
for(order=0; order<3; order++){ 
138 
eval= av_evaluate_lls(&m, var+1, order);

139 
printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n",

140 
var[0], order, eval, sqrt(m.variance[order] / (i+1)), 
141 
m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]); 
142 
} 
143 
} 
144 
return 0; 
145 
} 
146  
147 
#endif
