unit LeastSquares; {Contains the curve fitting routines based on the Simplex method described in the article " Fitting Curves to Data "} {in the May 1984 issue of Byte magazine, pages 340-362.} interface uses QuickDraw, Palettes, Printing, globals, Utilities, graphics; const nvpp = 2; maxn = 6; maxnp = 20; alpha = 1.0; beta = 0.5; gamma = 2.0; lw = 5; root2 = 1.414214; MaxError = 1e-7; type ColumnVector = array[1..maxnp] of extended; vector = array[1..maxn] of extended; datarow = array[1..nvpp] of extended; index = 0..255; var m, n: integer; done: boolean; maxx, maxy: extended; i, j: index; h, l: array[1..maxn] of index; np, npmax, niter, maxiter: integer; next, center, smean, error, maxerr, p, q, step: vector; simp: array[1..maxn] of vector; data: array[1..maxnp] of datarow; filename, newname: string; yoffset: integer; procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector); implementation function f (p: vector; d: datarow): extended; var x, y, ex: extended; begin x := d[1]; case info^.fit of StraightLine: f := p[1] + p[2] * x; Poly2: f := p[1] + p[2] * x + p[3] * x * x; Poly3: f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x; Poly4: f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x + p[5] * x * x * x * x; ExpoFit: f := p[1] * exp(p[2] * x); PowerFit: if x = 0.0 then f := 0.0 else f := p[1] * exp(p[2] * ln(x)); {y=ax^b} LogFit: begin if x = 0.0 then x := 0.5; f := p[1] * ln(p[2] * x) end; RodbardFit: begin if x = 0.0 then ex := 0.0 else ex := exp(ln(x / p[3]) * p[2]); y := p[1] - p[4]; y := y / (1 + ex); f := y + p[4]; end; {Rodbard fit} end; {case} end; procedure order; var i, j: index; begin for j := 1 to n do begin for i := 1 to n do begin if simp[i, j] < simp[l[j], j] then l[j] := i; if simp[i, j] > simp[h[j], j] then h[j] := i end end end; procedure sum_of_residuals (var x: vector); var i: index; begin x[n] := 0.0; for i := 1 to np do x[n] := x[n] + sqr(f(x, data[i]) - data[i, 2]) end; procedure Initialize; var i, j: index; firstx, firsty, lastx, lasty, xmean, ymean, slope, yintercept: extended; begin firstx := data[1, 1]; firsty := data[1, 2]; lastx := data[np, 1]; lasty := data[np, 2]; xmean := (firstx + lastx) / 2.0; ymean := (firsty + lasty) / 2.0; if (lastx - firstx) <> 0.0 then slope := (lasty - firsty) / (lastx - firstx) else slope := 1.0; yintercept := firsty - slope * firstx; case info^.fit of StraightLine: begin simp[1, 1] := yintercept; simp[1, 2] := slope; end; Poly2: begin simp[1, 1] := yintercept; simp[1, 2] := slope; simp[1, 3] := 0.0; end; Poly3: begin simp[1, 1] := yintercept; simp[1, 2] := slope; simp[1, 3] := 0.0; simp[1, 4] := 0.0; end; Poly4: begin simp[1, 1] := yintercept; simp[1, 2] := slope; simp[1, 3] := 0.0; simp[1, 4] := 0.0; simp[1, 5] := 0.0; end; ExpoFit: begin simp[1, 1] := 0.1; simp[1, 2] := 0.01; end; PowerFit: begin simp[1, 1] := 0.0; simp[1, 2] := 1.0; end; LogFit: begin simp[1, 1] := 0.5; simp[1, 2] := 0.05; end; RodbardFit: begin simp[1, 1] := firsty; simp[1, 2] := 1.0; simp[1, 3] := xmean; simp[1, 4] := lasty; end; end; maxiter := 100 * m * m; n := m + 1; for i := 1 to m do begin step[i] := simp[1, i] / 2.0; if step[i] = 0.0 then step[i] := 0.01; end; for i := 1 to n do maxerr[i] := MaxError; sum_of_residuals(simp[1]); for i := 1 to m do begin p[i] := step[i] * (sqrt(n) + m - 1) / (m * root2); q[i] := step[i] * (sqrt(n) - 1) / (m * root2) end; for i := 2 to n do begin for j := 1 to m do simp[i, j] := simp[i - 1, j] + q[j]; simp[i, i - 1] := simp[i, i - 1] + p[i - 1]; sum_of_residuals(simp[i]) end; for i := 1 to n do begin l[i] := 1; h[i] := 1 end; order; maxx := 255; end; procedure new_vertex; var i: index; begin for i := 1 to n do simp[h[n], i] := next[i] end; procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector); var i, j: integer; d: datarow; begin np := nStandards; m := nCoefficients; for i := 1 to np do begin data[i, 1] := xdata[i]; data[i, 2] := ydata[i]; end; Initialize; niter := 0; repeat done := true; niter := succ(niter); for i := 1 to n do center[i] := 0.0; for i := 1 to n do if i <> h[n] then for j := 1 to m do center[j] := center[j] + simp[i, j]; for i := 1 to n do begin center[i] := center[i] / m; next[i] := (1.0 + alpha) * center[i] - alpha * simp[h[n], i] end; sum_of_residuals(next); if next[n] <= simp[l[n], n] then begin new_vertex; for i := 1 to m do next[i] := gamma * simp[h[n], i] + (1.0 - gamma) * center[i]; sum_of_residuals(next); if next[n] <= simp[l[n], n] then new_vertex end else begin if next[n] <= simp[h[n], n] then new_vertex else begin for i := 1 to m do next[i] := beta * simp[h[n], i] + (1.0 - beta) * center[i]; sum_of_residuals(next); if (next[n] <= simp[h[n], n]) then new_vertex else begin for i := 1 to n do begin for j := 1 to m do simp[i, j] := (simp[i, j] + simp[l[n], j]) * beta; sum_of_residuals(simp[i]) end end end end; order; for j := 1 to n do begin if (simp[h[j], j] - simp[l[j], j]) = 0 then error[j] := 0 else if simp[h[j], j] <> 0 then error[j] := (simp[h[j], j] - simp[l[j], j]) / simp[h[j], j] else error[j] := (simp[h[j], j] - simp[l[j], j]) / simp[l[j], j]; if done then if abs(error[j]) > maxerr[j] then done := false end; until (done or (niter = maxiter)); ShowMessage(concat('interations=', long2str(niter), crStr, 'max interations=', long2str(maxiter))); for i := 1 to n do begin smean[i] := 0; for j := 1 to n do smean[i] := smean[i] + simp[j, i]; smean[i] := smean[i] / n; end; for i := 1 to m do Coefficients[i] := smean[i]; for i := 1 to nstandards do begin d[1] := xdata[i]; Residuals[i] := ydata[i] - f(smean, d); end; end; end.