
unit HookeJeevesMath;
{$R-}
interface

uses uMatToolsDynSingle33;

type tfunc=function(var x:vector):real;
// function to be minimized with respect to the parameters in x; define elsewhere

procedure HookeJeeves(var v1,delta:vector; function1:tfunc; var ff1,eps:real);
{ define function "function1(v1)" elsewhere; input parameters v1, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output ff1 and v1;
delta continuously reduced.
}

procedure HookeJeeves1(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta continuously reduced.
}

procedure NelderMead1(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta untouched.
}

procedure NelderMead3(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta untouched.
} /// !!! 3 calls to NelderMead !!!

implementation





procedure Explor1(nk:integer; var ff:real; var v,delta:vector; function1:tfunc);
// ff:function value at v preevaluated
// delta>0 is step size for each parameter in v
// on output, maybe new ff and v if minimum found
// input nk:=v.n;
var k:integer; f1:real;
begin
  for k:=1 to nk do begin
    v.v[k]:=v.v[k]+delta.v[k];
    f1:=function1(v);
    if f1<ff then
      ff:=f1
    else begin
      v.v[k]:=v.v[k]-2*delta.v[k];
      f1:=function1(v);
      if f1<ff then
        ff:=f1
      else
        v.v[k]:=v.v[k]+delta.v[k];
    end;
  end;
end;




procedure HookeJeeves(var v1,delta:vector; function1:tfunc; var ff1,eps:real);
{ define function "function1(v1)" elsewhere; input parameters v1, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output ff1 and v1;
delta continuously reduced.
}
label 1,2,3;
var k,nk:integer; ff,theta,mx,r:real; v:vector;
begin
  nk:=v1.n;
  create1(nk,v);

  ff1:=function1(v1);

  1:;
  ff:=ff1;
  equate1(v,v1);
  Explor1(nk,ff,v,delta,function1);

  if ff<ff1 then begin
    2:;
    for k:=1 to nk do begin
      theta:=v1.v[k];
      v1.v[k]:=v.v[k];
      v.v[k]:=2*v.v[k]-theta
    end;
    ff1:=ff;
    ff:=function1(v);
    if ff<ff1 then
      goto 2
    else
      goto 1;
  end;

  3:;
  mx:=-1e10;
  for k:=1 to nk do begin
    r:=delta.v[k]/2;
    delta.v[k]:=r;
    if mx<r then
      mx:=r;
  end;

  if mx<eps then
    exit
  else
    goto 1;

  destroy1(v);
end;




procedure HookeJeeves1(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta continuously reduced.
}
var k,nk:integer; f1,f2,r,mx:real; x1,x2:vector;
begin
  nk:=x0.n;
  create1(nk,x1);
  create1(nk,x2);
  f0:=function1(x0);
  repeat
    equate1(x1,x0);
    f1:=f0;
    Explor1(nk,f1,x1,delta,function1);
    if f1<f0 then begin
      for k:=1 to nk do
        x2.v[k]:=2*x1.v[k]-x0.v[k];
      f2:=function1(x2);
      if f2<f1 then begin
        f0:=f2;
        equate1(x0,x2)
      end else begin
        f0:=f1;
        equate1(x0,x1);
      end;
      mx:=1e20;
    end else begin
      mx:=-1e10;
      for k:=1 to nk do begin
        r:=delta.v[k]/2;
        delta.v[k]:=r;
        if mx<r then
          mx:=r;
      end;
    end;
  until mx<=eps;
  destroy1(x1);
  destroy1(x2);
end;




PROCEDURE NelderMead0(VAR p: matrix; // (ndim+1)*ndim; ndim+1 initial points in ndim space
                      VAR y: vector; // (ndim+1) function values at the initial points
                      ndim: integer;
                      ftol: real; // small value for covvergence
                      VAR nfunc: integer; // number of function evals
                      func:tfunc
                      );
LABEL 99;
CONST
   alpha = 1.0;
   beta = 0.5;
   gamma = 2.0;
VAR
   mpts,j,inhi,ilo,ihi,i: integer;
   ytry,ysave,sum,rtol: real;
   psum,ptry: vector;

FUNCTION NMtry(VAR p: matrix;
                VAR y,sum,ptry: vector;
                ndim,ihi: integer;
                VAR nfunc: integer;
                fac: real;
                func:tfunc
                ): real;
VAR
   j: integer;
   fac1,fac2,ytry: real;
BEGIN
   fac1 := (1.0-fac)/ndim;
   fac2 := fac1-fac;
   FOR j := 1 TO ndim DO
      ptry.v[j] := sum.v[j]*fac1-p.m[ihi,j]*fac2;
   ytry := func(ptry);
   nfunc := nfunc+1;
   IF ytry < y.v[ihi] THEN BEGIN
      y.v[ihi] := ytry;
      FOR j := 1 TO ndim DO BEGIN
         sum.v[j] := sum.v[j]+ptry.v[j]-p.m[ihi,j];
         p.m[ihi,j] := ptry.v[j]
      END
   END;
   NMtry := ytry;
END;

BEGIN
  create1(ndim,psum);
  create1(ndim,ptry);
   mpts := ndim+1;
   nfunc := 0;
   FOR j := 1 TO ndim DO BEGIN
      sum := 0.0;
      FOR i := 1 TO mpts DO
         sum := sum+p.m[i,j];
      psum.v[j] := sum
   END;
   WHILE true DO BEGIN
      ilo := 1;
      IF y.v[1] > y.v[kk2] THEN BEGIN
         ihi := 1;
         inhi := 2
      END
      ELSE BEGIN
         ihi := 2;
         inhi := 1
      END;
      FOR i := 1 TO mpts DO BEGIN
         IF y.v[i] < y.v[ilo] THEN ilo := i;
         IF y.v[i] > y.v[ihi] THEN BEGIN
            inhi := ihi;
            ihi := i
         END
         ELSE IF y.v[i] > y.v[inhi] THEN
            IF i <> ihi THEN inhi := i
      END;
      rtol := 2.0*abs(y.v[ihi]-y.v[ilo])/(abs(y.v[ihi])+abs(y.v[ilo]));
      IF rtol < ftol THEN GOTO 99;
      ytry := NMtry(p,y,psum,ptry,ndim,ihi,nfunc,-alpha,func);
      IF ytry <= y.v[ilo] THEN
         /// !!! never used says compiler: ytry := NMtry(p,y,psum,ptry,ndim,ihi,nfunc,gamma,func)
      ELSE IF ytry >= y.v[inhi] THEN BEGIN
         ysave := y.v[ihi];
         ytry := NMtry(p,y,psum,ptry,ndim,ihi,nfunc,beta,func);
         IF ytry >= ysave THEN BEGIN
            FOR i := 1 TO mpts DO
               IF i <> ilo THEN BEGIN
                  FOR j := 1 TO ndim DO BEGIN
                     psum.v[j] := 0.5*(p.m[i,j]+p.m[ilo,j]);
                     p.m[i,j] := psum.v[j]
                  END;
                  y.v[i] := func(psum)
               END;
            nfunc := nfunc+ndim;
            FOR j := 1 TO ndim DO BEGIN
               sum := 0.0;
               FOR i := 1 TO mpts DO
                  sum := sum+p.m[i,j];
               psum.v[j] := sum
            END
         END
      END
   END;
99:
   destroy1(psum);
   destroy1(ptry);
END;




procedure NelderMead1(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta untouched.
}
var p:matrix; i,j,ndim,nfunc:integer; y,xx:vector;
begin
  ndim:=x0.n;
  create1(ndim+1,ndim,p);
  create1(ndim+1,y);
  create1(ndim,xx);

  for j:=1 to ndim+1 do
    for i:=1 to ndim do
      p.m[j,i]:=x0.v[i];
  for j:=2 to ndim+1 do
    if random<0.5 then
      p.m[j,j-1]:=x0.v[j-1]+delta.v[j-1]
    else
      p.m[j,j-1]:=x0.v[j-1]-delta.v[j-1];

  for j:=1 to ndim+1 do begin
    for i:=1 to ndim do
      xx.v[i]:=p.m[j,i];
    y.v[j]:=function1(xx);
  end;

  NelderMead0(p,y,ndim,eps,nfunc,function1);
  for i:=1 to ndim do
    x0.v[i]:=p.m[1,i];
  f0:=y.v[1];

  destroy1(p);
  destroy1(y);
  destroy1(xx);
end;




procedure NelderMead3(var x0,delta:vector; function1:tfunc; var f0,eps:real);
{ define function "function1(x0)" elsewhere; input parameters x0, deltas>0 for
exploratory moves, eps which is compared to the largest delta; output f0 and x0;
delta untouched.
} /// !!! 3 calls to NelderMead !!!
begin
  NelderMead1(x0,delta,function1,f0,eps);
  NelderMead1(x0,delta,function1,f0,eps);
  NelderMead1(x0,delta,function1,f0,eps);
end;




procedure ScanForMin1(var x0:vector; function1:tfunc; var f0:real);
var p,i:integer; v:vector;
begin
  p:=x0.n;
  create1(p,v);

  for i:=-3 to 2 do

  destroy1(v);
end;




end.

