Libav
|
00001 /* 00002 * linear least squares model 00003 * 00004 * Copyright (c) 2006 Michael Niedermayer <michaelni@gmx.at> 00005 * 00006 * This file is part of FFmpeg. 00007 * 00008 * FFmpeg is free software; you can redistribute it and/or 00009 * modify it under the terms of the GNU Lesser General Public 00010 * License as published by the Free Software Foundation; either 00011 * version 2.1 of the License, or (at your option) any later version. 00012 * 00013 * FFmpeg is distributed in the hope that it will be useful, 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00016 * Lesser General Public License for more details. 00017 * 00018 * You should have received a copy of the GNU Lesser General Public 00019 * License along with FFmpeg; if not, write to the Free Software 00020 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00021 */ 00022 00028 #include <math.h> 00029 #include <string.h> 00030 00031 #include "lls.h" 00032 00033 void av_init_lls(LLSModel *m, int indep_count){ 00034 memset(m, 0, sizeof(LLSModel)); 00035 00036 m->indep_count= indep_count; 00037 } 00038 00039 void av_update_lls(LLSModel *m, double *var, double decay){ 00040 int i,j; 00041 00042 for(i=0; i<=m->indep_count; i++){ 00043 for(j=i; j<=m->indep_count; j++){ 00044 m->covariance[i][j] *= decay; 00045 m->covariance[i][j] += var[i]*var[j]; 00046 } 00047 } 00048 } 00049 00050 void av_solve_lls(LLSModel *m, double threshold, int min_order){ 00051 int i,j,k; 00052 double (*factor)[MAX_VARS+1]= (void*)&m->covariance[1][0]; 00053 double (*covar )[MAX_VARS+1]= (void*)&m->covariance[1][1]; 00054 double *covar_y = m->covariance[0]; 00055 int count= m->indep_count; 00056 00057 for(i=0; i<count; i++){ 00058 for(j=i; j<count; j++){ 00059 double sum= covar[i][j]; 00060 00061 for(k=i-1; k>=0; k--) 00062 sum -= factor[i][k]*factor[j][k]; 00063 00064 if(i==j){ 00065 if(sum < threshold) 00066 sum= 1.0; 00067 factor[i][i]= sqrt(sum); 00068 }else 00069 factor[j][i]= sum / factor[i][i]; 00070 } 00071 } 00072 for(i=0; i<count; i++){ 00073 double sum= covar_y[i+1]; 00074 for(k=i-1; k>=0; k--) 00075 sum -= factor[i][k]*m->coeff[0][k]; 00076 m->coeff[0][i]= sum / factor[i][i]; 00077 } 00078 00079 for(j=count-1; j>=min_order; j--){ 00080 for(i=j; i>=0; i--){ 00081 double sum= m->coeff[0][i]; 00082 for(k=i+1; k<=j; k++) 00083 sum -= factor[k][i]*m->coeff[j][k]; 00084 m->coeff[j][i]= sum / factor[i][i]; 00085 } 00086 00087 m->variance[j]= covar_y[0]; 00088 for(i=0; i<=j; i++){ 00089 double sum= m->coeff[j][i]*covar[i][i] - 2*covar_y[i+1]; 00090 for(k=0; k<i; k++) 00091 sum += 2*m->coeff[j][k]*covar[k][i]; 00092 m->variance[j] += m->coeff[j][i]*sum; 00093 } 00094 } 00095 } 00096 00097 double av_evaluate_lls(LLSModel *m, double *param, int order){ 00098 int i; 00099 double out= 0; 00100 00101 for(i=0; i<=order; i++) 00102 out+= param[i]*m->coeff[order][i]; 00103 00104 return out; 00105 } 00106 00107 #ifdef TEST 00108 00109 #include <stdlib.h> 00110 #include <stdio.h> 00111 00112 int main(void){ 00113 LLSModel m; 00114 int i, order; 00115 00116 av_init_lls(&m, 3); 00117 00118 for(i=0; i<100; i++){ 00119 double var[4]; 00120 double eval; 00121 var[0] = (rand() / (double)RAND_MAX - 0.5)*2; 00122 var[1] = var[0] + rand() / (double)RAND_MAX - 0.5; 00123 var[2] = var[1] + rand() / (double)RAND_MAX - 0.5; 00124 var[3] = var[2] + rand() / (double)RAND_MAX - 0.5; 00125 av_update_lls(&m, var, 0.99); 00126 av_solve_lls(&m, 0.001, 0); 00127 for(order=0; order<3; order++){ 00128 eval= av_evaluate_lls(&m, var+1, order); 00129 printf("real:%9f order:%d pred:%9f var:%f coeffs:%f %9f %9f\n", 00130 var[0], order, eval, sqrt(m.variance[order] / (i+1)), 00131 m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]); 00132 } 00133 } 00134 return 0; 00135 } 00136 00137 #endif