
unit uHookeJeeves;

{$MODE Delphi}

{$R-}
interface

uses uVecMat;

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

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

implementation




procedure Explor1(nk:integer; var ff:real; var v,delta:vvector; 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[k]:=v[k]+delta[k];
    f1:=function1(v);
    if f1<ff then
      ff:=f1
    else begin
      v[k]:=v[k]-2*delta[k];
      f1:=function1(v);
      if f1<ff then
        ff:=f1
      else
        v[k]:=v[k]+delta[k];
    end;
  end;
end;




var x1,x2:vvector;




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




end.

