yading@11: /* yading@11: * linear least squares model yading@11: * yading@11: * Copyright (c) 2006 Michael Niedermayer yading@11: * yading@11: * This file is part of FFmpeg. yading@11: * yading@11: * FFmpeg is free software; you can redistribute it and/or yading@11: * modify it under the terms of the GNU Lesser General Public yading@11: * License as published by the Free Software Foundation; either yading@11: * version 2.1 of the License, or (at your option) any later version. yading@11: * yading@11: * FFmpeg is distributed in the hope that it will be useful, yading@11: * but WITHOUT ANY WARRANTY; without even the implied warranty of yading@11: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU yading@11: * Lesser General Public License for more details. yading@11: * yading@11: * You should have received a copy of the GNU Lesser General Public yading@11: * License along with FFmpeg; if not, write to the Free Software yading@11: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA yading@11: */ yading@11: yading@11: /** yading@11: * @file yading@11: * linear least squares model yading@11: */ yading@11: yading@11: #include yading@11: #include yading@11: yading@11: #include "version.h" yading@11: #include "lls.h" yading@11: yading@11: void avpriv_init_lls(LLSModel *m, int indep_count) yading@11: { yading@11: memset(m, 0, sizeof(LLSModel)); yading@11: m->indep_count = indep_count; yading@11: } yading@11: yading@11: void avpriv_update_lls(LLSModel *m, double *var, double decay) yading@11: { yading@11: int i, j; yading@11: yading@11: for (i = 0; i <= m->indep_count; i++) { yading@11: for (j = i; j <= m->indep_count; j++) { yading@11: m->covariance[i][j] *= decay; yading@11: m->covariance[i][j] += var[i] * var[j]; yading@11: } yading@11: } yading@11: } yading@11: yading@11: void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order) yading@11: { yading@11: int i, j, k; yading@11: double (*factor)[MAX_VARS + 1] = (void *) &m->covariance[1][0]; yading@11: double (*covar) [MAX_VARS + 1] = (void *) &m->covariance[1][1]; yading@11: double *covar_y = m->covariance[0]; yading@11: int count = m->indep_count; yading@11: yading@11: for (i = 0; i < count; i++) { yading@11: for (j = i; j < count; j++) { yading@11: double sum = covar[i][j]; yading@11: yading@11: for (k = i - 1; k >= 0; k--) yading@11: sum -= factor[i][k] * factor[j][k]; yading@11: yading@11: if (i == j) { yading@11: if (sum < threshold) yading@11: sum = 1.0; yading@11: factor[i][i] = sqrt(sum); yading@11: } else { yading@11: factor[j][i] = sum / factor[i][i]; yading@11: } yading@11: } yading@11: } yading@11: yading@11: for (i = 0; i < count; i++) { yading@11: double sum = covar_y[i + 1]; yading@11: yading@11: for (k = i - 1; k >= 0; k--) yading@11: sum -= factor[i][k] * m->coeff[0][k]; yading@11: yading@11: m->coeff[0][i] = sum / factor[i][i]; yading@11: } yading@11: yading@11: for (j = count - 1; j >= min_order; j--) { yading@11: for (i = j; i >= 0; i--) { yading@11: double sum = m->coeff[0][i]; yading@11: yading@11: for (k = i + 1; k <= j; k++) yading@11: sum -= factor[k][i] * m->coeff[j][k]; yading@11: yading@11: m->coeff[j][i] = sum / factor[i][i]; yading@11: } yading@11: yading@11: m->variance[j] = covar_y[0]; yading@11: yading@11: for (i = 0; i <= j; i++) { yading@11: double sum = m->coeff[j][i] * covar[i][i] - 2 * covar_y[i + 1]; yading@11: yading@11: for (k = 0; k < i; k++) yading@11: sum += 2 * m->coeff[j][k] * covar[k][i]; yading@11: yading@11: m->variance[j] += m->coeff[j][i] * sum; yading@11: } yading@11: } yading@11: } yading@11: yading@11: double avpriv_evaluate_lls(LLSModel *m, double *param, int order) yading@11: { yading@11: int i; yading@11: double out = 0; yading@11: yading@11: for (i = 0; i <= order; i++) yading@11: out += param[i] * m->coeff[order][i]; yading@11: yading@11: return out; yading@11: } yading@11: yading@11: #if FF_API_LLS_PRIVATE yading@11: void av_init_lls(LLSModel *m, int indep_count) yading@11: { yading@11: avpriv_init_lls(m, indep_count); yading@11: } yading@11: void av_update_lls(LLSModel *m, double *param, double decay) yading@11: { yading@11: avpriv_update_lls(m, param, decay); yading@11: } yading@11: void av_solve_lls(LLSModel *m, double threshold, int min_order) yading@11: { yading@11: avpriv_solve_lls(m, threshold, min_order); yading@11: } yading@11: double av_evaluate_lls(LLSModel *m, double *param, int order) yading@11: { yading@11: return avpriv_evaluate_lls(m, param, order); yading@11: } yading@11: #endif /* FF_API_LLS_PRIVATE */ yading@11: yading@11: #ifdef TEST yading@11: yading@11: #include yading@11: #include yading@11: #include "lfg.h" yading@11: yading@11: int main(void) yading@11: { yading@11: LLSModel m; yading@11: int i, order; yading@11: AVLFG lfg; yading@11: yading@11: av_lfg_init(&lfg, 1); yading@11: avpriv_init_lls(&m, 3); yading@11: yading@11: for (i = 0; i < 100; i++) { yading@11: double var[4]; yading@11: double eval; yading@11: yading@11: var[0] = (av_lfg_get(&lfg) / (double) UINT_MAX - 0.5) * 2; yading@11: var[1] = var[0] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; yading@11: var[2] = var[1] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; yading@11: var[3] = var[2] + av_lfg_get(&lfg) / (double) UINT_MAX - 0.5; yading@11: avpriv_update_lls(&m, var, 0.99); yading@11: avpriv_solve_lls(&m, 0.001, 0); yading@11: for (order = 0; order < 3; order++) { yading@11: eval = avpriv_evaluate_lls(&m, var + 1, order); yading@11: printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n", yading@11: var[0], order, eval, sqrt(m.variance[order] / (i + 1)), yading@11: m.coeff[order][0], m.coeff[order][1], yading@11: m.coeff[order][2]); yading@11: } yading@11: } yading@11: return 0; yading@11: } yading@11: yading@11: #endif