1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
|
/*
* Gnome Wave Cleaner Version 0.19
* Copyright (C) 2001 Jeffrey J. Welty
FILE stat.c
PURPOSE - linear regression
CONTENTS -
*/
#include <stdio.h>
#include "stat.h"
void matrix_solve(MAT *) ;
#define abs(x) ((x) > 0 ? (x) : -(x))
#define tiny 1e-200
static MAT *coef = MNULL ;
static VEC *b = VNULL ;
static VEC *answer = VNULL ;
static int row, col, N, i ;
static int failed ;
/* LUsolve -- given an LU factorisation in A, solve Ax=b */
VEC *myLUsolve(A,pivot,b,x)
MAT *A;
PERM *pivot;
VEC *b,*x;
{
if ( A==(MAT *)NULL || b==(VEC *)NULL || pivot==(PERM *)NULL )
error(E_NULL,"LUsolve");
if ( A->m != A->n || A->n != b->dim )
error(E_SIZES,"LUsolve");
x = v_resize(x,b->dim);
px_vec(pivot,b,x); /* x := P.b */
Lsolve(A,x,x,1.0); /* implicit diagonal = 1 */
Usolve(A,x,x,0.0); /* explicit diagonal */
return (x);
}
void matrix_solve(MAT *A)
{
PERM *pivot ;
failed = 0 ;
count_errs(0) ;
pivot = px_get(A->m) ;
catchall(
LUfactor(coef, pivot) ;
answer = v_get(b->dim);
answer = myLUsolve(A,pivot,b,answer) ;
, failed=1 );
if(failed == 1) {
answer = v_resize(answer,N) ;
v_zero(answer) ;
}
PX_FREE(pivot) ;
}
void init_reg(int n)
{
N = n ;
coef = m_resize(coef, N, N) ;
b = v_resize(b, N) ;
/** zero the coef array which will hold the sums **/
m_zero(coef) ;
v_zero(b) ;
}
/* Paul Sanders' speedup 1/12/2007 works a little better */
#define PAUL_SANDERS
#ifdef JW_NEW
void sum_reg(double x[], double y)
{
for(row = 0 ; row < N ; row++) {
for(col = 0 ; col < N ; col++)
coef->me[row][col] += x[row] * x[col] ;
b->ve[row] += y * x[row] ;
}
}
#endif
#ifdef PAUL_SANDERS
void sum_reg(double x[], double y)
{
double *pbve = b->ve ;
for(row = 0 ; row < N ; row++) {
double *coef_row = coef->me[row] ;
for(col = 0 ; col < N ; col++)
coef_row[col] += x[row] * x[col] ;
*pbve += x[row] * y ;
pbve++ ;
}
}
#endif
int estimate_reg(double *b_solution)
{
matrix_solve(coef) ;
for(row = 0 ; row < N ; row++)
b_solution[row] = answer->ve[row] ;
V_FREE(answer) ;
M_FREE(coef) ;
V_FREE(b) ;
return failed ;
}
|