unit LeastSquares;
{Public-domain code for solving least-squares problems.}

{ Compatible with IMAGE (an image-measurement program by Wayne Rasband at NIH)}
{ and a functional substitute for the TurboPascalª Numerical Methods Toolbox in}
{ that application.  Not all calling-options of the Numerical Methods package}
{ are implemented.}

{ Written September 1989 under MPW Pascal by}
{    Dr. John C. Russ}
{    North Carolina State University}
{    P.O. Box 7907}
{    Raleigh, NC 27695}

{ Source code for performing non-linear least squares analysis added January 1990}
{ to implement 4 parameter general curve fitting (Rodbard equation) }
{ Written for TML Pascal by}
{     Cary N. Mariash, M.D.}
{     Associate Professor of Medicine}
{     Box 91 Mayo Memorial Hospital}
{     University of Minnesota}
{     Minneapolis, MN 55455}
{     Compuserve: 71715,1331}
{     BitNet: mariasc@UMNHSNVE.BITNET}

interface


	const
		TNColumnSize = 30;
		TNRowSize = 30;
		maxmsize = 7;
		maxcsize = 30;

	type
		{ TNColumnVector = array[1..TNColumnSize] of EXTENDED;}
    {TNRowVector = array[1..TNRowSize] of EXTENDED; }
		FitType = (expo, fourier, log, poly, power, user, rodbard);
		vector = array[1..maxcsize] of extended;
		matrix = array[1..maxmsize, 1..maxmsize] of extended;
		ivector = array[1..maxmsize] of integer;
		matrixHndl = ^matrixPtr;
		matrixPtr = ^matrix;
		TNColumnVector = vector;
		TNRowVector = vector;
		TNmatrix = array[1..TNColumnSize] of TNRowVector;
		TNSquareMatrix = array[1..TNRowSize] of TNRowVector;
		TNString40 = string[40];

	var
		ochisq: extended;
		da, atry, beta, sig: vector;
		oneda: matrixHndl;
		ma, mfit: integer;
		alamda: extended;
		lista: ivector;

	procedure LeastSquares (nstandards: INTEGER; xdata, ydata: TNColumnVector; ncoefficients: INTEGER; var solution: TNRowVector; var yfit, residuals: TNColumnVector; var fitsd, variance: EXTENDED; var err: BYTE; typeoffit: FitType); {Poly,Expo,Power,Log}


implementation

	procedure gaussj (var a: matrix; n: integer; var b: matrix; m: integer);
	forward;

	procedure covsrt (var covar: matrix; ma: integer; lista: ivector; mfit: integer);
	forward;

	procedure ReportError (s: string);
	{ Implemented only when needed for debugging purposes }
	begin
	end;

	procedure fRodbard (x: extended; a: vector; var y: extended; var dyda: vector; na: integer);
	{ This is the 4 parameter general curve fit function as proposed by }
	{ Dr. David Rodbard at the NIH }
	{ The form of the equation is:  y = (D) + (A - D)/(1 + (x/C)^B) }
	{ where A = response at 0 dose, B = slope of curve (Hill coefficient) }
	{ C = ED50, and D = response at infinite dose }

	{ a is the vector of parameters, dyda is a vector of the 1st dervivates for each parameter }

		var
			ex: extended;
	begin
		if (x = 0.0) then
			ex := 0.0
		else
			ex := exp(ln(x / a[3]) * a[2]);  {ex:=(x/a[3] )**a[2]}

		y := a[1] - a[4];
		y := y / (1 + ex);
		y := y + a[4];

		dyda[1] := 1.0 / (1.0 + ex);
		if (x = 0.0) then
			dyda[2] := 0.0
		else begin
				dyda[2] := -(a[1] - a[4]) * ex * ln(x / a[3]);
				dyda[2] := dyda[2] / sqr(1.0 + ex);
			end;
		dyda[3] := a[2] * (a[1] - a[4]) * ex / a[3];
		dyda[3] := dyda[3] / sqr(1 + ex);
		dyda[4] := 1.0 - dyda[1];
	end;


	procedure mrqcof (x, y, sig: vector; ndata: integer; var a: vector; ma: integer; lista: ivector; mfit: integer; var alpha: matrix; var beta: vector; var chisq: extended; procedure funcs (x: extended; a: vector; var ymod: extended; var dyda: vector; ma: integer));

		var
			k, j, i: integer;
			ymod, wt, sig2i, dy: extended;
			dyda: vector;

	begin
		for j := 1 to mfit do begin
				for k := 1 to j do
					alpha[j][k] := 0.0;
				beta[j] := 0.0;
			end;
		chisq := 0.0;
		for i := 1 to ndata do begin
				funcs(x[i], a, ymod, dyda, ma);
				sig2i := 1.0 / (sig[i] * sig[i]);
				dy := y[i] - ymod;
				for j := 1 to mfit do begin
						wt := dyda[lista[j]] * sig2i;
						for k := 1 to j do
							alpha[j][k] := alpha[j][k] + wt * dyda[lista[k]];
						beta[j] := beta[j] + dy * wt;
					end;
				chisq := chisq + dy * dy * sig2i;
			end;
		for j := 2 to mfit do
			for k := 1 to j - 1 do
				alpha[k][j] := alpha[j][k];
	end;


{ Levenber-Marquardt method, attempting to reduce the value chisq of a fit between a set of}
{ points x[1..ndata],y[1..ndata] with individual standard deviations of sig[1..ndata],}
{ and a nonlinear function dependent on coefficients a[1..ma].  The array lista[1..ma]}
{ numbers the parameters a such that the first mfit elements correspond to values actually}
{ being adjusted; the remaining ma-mfit parameters are held fixed at their input value.}
{ The program returns current best-fit values for th ma fit parameters a, and chisq.}
{ The [1..mfit][1..mfit] elements of the arrays covar[1..ma][1..ma], alpha[1..ma][1..ma]}
{ are used as working space during most iterations.  Supply a routine funcs(x,a,yfit,dyda,ma)}
{ that evaluates the fitting function yfit, and its derivatives dyda[1..ma] with respect to}
{ the fitting parameters a at x.  On the first call provide an initial guess for the}
{ parameters a, and set alamda<0 for initialization (which then sets alamda:=0.001).}
{ If a step succeeds chisq becomes smaller and alamda decreases by a factor of 10.}
{ If a step fails alamda grows by a factor of 10.  You must call this routine repeatedly}
{ until convergence is achieved.  Then, make one final call with alamda = 0, so that}
{ covar returns the covariance matrix, and alpha the curvature matrix. }

	procedure mrqmin (var x: vector; var y: vector; var sig: vector; ndata: integer; var a: vector; ma: integer; var lista: ivector; mfit: integer; var covar: matrix; var alpha: matrix; var chisq: extended; procedure funcs (x: extended; a: vector; var ymod: extended; var dyda: vector; ma: integer); var alamda: extended);

		var
			k, kk, j, ihit: integer;

	begin
		if alamda < 0.0 then  { Initialize }
			begin
				oneda := matrixHndl(NewHandle(SizeOf(matrix)));
				kk := mfit + 1;
				for j := 1 to ma do begin
						ihit := 0;
						for k := 1 to mfit do begin
								if lista[k] = j then
									ihit := ihit + 1;
							end;
						if ihit = 0 then begin
								lista[kk] := j;
								kk := kk + 1;
							end
						else if ihit > 1 then
							ReportError('Bad LISTA permutation in MRQMIN-1');
					end;
				if (kk <> (ma + 1)) then
					ReportError('Bad LISTA permuation in MRQMIN-2');
				alamda := 0.001;
				mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, chisq, funcs);
				ochisq := chisq;
			end;

            { alter linearized fitting matrix, by augmenting diagonal elements. }
		for j := 1 to mfit do begin
				for k := 1 to mfit do
					covar[j][k] := alpha[j][k];
				covar[j][j] := alpha[j][j] * (1.0 + alamda);
				oneda^^[j][1] := beta[j];
			end;

            { matrix solution }
		Hlock(Handle(oneda));
		gaussj(covar, mfit, oneda^^, 1);
		for j := 1 to mfit do
			da[j] := oneda^^[j][1];
		HUnlock(Handle(oneda));

            { once converged evaluate covariance matrix with alamda = 0.0 }
		if (alamda = 0.0) then begin
				covsrt(covar, ma, lista, mfit);
				DisposHandle(Handle(oneda));
				exit(mrqmin);
			end;
		for j := 1 to ma do
			atry[j] := a[j];
		for j := 1 to mfit do
			atry[lista[j]] := a[lista[j]] + da[j];
		mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, chisq, funcs);
		if (chisq < ochisq) then begin
				alamda := alamda * 0.1;
				ochisq := chisq;
				for j := 1 to mfit do begin
						for k := 1 to mfit do
							alpha[j][k] := covar[j][k];
						beta[j] := da[j];
						a[lista[j]] := atry[lista[j]];
					end;
			end
		else begin
				alamda := alamda * 10.0;
				chisq := ochisq;
			end;
	end;

	procedure gaussj (var a: matrix; n: integer; var b: matrix; m: integer);

{ Linear equation solution by Gauss-Jordon elimination. a[1..n][1..n] is}
{  input matrix of n by n elements, b[1..n][1..m] is input matrix of size}
{  n by m conatining the m right-hand side vectors.  On ouput a is replaced by}
{  its matrix inverse, and b is replaced by the corresponding set of solution}
{  vectors. }

		var
			indxc, indxr, ipiv: ivector;
			i, icol, irow, j, k, l, ll: integer;
			big, dum, pivinv: extended;
			temp: extended;

	begin
		for j := 1 to n do
			ipiv[j] := 0;
		for i := 1 to n do begin
				big := 0.0;
				for j := 1 to n do begin
						if (ipiv[j] <> 1) then
							for k := 1 to n do begin
									if (ipiv[k] = 0) then
										if (abs(a[j][k]) >= big) then begin
												big := abs(a[j][k]);
												irow := j;
												icol := k;
											end;
								end
						else if (ipiv[k] > 1) then
							ReportError('GAUSSJ: Singular Matrix-1');
					end;
				ipiv[icol] := ipiv[icol] + 1;
				if (irow <> icol) then begin
						for l := 1 to n do begin
								temp := a[irow][l];
								a[irow][l] := a[icol][l];
								a[icol][l] := temp;
							end;
						for l := 1 to m do begin
								temp := b[irow][l];
								b[irow][l] := b[icol][l];
								b[icol][l] := temp;
							end;
					end;
				indxr[i] := irow;
				indxc[i] := icol;
				if (a[icol][icol] = 0.0) then
					ReportError('GAUSSJ: Singular Matrix-2');
				pivinv := 1.0 / a[icol][icol];
				a[icol][icol] := 1.0;
				for l := 1 to n do
					a[icol][l] := a[icol][l] * pivinv;
				for l := 1 to m do
					b[icol][l] := b[icol][l] * pivinv;
				for ll := 1 to n do begin
						if (ll <> icol) then begin
								dum := a[ll][icol];
								a[ll][icol] := 0.0;
								for l := 1 to n do
									a[ll][l] := a[ll][l] - (a[icol][l]) * dum;
								for l := 1 to m do
									b[ll][l] := b[ll][l] - (b[icol][l]) * dum;
							end;
					end;
			end;
		for l := n downto 1 do begin
				if (indxr[l] <> indxc[l]) then
					for k := 1 to n do begin
							temp := a[k][indxr[l]];
							a[k][indxr[l]] := a[k][indxc[l]];
							a[k][indxc[l]] := temp;
						end;
			end;
	end;

	procedure covsrt (var covar: matrix; ma: integer; lista: ivector; mfit: integer);

		var
			i, j: integer;
			swap: extended;

	begin
		for j := 1 to ma - 1 do begin
				for i := j + 1 to ma do
					covar[i][j] := 0.0;
			end;
		for i := 1 to mfit - 1 do begin
				for j := i + 1 to mfit do begin
						if lista[j] > lista[i] then
							covar[lista[j]][lista[i]] := covar[i][j]
						else
							covar[lista[i]][lista[j]] := covar[i][j];
					end;
			end;
		swap := covar[1][1];
		for j := 1 to ma do begin
				covar[1][j] := covar[j][j];
				covar[j][j] := 0.0;
			end;
		covar[lista[1]][lista[1]] := swap;
		for j := 2 to mfit do
			covar[lista[j]][lista[j]] := covar[1][j];
		for j := 2 to ma do begin
				for i := 1 to j - 1 do
					covar[i][j] := covar[j][i];
			end;
	end;


	procedure LeastSquares;

		const
			MaxN = 30;

		var
			i, j, k, L, nc, nx, ct: INTEGER;
			a: array[1..6, 1..7] of EXTENDED;
			s, t, valu: array[1..7] of EXTENDED;
			z, p, r, tval, temp: EXTENDED;
			yv, xv: array[1..MaxN] of EXTENDED;
			xa: array[1..5, 1..MaxN] of EXTENDED;
			iter: integer;
			test, oldchisq, oldlamda: extended;
			covar, alpha: matrix;

	begin
		{check limits}
		nc := ncoefficients - 1;

		if (nstandards < 2) then begin
				err := 1;
				EXIT(LeastSquares);
			end;
		if (nc < 1) or (nc > 5) then begin
				err := 2;
				EXIT(LeastSquares);
			end;
		if (nc >= nstandards) then begin
				err := 3;
				EXIT(LeastSquares);
			end;

		ct := nstandards;
		if (ct > MaxN) then
			ct := MaxN;

		{initialize arrays}
		for j := 1 to 6 do
			for k := 1 to 7 do
				a[j, k] := 0;

		for j := 1 to 7 do begin
				s[j] := 0;
				t[j] := 0;
				valu[j] := 0;
			end;

		{copy incoming data and prepare columns of values}
		for i := 1 to ct do begin
				xv[i] := xdata[i];
				yv[i] := ydata[i];
				solution[i] := 0;
				yfit[i] := 0;
				residuals[i] := 0;
			end;
		fitsd := 0;

		if typeoffit = Rodbard then begin
					{ for i:=1 to nstandards do sig[i]:= (abs(yData[i]))/2; }
					{ the line above weights error by the value }
					{ the line below makes the fit unweighted }
				for i := 1 to nstandards do
					sig[i] := 1.0;
				ma := ncoefficients;
				alamda := -1.0;
				mfit := 4;
				for i := 1 to mfit do
					lista[i] := i;
				Solution[1] := YData[1];
				Solution[2] := 1.0;
				Solution[3] := XData[nstandards div 2];
				Solution[4] := YData[nstandards];
				iter := 1;
				mrqmin(XData, YData, sig, nstandards, Solution, ma, lista, mfit, covar, alpha, Variance, fRodbard, alamda);
				repeat
					oldchisq := Variance;
					oldlamda := alamda;
					mrqmin(XData, YData, sig, nstandards, Solution, ma, lista, mfit, covar, alpha, Variance, fRodbard, alamda);
					iter := iter + 1;
					test := (oldchisq - Variance) / oldchisq;
				until ((iter = 40) or ((test > 0) and (test < 0.00001)));
				alamda := 0.0;
				mrqmin(XData, YData, sig, nstandards, Solution, ma, lista, mfit, covar, alpha, Variance, fRodbard, alamda);
				for i := 1 to nstandards do begin
						fRodbard(XData[i], Solution, test, sig, ncoefficients);
						Residuals[i] := YData[i] - test;
					end;
				err := 0;
				exit(LeastSquares);
			end;{ if fit Rodbard }


		{if type of fit is poly, create x^2, x^3, ...}
		if (typeoffit = poly) then
			if nc > 1 then
				for i := 1 to ct do begin
						tval := xv[i];
						temp := tval;
						for j := 2 to nc do begin
								tval := tval * temp;
								xa[j, i] := tval;
							end;
					end;

		{if type of fit is log, y=a ln(bx) then create ln(x)}
{		 - solve as y= a*ln(b) + a*ln(x)}
		if typeoffit = log then
			for i := 1 to ct do begin
					tval := xv[i];
					if tval > 0 then
						xv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{if type of fit is power, y=ax^b then create ln(y) and ln(x)}
{		  - solve as ln(y)=ln(a)+b*ln(x)}
		if typeoffit = power then
			for i := 1 to ct do begin
					tval := xv[i];
					if tval > 0 then
						xv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
					tval := yv[i];
					if tval > 0 then
						yv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{if type of fit is expo, y=ae^bx then create ln(y)}
{		 - solve as ln(y)=ln(a)+bx}
		if typeoffit = expo then
			for i := 1 to ct do begin
					tval := yv[i];
					if tval > 0 then
						yv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{Build the matrix}
		nx := ncoefficients;
		for i := 1 to ct do begin
				valu[1] := 1.0;
				valu[2] := xv[i];
				if (typeoffit = poly) and (nc > 1) then
					for j := 2 to nc do
						valu[j + 1] := xa[j, i];
							{subscript order?}

				valu[nx + 1] := yv[i];
				for k := 1 to nx do
					for L := 1 to nx + 1 do begin
							a[k, L] := a[k, L] + valu[k] * valu[L];
							s[k] := a[k, nx + 1];
						end;

				s[nx + 1] := s[nx + 1] + valu[nx + 1] * valu[nx + 1];
			end; {for i}

		{Solve the matrix by inversion}
		for i := 2 to nx do
			t[i] := a[1, i];
		for i := 1 to nx do begin
				j := i - 1;
				repeat
					j := j + 1;
				until (j > nx) or (a[j, i] <> 0);

				if j > nx then begin
						err := 4;
						EXIT(LeastSquares);
					end;

				for k := 1 to nx + 1 do begin
						p := a[i, k];
						a[i, k] := a[j, k];
						a[j, k] := p;
					end;

				z := 1 / a[i, i];
				for k := 1 to nx + 1 do
					a[i, k] := z * a[i, k];

				for j := 1 to nx do
					if j <> i then begin
							z := -a[j, i];
							for k := 1 to nx + 1 do
								a[j, k] := a[j, k] + z * a[i, k];
						end;
			end; {for i=1..nx}

		{copy back the solution vector}
		for j := 1 to nx do
			solution[j] := a[j, nx + 1];

		if typeoffit = log then begin
				if solution[2] <= 0 then begin
						err := 4;
						EXIT(LeastSquares);
					end
				else begin
						p := solution[2];
						solution[2] := EXP(solution[1] / p);
						solution[1] := p;
					end;
			end;
		{solved as y= a*ln(b)+a*ln(x)}

		if typeoffit = power then
			solution[1] := EXP(solution[1]);
		{solved as ln(y)=ln(a)+b*ln(x)}

		if typeoffit = expo then
			solution[1] := EXP(solution[1]);
		{solved as ln(y)=ln(a)+b*ln(x)}

		{Create the predicted and residual columns}
		for i := 1 to ct do begin
				valu[1] := xv[i];
				if nc > 1 then
					for j := 2 to nc do
						valu[j] := xa[j, i];

				valu[nx] := yv[i];
				p := solution[1];
				for j := 1 to nx - 1 do
					p := p + solution[j + 1] * valu[j];

				if (typeoffit = power) or (typeoffit = expo) then
					p := EXP(p); {convert from ln(y)}
				yfit[i] := p;
				residuals[i] := valu[nx] - p;
			end;
		p := 0;
		for i := 2 to nx do
			p := p + a[i, nx + 1] * (s[i] - t[i] * s[1] / ct);
		r := s[nx + 1] - s[1] * s[1] / ct;
		z := r - p;
		L := ct - nx;

		{Coeff. of determination (R-squared)= p/r}
    {Std error of estimate=Sqrt(Abs(Z/L))}
		variance := z / L;
		fitsd := SQRT(ABS(z / L));
		err := 0;

	end; {LeastSquares}

end. {unit leastsq}