```Contributor: HOWARD KAPLAN

{
In response to a request from rfv@maties.sun.ac.za in sunny South
Africa, here is some code from snowy Toronto which implements the
simplex method:

{SIMPLEX.PAS -- Simplex routine adapted from FORTRAN version in
Numerical Recipes, by Press, Flannery, Teukolsy, and
Howard L. Kaplan, Addiction Research Foundation,
Toronto, Ontario in 1993.  If you have the right to
use the FORTRAN original, then Howard Kaplan grants
you the right to use this adaptation.}

{\$N+}
{if you want a standalone test, then do not \$define UnitMode}
{\$define UnitMode}

{\$ifdef UnitMode}
unit    SIMPLEX;

interface
{\$endif}

const   MaxSimplexDimensions=10;
MaxIterations=200;

type    ParameterVector=array[1..MaxSimplexDimensions] of single;
ParameterArray=array[1..MaxSimplexDimensions+1] of
ParameterVector;
EvaluationArray=array[1..MaxSimplexDimensions+1] of single;
String20=string[20];
Evaluator=function(R: ParameterVector;
const Why: OpenString): single;

procedure Amoeba(var P: ParameterArray; {starting vertices}
NDim: integer; {N of elements of ParameterVector used}
FTol: single;  {fractional convergence tolerance}
Func: Evaluator; {fit evaluation function}
var Iter: integer {number of iterations used});
{\$ifdef UnitMode}
implementation
{\$else}
forward;
{\$endif}

procedure Amoeba(var P: ParameterArray; {starting vertices}
NDim: integer; {N of elements of ParameterVector used}
FTol: single;  {fractional convergence tolerance}
Func: Evaluator; {fit evaluation function}
var Iter: integer {number of iterations used});
const Alpha=1.0;
Beta=0.5;
Gamma=2.0;
var   Y:     EvaluationArray;
PBar:  ParameterVector; {centroid without worst point}
PR,
PRR:   ParameterVector; {test points to be evaluated}
MPts,              {number of points in the simplex}
I,J,
IHi,               {index of highest (worst) evaluation}
INHi,              {index of next-highest evaluation}
ILo:   integer;    {index of lowest (best) evaluation}
YPr,               {function evaluated at PR}
YPrr,              {function evaluated at PRR}
STemp: single;
Why:   String[20];
{sub}function NString(N: integer): string20;
var   S: String[20];
begin
Str(N,S);
while (Length(S)<4) do
S:=' '+S;
NString:=S
end;
begin
MPts:=NDim+1; {number of points in the simplex}
Iter:=0;
FillChar(Y,SizeOf(Y),0); {facilitate debugging}
for J:=1 to MPts do
Y[J]:=Func(P[J],'Initial row '+NString(J)); {initial simplex}
repeat
{Find the worst, next worst, and best vertices so far}
if (Y[1]>Y[2]) then begin
IHi:=1;
INHi:=2
end
else begin
IHi:=2;
INHi:=1
end;
ILo:=1;
for I:=1 to MPts do begin
if (Y[I]Y[IHi]) then begin
INHi:=IHi;
IHi:=I
end
else if (Y[I]>Y[INHi]) then
if (I<>IHi) then
INHi:=I
end;
{If the worst is 0 or close, return}
if (Y[IHi]<=FTol) then
Exit;
{Compute the fractional range from worst to best, and return if
satisfactory}
if ((2*Abs(Y[IHi]-Y[ILo])/(Abs(Y[IHi])+Abs(Y[ILo])))IHi) then
STemp:=STemp+P[I,J];
PBar[J]:=STemp/NDim
end;
{Reflect the simplex from the worst point, and evaluate}
for J:=1 to NDim do
PR[J]:=(1+Alpha)*PBar[J]-Alpha*P[IHi,J];
YPr:=Func(PR,'Reflect worst');
if (YPr<=Y[ILo]) then begin
{This is better than the best so far, so try an additional
extrapolation factor of Gamma}
for J:=1 to NDim do
PRR[J]:=Gamma*Pr[J]+(1-Gamma)*PBar[J];
YPrr:=Func(PRR,'Extend reflection');
if (YPrr=Y[INHi]) then begin
{The new point is worse than the second highest, but it might
still be better than the highest}
if (YPrILo) then begin
for J:=1 to NDim do
P[I,J]:=(P[I,J]+P[ILo,J])/2;
Y[I]:=Func(P[I],'Contract')
end
end {PRR was worse than the second-highest}
else begin
{PR was better than the second-highest, so we replace the old
highest point}
P[IHi]:=PR;
Y[IHi]:=YPr
end
until (false)
end;

{\$ifndef UnitMode}

function  TestEvaluation(R: ParameterVector;
const S: OpenString): single;
{Try to fit Y=A*Sqr(X)+B*X+C, where A=20, B=3, C=16, for -30<=X<=30}
var   N:     integer;
Sigma: single;
begin
Write('Evaluating A=',R[1]:6:3,', B=',R[2]:6:3,', C=',R[3]:6:3);
Sigma:=0;
for N:=-30 to 30 do
Sigma:=Sigma+Sqr((R[1]*N*N+R[2]*N+R[3])-(20*N*N+3*N+16));
WriteLn(' : ',Sigma:10:2);
TestEvaluation:=Sigma
end;

procedure TestAmoeba;
const Dimensions=3;
var   Simplex: ParameterArray;
Iter:    integer;
begin
WriteLn;
Simplex[1,1]:=1;
Simplex[1,2]:=1;
Simplex[1,3]:=1;
Simplex[2,1]:=-1;
Simplex[2,2]:=-1;
Simplex[2,3]:=-1;
Simplex[3,1]:=1;
Simplex[3,2]:=3;
Simplex[3,3]:=5;
Simplex[4,1]:=6;
Simplex[4,2]:=4;
Simplex[4,3]:=7;
Amoeba(Simplex,3,0.00001,TestEvaluation,Iter);
WriteLn('Converged to A=',Simplex[1,1]:6:3,', B=',Simplex[1,2]:6:3,
', C=',Simplex[1,3]:6:3,' after ',Iter,' iterations')
end;
{\$endif}

begin
{\$ifndef UnitMode}
TestAmoeba
{\$endif}
end.
```