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.