[Back to MATH SWAG index] [Back to Main SWAG index] [Original]
{
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[ILo]) then
ILo:=I;
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])))<FTol) then
Exit;
{If we are not allowed to do any more work, return although
unsatisfactory}
if (Iter=MaxIterations) then
Exit;
{Do another iteration}
Inc(Iter);
{Compute the centroid of the face that leaves out the
one worst point}
for J:=1 to NDim do begin
STemp:=0;
for I:=1 to MPts do
if (I<>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[ILo]) then begin
{replace the highest point with PRR}
P[IHi]:=PRR;
Y[IHi]:=YPrr
end
else begin
{replace the highest point with PR}
P[IHI]:=PR;
Y[IHi]:=YPr
end
end {YPr<YLo}
else if (YPr>=Y[INHi]) then begin
{The new point is worse than the second highest, but it might
still be better than the highest}
if (YPr<Y[IHi]) then begin
P[IHi]:=PR;
Y[IHi]:=YPr
end;
{Whether or not the new point replaced the highest, see
whether a point of the interior, interpolating between
the (possibly new) highest point and the old centroid,
is better than the highest point, and if so, replace the
highest point; if not, contract around the best point so
far.}
for J:=1 to NDim do
PRR[J]:=Beta*P[IHi,J]+(1-Beta)*PBar[J];
YPrr:=Func(PRR,'Interp. reflection');
if (YPrr<Y[IHi]) then begin {replace}
P[IHi]:=PRR;
Y[IHi]:=YPrr
end
else {contract}
for I:=1 to MPts do
if (I<>ILo) 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.
[Back to MATH SWAG index] [Back to Main SWAG index] [Original]