lls.c 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137
  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. * @file
  24. * linear least squares model
  25. */
  26. #include <math.h>
  27. #include <string.h>
  28. #include "lls.h"
  29. void av_init_lls(LLSModel *m, int indep_count){
  30. memset(m, 0, sizeof(LLSModel));
  31. m->indep_count= indep_count;
  32. }
  33. void av_update_lls(LLSModel *m, double *var, double decay){
  34. int i,j;
  35. for(i=0; i<=m->indep_count; i++){
  36. for(j=i; j<=m->indep_count; j++){
  37. m->covariance[i][j] *= decay;
  38. m->covariance[i][j] += var[i]*var[j];
  39. }
  40. }
  41. }
  42. void av_solve_lls(LLSModel *m, double threshold, int min_order){
  43. int i,j,k;
  44. double (*factor)[MAX_VARS+1]= (void*)&m->covariance[1][0];
  45. double (*covar )[MAX_VARS+1]= (void*)&m->covariance[1][1];
  46. double *covar_y = m->covariance[0];
  47. int count= m->indep_count;
  48. for(i=0; i<count; i++){
  49. for(j=i; j<count; j++){
  50. double sum= covar[i][j];
  51. for(k=i-1; k>=0; k--)
  52. sum -= factor[i][k]*factor[j][k];
  53. if(i==j){
  54. if(sum < threshold)
  55. sum= 1.0;
  56. factor[i][i]= sqrt(sum);
  57. }else
  58. factor[j][i]= sum / factor[i][i];
  59. }
  60. }
  61. for(i=0; i<count; i++){
  62. double sum= covar_y[i+1];
  63. for(k=i-1; k>=0; k--)
  64. sum -= factor[i][k]*m->coeff[0][k];
  65. m->coeff[0][i]= sum / factor[i][i];
  66. }
  67. for(j=count-1; j>=min_order; j--){
  68. for(i=j; i>=0; i--){
  69. double sum= m->coeff[0][i];
  70. for(k=i+1; k<=j; k++)
  71. sum -= factor[k][i]*m->coeff[j][k];
  72. m->coeff[j][i]= sum / factor[i][i];
  73. }
  74. m->variance[j]= covar_y[0];
  75. for(i=0; i<=j; i++){
  76. double sum= m->coeff[j][i]*covar[i][i] - 2*covar_y[i+1];
  77. for(k=0; k<i; k++)
  78. sum += 2*m->coeff[j][k]*covar[k][i];
  79. m->variance[j] += m->coeff[j][i]*sum;
  80. }
  81. }
  82. }
  83. double av_evaluate_lls(LLSModel *m, double *param, int order){
  84. int i;
  85. double out= 0;
  86. for(i=0; i<=order; i++)
  87. out+= param[i]*m->coeff[order][i];
  88. return out;
  89. }
  90. #ifdef TEST
  91. #include <stdlib.h>
  92. #include <stdio.h>
  93. int main(void){
  94. LLSModel m;
  95. int i, order;
  96. av_init_lls(&m, 3);
  97. for(i=0; i<100; i++){
  98. double var[4];
  99. double eval;
  100. var[0] = (rand() / (double)RAND_MAX - 0.5)*2;
  101. var[1] = var[0] + rand() / (double)RAND_MAX - 0.5;
  102. var[2] = var[1] + rand() / (double)RAND_MAX - 0.5;
  103. var[3] = var[2] + rand() / (double)RAND_MAX - 0.5;
  104. av_update_lls(&m, var, 0.99);
  105. av_solve_lls(&m, 0.001, 0);
  106. for(order=0; order<3; order++){
  107. eval= av_evaluate_lls(&m, var+1, order);
  108. printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n",
  109. var[0], order, eval, sqrt(m.variance[order] / (i+1)),
  110. m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]);
  111. }
  112. }
  113. return 0;
  114. }
  115. #endif