1# fit a polynomial to a set of points 2# fit -dn [-v] 3# where n is the degree of the polynomial 4 5implement Fit; 6 7include "sys.m"; 8 sys: Sys; 9include "draw.m"; 10include "math.m"; 11 maths: Math; 12include "bufio.m"; 13 bufio: Bufio; 14include "arg.m"; 15 16Fit: module 17{ 18 init: fn(nil: ref Draw->Context, argv: list of string); 19}; 20 21MAXPTS: con 512; 22MAXDEG: con 16; 23EPS: con 0.0000005; 24 25init(nil: ref Draw->Context, argv: list of string) 26{ 27 sys = load Sys Sys->PATH; 28 maths = load Math Math->PATH; 29 if(maths == nil) 30 fatal(sys->sprint("cannot load maths library")); 31 bufio = load Bufio Bufio->PATH; 32 if(bufio == nil) 33 fatal(sys->sprint("cannot load bufio")); 34 main(argv); 35} 36 37isn(r: real, n: int): int 38{ 39 s := r - real n; 40 if(s < 0.0) 41 s = -s; 42 return s < EPS; 43} 44 45fact(n: int): real 46{ 47 f := 1.0; 48 for(i := 1; i <= n; i++) 49 f *= real i; 50 return f; 51} 52 53comb(n: int, r: int): real 54{ 55 f := 1.0; 56 for(i := 0; i < r; i++) 57 f *= real (n-i); 58 return f/fact(r); 59} 60 61matalloc(n: int): array of array of real 62{ 63 mat := array[n] of array of real; 64 for(i := 0; i < n; i++) 65 mat[i] = array[n] of real; 66 return mat; 67} 68 69matsalloc(n: int): array of array of array of real 70{ 71 mats := array[n+1] of array of array of real; 72 for(i := 0; i <= n; i++) 73 mats[i] = matalloc(i); 74 return mats; 75} 76 77det(mat: array of array of real, n: int, mats: array of array of array of real): real 78{ 79 # easy cases first 80 if(n == 0) 81 return 1.0; 82 if(n == 1) 83 return mat[0][0]; 84 if(n == 2) 85 return mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0]; 86 d := 0.0; 87 s := 1; 88 m := mats[n-1]; 89 for(k := 0; k < n; k++){ 90 for(i := 0; i < n-1; i++){ 91 for(j := 0; j < n-1; j++){ 92 if(j < k) 93 m[i][j] = mat[i+1][j]; 94 else 95 m[i][j] = mat[i+1][j+1]; 96 } 97 } 98 d += (real s)*mat[0][k]*det(m, n-1, mats); 99 s = -s; 100 } 101 return d; 102} 103 104main(argv: list of string) 105{ 106 i, j: int; 107 x, y, z: real; 108 fb: ref Bufio->Iobuf; 109 110 n := 0; 111 p := 1; 112 arg := load Arg Arg->PATH; 113 if(arg == nil) 114 fatal(sys->sprint("cannot load %s: %r", Arg->PATH)); 115 arg->init(argv); 116 verbose := 0; 117 while((o := arg->opt()) != 0) 118 case o{ 119 'd' => 120 p = int arg->arg(); 121 'v' => 122 verbose = 1; 123 * => 124 fatal(sys->sprint("bad option %c", o)); 125 } 126 args := arg->argv(); 127 arg = nil; 128 if(args != nil){ 129 s := hd args; 130 fb = bufio->open(s, bufio->OREAD); 131 if(fb == nil) 132 fatal(sys->sprint("cannot open %s", s)); 133 } 134 else{ 135 fb = bufio->open("/dev/cons", bufio->OREAD); 136 if(fb == nil) 137 fatal(sys->sprint("missing data file name")); 138 } 139 a := array[p+1] of real; 140 b := array[p+1] of real; 141 sx := array[2*p+1] of real; 142 sxy := array[p+1] of real; 143 xd := array[MAXPTS] of real; 144 yd := array[MAXPTS] of real; 145 while(1){ 146 xs := ss(bufio->fb.gett(" \t\r\n")); 147 if(xs == nil) 148 break; 149 ys := ss(bufio->fb.gett(" \t\r\n")); 150 if(ys == nil) 151 fatal(sys->sprint("missing value")); 152 if(n >= MAXPTS) 153 fatal(sys->sprint("too many points")); 154 xd[n] = real xs; 155 yd[n] = real ys; 156 n++; 157 } 158 if(p < 0) 159 fatal(sys->sprint("negative power")); 160 if(p > MAXDEG) 161 fatal(sys->sprint("power too large")); 162 if(n < p+1) 163 fatal(sys->sprint("not enough points")); 164 # use x-xbar, y-ybar to avoid overflow 165 for(i = 0; i <= p; i++) 166 sxy[i] = 0.0; 167 for(i = 0; i <= 2*p; i++) 168 sx[i] = 0.0; 169 xbar := ybar := 0.0; 170 for(i = 0; i < n; i++){ 171 xbar += xd[i]; 172 ybar += yd[i]; 173 } 174 xbar = xbar/(real n); 175 ybar = ybar/(real n); 176 for(i = 0; i < n; i++){ 177 x = xd[i]-xbar; 178 y = yd[i]-ybar; 179 for(j = 0; j <= p; j++) 180 sxy[j] += y*x**j; 181 for(j = 0; j <= 2*p; j++) 182 sx[j] += x**j; 183 } 184 mats := matsalloc(p+1); 185 mat := mats[p+1]; 186 for(i = 0; i <= p; i++) 187 for(j = 0; j <= p; j++) 188 mat[i][j] = sx[i+j]; 189 d := det(mat, p+1, mats); 190 if(isn(d, 0)) 191 fatal(sys->sprint("points not independent")); 192 for(j = 0; j <= p; j++){ 193 for(i = 0; i <= p; i++) 194 mat[i][j] = sxy[i]; 195 a[j] = det(mat, p+1, mats)/d; 196 for(i = 0; i <= p; i++) 197 mat[i][j] = sx[i+j]; 198 } 199 if(verbose) 200 sys->print("\npt actual x actual y predicted y\n"); 201 e := 0.0; 202 for(i = 0; i < n; i++){ 203 x = xd[i]-xbar; 204 y = yd[i]-ybar; 205 z = 0.0; 206 for(j = 0; j <= p; j++) 207 z += a[j]*x**j; 208 z += ybar; 209 e += (z-yd[i])*(z-yd[i]); 210 if(verbose) 211 sys->print("%d. %f %f %f\n", i+1, xd[i], yd[i], z); 212 } 213 if(verbose) 214 sys->print("root mean squared error = %f\n", maths->sqrt(e/(real n))); 215 for(i = 0; i <= p; i++) 216 b[i] = 0.0; 217 b[0] += ybar; 218 for(i = 0; i <= p; i++) 219 for(j = 0; j <= i; j++) 220 b[j] += a[i]*comb(i, j)*(-xbar)**(i-j); 221 pr := 0; 222 sys->print("y = "); 223 for(i = p; i >= 0; i--){ 224 if(!isn(b[i], 0) || (i == 0 && pr == 0)){ 225 if(b[i] < 0.0){ 226 sys->print("-"); 227 b[i] = -b[i]; 228 } 229 else if(pr) 230 sys->print("+"); 231 pr = 1; 232 if(i == 0) 233 sys->print("%f", b[i]); 234 else{ 235 if(!isn(b[i], 1)) 236 sys->print("%f*", b[i]); 237 sys->print("x"); 238 if(i > 1) 239 sys->print("^%d", i); 240 } 241 } 242 } 243 sys->print("\n"); 244} 245 246ss(s: string): string 247{ 248 l := len s; 249 while(l > 0 && (s[0] == ' ' || s[0] == '\t' || s[0] == '\r' || s[0] == '\n')){ 250 s = s[1: ]; 251 l--; 252 } 253 while(l > 0 && (s[l-1] == ' ' || s[l-1] == '\t' || s[l-1] == '\r' || s[l-1] == '\n')){ 254 s = s[0: l-1]; 255 l--; 256 } 257 return s; 258} 259 260fatal(s: string) 261{ 262 sys->fprint(sys->fildes(2), "fit: %s\n", s); 263 exit; 264} 265