blob: 0560b6a79caf76ded3c696d90403da8c84555332
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 "internal.h" |
33 | #include "version.h" |
34 | #include "lls.h" |
35 | |
36 | static void update_lls(LLSModel *m, const double *var) |
37 | { |
38 | int i, j; |
39 | |
40 | for (i = 0; i <= m->indep_count; i++) { |
41 | for (j = i; j <= m->indep_count; j++) { |
42 | m->covariance[i][j] += var[i] * var[j]; |
43 | } |
44 | } |
45 | } |
46 | |
47 | void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order) |
48 | { |
49 | int i, j, k; |
50 | double (*factor)[MAX_VARS_ALIGN] = (void *) &m->covariance[1][0]; |
51 | double (*covar) [MAX_VARS_ALIGN] = (void *) &m->covariance[1][1]; |
52 | double *covar_y = m->covariance[0]; |
53 | int count = m->indep_count; |
54 | |
55 | for (i = 0; i < count; i++) { |
56 | for (j = i; j < count; j++) { |
57 | double sum = covar[i][j]; |
58 | |
59 | for (k = 0; k <= i-1; k++) |
60 | sum -= factor[i][k] * factor[j][k]; |
61 | |
62 | if (i == j) { |
63 | if (sum < threshold) |
64 | sum = 1.0; |
65 | factor[i][i] = sqrt(sum); |
66 | } else { |
67 | factor[j][i] = sum / factor[i][i]; |
68 | } |
69 | } |
70 | } |
71 | |
72 | for (i = 0; i < count; i++) { |
73 | double sum = covar_y[i + 1]; |
74 | |
75 | for (k = 0; k <= i-1; k++) |
76 | sum -= factor[i][k] * m->coeff[0][k]; |
77 | |
78 | m->coeff[0][i] = sum / factor[i][i]; |
79 | } |
80 | |
81 | for (j = count - 1; j >= min_order; j--) { |
82 | for (i = j; i >= 0; i--) { |
83 | double sum = m->coeff[0][i]; |
84 | |
85 | for (k = i + 1; k <= j; k++) |
86 | sum -= factor[k][i] * m->coeff[j][k]; |
87 | |
88 | m->coeff[j][i] = sum / factor[i][i]; |
89 | } |
90 | |
91 | m->variance[j] = covar_y[0]; |
92 | |
93 | for (i = 0; i <= j; i++) { |
94 | double sum = m->coeff[j][i] * covar[i][i] - 2 * covar_y[i + 1]; |
95 | |
96 | for (k = 0; k < i; k++) |
97 | sum += 2 * m->coeff[j][k] * covar[k][i]; |
98 | |
99 | m->variance[j] += m->coeff[j][i] * sum; |
100 | } |
101 | } |
102 | } |
103 | |
104 | static double evaluate_lls(LLSModel *m, const double *param, int order) |
105 | { |
106 | int i; |
107 | double out = 0; |
108 | |
109 | for (i = 0; i <= order; i++) |
110 | out += param[i] * m->coeff[order][i]; |
111 | |
112 | return out; |
113 | } |
114 | |
115 | av_cold void avpriv_init_lls(LLSModel *m, int indep_count) |
116 | { |
117 | memset(m, 0, sizeof(LLSModel)); |
118 | m->indep_count = indep_count; |
119 | m->update_lls = update_lls; |
120 | m->evaluate_lls = evaluate_lls; |
121 | if (ARCH_X86) |
122 | ff_init_lls_x86(m); |
123 | } |
124 |