• Main Page
  • Related Pages
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

libavutil/lls.c

Go to the documentation of this file.
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

Generated on Fri Sep 16 2011 17:17:51 for FFmpeg by  doxygen 1.7.1