Commit | Line | Data |
---|---|---|
2ba45a60 DM |
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 02110-1301 USA | |
21 | */ | |
22 | ||
23 | /** | |
24 | * @file | |
25 | * linear least squares model | |
26 | */ | |
27 | ||
28 | #include <math.h> | |
29 | #include <string.h> | |
30 | ||
31 | #include "attributes.h" | |
32 | #include "version.h" | |
33 | #include "lls.h" | |
34 | ||
f6fa7814 | 35 | static void update_lls(LLSModel *m, const double *var) |
2ba45a60 DM |
36 | { |
37 | int i, j; | |
38 | ||
39 | for (i = 0; i <= m->indep_count; i++) { | |
40 | for (j = i; j <= m->indep_count; j++) { | |
41 | m->covariance[i][j] += var[i] * var[j]; | |
42 | } | |
43 | } | |
44 | } | |
45 | ||
46 | void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order) | |
47 | { | |
48 | int i, j, k; | |
49 | double (*factor)[MAX_VARS_ALIGN] = (void *) &m->covariance[1][0]; | |
50 | double (*covar) [MAX_VARS_ALIGN] = (void *) &m->covariance[1][1]; | |
51 | double *covar_y = m->covariance[0]; | |
52 | int count = m->indep_count; | |
53 | ||
54 | for (i = 0; i < count; i++) { | |
55 | for (j = i; j < count; j++) { | |
56 | double sum = covar[i][j]; | |
57 | ||
58 | for (k = i - 1; k >= 0; k--) | |
59 | sum -= factor[i][k] * factor[j][k]; | |
60 | ||
61 | if (i == j) { | |
62 | if (sum < threshold) | |
63 | sum = 1.0; | |
64 | factor[i][i] = sqrt(sum); | |
65 | } else { | |
66 | factor[j][i] = sum / factor[i][i]; | |
67 | } | |
68 | } | |
69 | } | |
70 | ||
71 | for (i = 0; i < count; i++) { | |
72 | double sum = covar_y[i + 1]; | |
73 | ||
74 | for (k = i - 1; k >= 0; k--) | |
75 | sum -= factor[i][k] * m->coeff[0][k]; | |
76 | ||
77 | m->coeff[0][i] = sum / factor[i][i]; | |
78 | } | |
79 | ||
80 | for (j = count - 1; j >= min_order; j--) { | |
81 | for (i = j; i >= 0; i--) { | |
82 | double sum = m->coeff[0][i]; | |
83 | ||
84 | for (k = i + 1; k <= j; k++) | |
85 | sum -= factor[k][i] * m->coeff[j][k]; | |
86 | ||
87 | m->coeff[j][i] = sum / factor[i][i]; | |
88 | } | |
89 | ||
90 | m->variance[j] = covar_y[0]; | |
91 | ||
92 | for (i = 0; i <= j; i++) { | |
93 | double sum = m->coeff[j][i] * covar[i][i] - 2 * covar_y[i + 1]; | |
94 | ||
95 | for (k = 0; k < i; k++) | |
96 | sum += 2 * m->coeff[j][k] * covar[k][i]; | |
97 | ||
98 | m->variance[j] += m->coeff[j][i] * sum; | |
99 | } | |
100 | } | |
101 | } | |
102 | ||
f6fa7814 | 103 | static double evaluate_lls(LLSModel *m, const double *param, int order) |
2ba45a60 DM |
104 | { |
105 | int i; | |
106 | double out = 0; | |
107 | ||
108 | for (i = 0; i <= order; i++) | |
109 | out += param[i] * m->coeff[order][i]; | |
110 | ||
111 | return out; | |
112 | } | |
113 | ||
114 | av_cold void avpriv_init_lls(LLSModel *m, int indep_count) | |
115 | { | |
116 | memset(m, 0, sizeof(LLSModel)); | |
117 | m->indep_count = indep_count; | |
118 | m->update_lls = update_lls; | |
119 | m->evaluate_lls = evaluate_lls; | |
120 | if (ARCH_X86) | |
121 | ff_init_lls_x86(m); | |
122 | } | |
123 | ||
124 | #ifdef TEST | |
125 | ||
126 | #include <stdio.h> | |
127 | #include <limits.h> | |
128 | #include "lfg.h" | |
129 | ||
130 | int main(void) | |
131 | { | |
132 | LLSModel m; | |
133 | int i, order; | |
134 | AVLFG lfg; | |
135 | ||
136 | av_lfg_init(&lfg, 1); | |
137 | avpriv_init_lls(&m, 3); | |
138 | ||
139 | for (i = 0; i < 100; i++) { | |
140 | LOCAL_ALIGNED(32, double, var, [4]); | |
141 | double eval; | |
142 | ||
143 | var[0] = (av_lfg_get(&lfg) / (double) UINT_MAX - 0.5) * 2; | |
144 | var[1] = var[0] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; | |
145 | var[2] = var[1] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; | |
146 | var[3] = var[2] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; | |
147 | m.update_lls(&m, var); | |
148 | avpriv_solve_lls(&m, 0.001, 0); | |
149 | for (order = 0; order < 3; order++) { | |
150 | eval = m.evaluate_lls(&m, var + 1, order); | |
151 | printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n", | |
152 | var[0], order, eval, sqrt(m.variance[order] / (i + 1)), | |
153 | m.coeff[order][0], m.coeff[order][1], | |
154 | m.coeff[order][2]); | |
155 | } | |
156 | } | |
157 | return 0; | |
158 | } | |
159 | ||
160 | #endif |