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
                Vetterling, 1986 edition.  Adaptation made by
                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.