00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00028 #include <math.h>
00029 #include <string.h>
00030
00031 #include "lls.h"
00032
00033 #ifdef TEST
00034 #define av_log(a,b,...) printf(__VA_ARGS__)
00035 #endif
00036
00037 void av_init_lls(LLSModel *m, int indep_count){
00038 memset(m, 0, sizeof(LLSModel));
00039
00040 m->indep_count= indep_count;
00041 }
00042
00043 void av_update_lls(LLSModel *m, double *var, double decay){
00044 int i,j;
00045
00046 for(i=0; i<=m->indep_count; i++){
00047 for(j=i; j<=m->indep_count; j++){
00048 m->covariance[i][j] *= decay;
00049 m->covariance[i][j] += var[i]*var[j];
00050 }
00051 }
00052 }
00053
00054 void av_solve_lls(LLSModel *m, double threshold, int min_order){
00055 int i,j,k;
00056 double (*factor)[MAX_VARS+1]= (void*)&m->covariance[1][0];
00057 double (*covar )[MAX_VARS+1]= (void*)&m->covariance[1][1];
00058 double *covar_y = m->covariance[0];
00059 int count= m->indep_count;
00060
00061 for(i=0; i<count; i++){
00062 for(j=i; j<count; j++){
00063 double sum= covar[i][j];
00064
00065 for(k=i-1; k>=0; k--)
00066 sum -= factor[i][k]*factor[j][k];
00067
00068 if(i==j){
00069 if(sum < threshold)
00070 sum= 1.0;
00071 factor[i][i]= sqrt(sum);
00072 }else
00073 factor[j][i]= sum / factor[i][i];
00074 }
00075 }
00076 for(i=0; i<count; i++){
00077 double sum= covar_y[i+1];
00078 for(k=i-1; k>=0; k--)
00079 sum -= factor[i][k]*m->coeff[0][k];
00080 m->coeff[0][i]= sum / factor[i][i];
00081 }
00082
00083 for(j=count-1; j>=min_order; j--){
00084 for(i=j; i>=0; i--){
00085 double sum= m->coeff[0][i];
00086 for(k=i+1; k<=j; k++)
00087 sum -= factor[k][i]*m->coeff[j][k];
00088 m->coeff[j][i]= sum / factor[i][i];
00089 }
00090
00091 m->variance[j]= covar_y[0];
00092 for(i=0; i<=j; i++){
00093 double sum= m->coeff[j][i]*covar[i][i] - 2*covar_y[i+1];
00094 for(k=0; k<i; k++)
00095 sum += 2*m->coeff[j][k]*covar[k][i];
00096 m->variance[j] += m->coeff[j][i]*sum;
00097 }
00098 }
00099 }
00100
00101 double av_evaluate_lls(LLSModel *m, double *param, int order){
00102 int i;
00103 double out= 0;
00104
00105 for(i=0; i<=order; i++)
00106 out+= param[i]*m->coeff[order][i];
00107
00108 return out;
00109 }
00110
00111 #ifdef TEST
00112
00113 #include <stdlib.h>
00114 #include <stdio.h>
00115
00116 int main(void){
00117 LLSModel m;
00118 int i, order;
00119
00120 av_init_lls(&m, 3);
00121
00122 for(i=0; i<100; i++){
00123 double var[4];
00124 double eval;
00125 #if 0
00126 var[1] = rand() / (double)RAND_MAX;
00127 var[2] = rand() / (double)RAND_MAX;
00128 var[3] = rand() / (double)RAND_MAX;
00129
00130 var[2]= var[1] + var[3]/2;
00131
00132 var[0] = var[1] + var[2] + var[3] + var[1]*var[2]/100;
00133 #else
00134 var[0] = (rand() / (double)RAND_MAX - 0.5)*2;
00135 var[1] = var[0] + rand() / (double)RAND_MAX - 0.5;
00136 var[2] = var[1] + rand() / (double)RAND_MAX - 0.5;
00137 var[3] = var[2] + rand() / (double)RAND_MAX - 0.5;
00138 #endif
00139 av_update_lls(&m, var, 0.99);
00140 av_solve_lls(&m, 0.001, 0);
00141 for(order=0; order<3; order++){
00142 eval= av_evaluate_lls(&m, var+1, order);
00143 av_log(NULL, AV_LOG_DEBUG, "real:%f order:%d pred:%f var:%f coeffs:%f %f %f\n",
00144 var[0], order, eval, sqrt(m.variance[order] / (i+1)),
00145 m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]);
00146 }
00147 }
00148 return 0;
00149 }
00150
00151 #endif