diff options
Diffstat (limited to 'appl/math/fit.b')
| -rw-r--r-- | appl/math/fit.b | 264 |
1 files changed, 264 insertions, 0 deletions
diff --git a/appl/math/fit.b b/appl/math/fit.b new file mode 100644 index 00000000..b512c13b --- /dev/null +++ b/appl/math/fit.b @@ -0,0 +1,264 @@ +# fit a polynomial to a set of points +# fit -dn [-v] +# where n is the degree of the polynomial + +implement Fit; + +include "sys.m"; + sys: Sys; +include "draw.m"; +include "math.m"; + maths: Math; +include "bufio.m"; + bufio: Bufio; +include "arg.m"; + +Fit: module +{ + init: fn(nil: ref Draw->Context, argv: list of string); +}; + +MAXPTS: con 512; +MAXDEG: con 16; +EPS: con 0.0000005; + +init(nil: ref Draw->Context, argv: list of string) +{ + sys = load Sys Sys->PATH; + maths = load Math Math->PATH; + if(maths == nil) + fatal(sys->sprint("cannot load maths library")); + bufio = load Bufio Bufio->PATH; + if(bufio == nil) + fatal(sys->sprint("cannot load bufio")); + main(argv); +} + +isn(r: real, n: int): int +{ + s := r - real n; + if(s < 0.0) + s = -s; + return s < EPS; +} + +fact(n: int): real +{ + f := 1.0; + for(i := 1; i <= n; i++) + f *= real i; + return f; +} + +comb(n: int, r: int): real +{ + f := 1.0; + for(i := 0; i < r; i++) + f *= real (n-i); + return f/fact(r); +} + +matalloc(n: int): array of array of real +{ + mat := array[n] of array of real; + for(i := 0; i < n; i++) + mat[i] = array[n] of real; + return mat; +} + +matsalloc(n: int): array of array of array of real +{ + mats := array[n+1] of array of array of real; + for(i := 0; i <= n; i++) + mats[i] = matalloc(i); + return mats; +} + +det(mat: array of array of real, n: int, mats: array of array of array of real): real +{ + # easy cases first + if(n == 0) + return 1.0; + if(n == 1) + return mat[0][0]; + if(n == 2) + return mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0]; + d := 0.0; + s := 1; + m := mats[n-1]; + for(k := 0; k < n; k++){ + for(i := 0; i < n-1; i++){ + for(j := 0; j < n-1; j++){ + if(j < k) + m[i][j] = mat[i+1][j]; + else + m[i][j] = mat[i+1][j+1]; + } + } + d += (real s)*mat[0][k]*det(m, n-1, mats); + s = -s; + } + return d; +} + +main(argv: list of string) +{ + i, j: int; + x, y, z: real; + fb: ref Bufio->Iobuf; + + n := 0; + p := 1; + arg := load Arg Arg->PATH; + if(arg == nil) + fatal(sys->sprint("cannot load %s: %r", Arg->PATH)); + arg->init(argv); + verbose := 0; + while((o := arg->opt()) != 0) + case o{ + 'd' => + p = int arg->arg(); + 'v' => + verbose = 1; + * => + fatal(sys->sprint("bad option %c", o)); + } + args := arg->argv(); + arg = nil; + if(args != nil){ + s := hd args; + fb = bufio->open(s, bufio->OREAD); + if(fb == nil) + fatal(sys->sprint("cannot open %s", s)); + } + else{ + fb = bufio->open("/dev/cons", bufio->OREAD); + if(fb == nil) + fatal(sys->sprint("missing data file name")); + } + a := array[p+1] of real; + b := array[p+1] of real; + sx := array[2*p+1] of real; + sxy := array[p+1] of real; + xd := array[MAXPTS] of real; + yd := array[MAXPTS] of real; + while(1){ + xs := ss(bufio->fb.gett(" \t\r\n")); + if(xs == nil) + break; + ys := ss(bufio->fb.gett(" \t\r\n")); + if(ys == nil) + fatal(sys->sprint("missing value")); + if(n >= MAXPTS) + fatal(sys->sprint("too many points")); + xd[n] = real xs; + yd[n] = real ys; + n++; + } + if(p < 0) + fatal(sys->sprint("negative power")); + if(p > MAXDEG) + fatal(sys->sprint("power too large")); + if(n < p+1) + fatal(sys->sprint("not enough points")); + # use x-xbar, y-ybar to avoid overflow + for(i = 0; i <= p; i++) + sxy[i] = 0.0; + for(i = 0; i <= 2*p; i++) + sx[i] = 0.0; + xbar := ybar := 0.0; + for(i = 0; i < n; i++){ + xbar += xd[i]; + ybar += yd[i]; + } + xbar = xbar/(real n); + ybar = ybar/(real n); + for(i = 0; i < n; i++){ + x = xd[i]-xbar; + y = yd[i]-ybar; + for(j = 0; j <= p; j++) + sxy[j] += y*x**j; + for(j = 0; j <= 2*p; j++) + sx[j] += x**j; + } + mats := matsalloc(p+1); + mat := mats[p+1]; + for(i = 0; i <= p; i++) + for(j = 0; j <= p; j++) + mat[i][j] = sx[i+j]; + d := det(mat, p+1, mats); + if(isn(d, 0)) + fatal(sys->sprint("points not independent")); + for(j = 0; j <= p; j++){ + for(i = 0; i <= p; i++) + mat[i][j] = sxy[i]; + a[j] = det(mat, p+1, mats)/d; + for(i = 0; i <= p; i++) + mat[i][j] = sx[i+j]; + } + if(verbose) + sys->print("\npt actual x actual y predicted y\n"); + e := 0.0; + for(i = 0; i < n; i++){ + x = xd[i]-xbar; + y = yd[i]-ybar; + z = 0.0; + for(j = 0; j <= p; j++) + z += a[j]*x**j; + z += ybar; + e += (z-yd[i])*(z-yd[i]); + if(verbose) + sys->print("%d. %f %f %f\n", i+1, xd[i], yd[i], z); + } + if(verbose) + sys->print("root mean squared error = %f\n", maths->sqrt(e/(real n))); + for(i = 0; i <= p; i++) + b[i] = 0.0; + b[0] += ybar; + for(i = 0; i <= p; i++) + for(j = 0; j <= i; j++) + b[j] += a[i]*comb(i, j)*(-xbar)**(i-j); + pr := 0; + sys->print("y = "); + for(i = p; i >= 0; i--){ + if(!isn(b[i], 0) || (i == 0 && pr == 0)){ + if(b[i] < 0.0){ + sys->print("-"); + b[i] = -b[i]; + } + else if(pr) + sys->print("+"); + pr = 1; + if(i == 0) + sys->print("%f", b[i]); + else{ + if(!isn(b[i], 1)) + sys->print("%f*", b[i]); + sys->print("x"); + if(i > 1) + sys->print("^%d", i); + } + } + } + sys->print("\n"); +} + +ss(s: string): string +{ + l := len s; + while(l > 0 && (s[0] == ' ' || s[0] == '\t' || s[0] == '\r' || s[0] == '\n')){ + s = s[1: ]; + l--; + } + while(l > 0 && (s[l-1] == ' ' || s[l-1] == '\t' || s[l-1] == '\r' || s[l-1] == '\n')){ + s = s[0: l-1]; + l--; + } + return s; +} + +fatal(s: string) +{ + sys->fprint(sys->fildes(2), "fit: %s\n", s); + exit; +} |