{ **********************************************************************
  *                            Unit OPTIM.PAS                          *
  *                             Version 2.4d                           *
  *                     (c) J. Debord, February 2003                   *
  **********************************************************************
  This unit implements the following methods for function minimization:

    * Golden search for a function of one variable
    * Simplex, Marquardt, BFGS for a function of several variables
  **********************************************************************
  References:
  1) 'Numerical Recipes' by Press et al.
  2) D. W. MARQUARDT, J. Soc. Indust. Appl. Math., 1963, 11, 431-441
  3) J. A. NELDER & R. MEAD, Comput. J., 1964, 7, 308-313
  4) R. O'NEILL, Appl. Statist., 1971, 20, 338-345
  ********************************************************************** }

unit optim;

interface

uses
  fmath, matrices;

{ **********************************************************************
  Error codes
  ********************************************************************** }

const
  OPT_OK         =   0;  { No error }
  OPT_SING       = - 1;  { Singular hessian matrix }
  OPT_BIG_LAMBDA = - 2;  { Too high Marquardt's parameter }
  OPT_NON_CONV   = - 3;  { Non-convergence }

{ **********************************************************************
  Functional types
  ********************************************************************** }

type
  { Procedure to compute gradient vector }
  TGradient = procedure(Func           : TFuncNVar;
                        X              : TVector;
                        Lbound, Ubound : Integer;
                        G              : TVector);

  { Procedure to compute gradient vector and hessian matrix }
  THessGrad = procedure(Func           : TFuncNVar;
                        X              : TVector;
                        Lbound, Ubound : Integer;
                        G              : TVector;
                        H              : TMatrix);

{ **********************************************************************
  Log file
  ********************************************************************** }

const
  WriteLogFile : Boolean = False;        { Write iteration info to log file }
  LogFileName  : String  = 'optim.log';  { Name of log file }

{ **********************************************************************
  Minimization routines
  ********************************************************************** }

function GoldSearch(Func           : TFunc;
                    A, B           : Float;
                    MaxIter        : Integer;
                    Tol            : Float;
                    var Xmin, Ymin : Float) : Integer;
{ ----------------------------------------------------------------------
  Performs a golden search for the minimum of function Func
  ----------------------------------------------------------------------
  Input parameters  : Func    = objective function
                      A, B    = two points near the minimum
                      MaxIter = maximum number of iterations
                      Tol     = required precision (should not be less than
                                the square root of the machine precision)
  ----------------------------------------------------------------------
  Output parameters : Xmin, Ymin = coordinates of minimum
  ----------------------------------------------------------------------
  Possible results  : OPT_OK
                      OPT_NON_CONV
  ---------------------------------------------------------------------- }

function LinMin(Func           : TFuncNVar;
                X, DeltaX      : TVector;
                Lbound, Ubound : Integer;
                MaxIter        : Integer;
                Tol            : Float;
                var F_min      : Float) : Integer;

{ ----------------------------------------------------------------------
  Minimizes function Func from point X in the direction specified by
  DeltaX
  ----------------------------------------------------------------------
  Input parameters  : Func    = objective function
                      X       = initial minimum coordinates
                      DeltaX  = direction in which minimum is searched
                      Lbound,
                      Ubound  = indices of first and last variables
                      MaxIter = maximum number of iterations
                      Tol     = required precision
  ----------------------------------------------------------------------
  Output parameters : X     = refined minimum coordinates
                      F_min = function value at minimum
  ----------------------------------------------------------------------
  Possible results  : OPT_OK
                      OPT_NON_CONV
  ---------------------------------------------------------------------- }

function Simplex(Func           : TFuncNVar;
                 X              : TVector;
                 Lbound, Ubound : Integer;
                 MaxIter        : Integer;
                 Tol            : Float;
                 var F_min      : Float) : Integer;
{ ----------------------------------------------------------------------
  Minimization of a function of several variables by the simplex method
  of Nelder and Mead
  ----------------------------------------------------------------------
  Input parameters  : Func    = objective function
                      X       = initial minimum coordinates
                      Lbound,
                      Ubound  = indices of first and last variables
                      MaxIter = maximum number of iterations
                      Tol     = required precision
  ----------------------------------------------------------------------
  Output parameters : X     = refined minimum coordinates
                      F_min = function value at minimum
  ----------------------------------------------------------------------
  Possible results : OPT_OK
                     OPT_NON_CONV
  ---------------------------------------------------------------------- }

procedure NumGradient(Func           : TFuncNVar;
                      X              : TVector;
                      Lbound, Ubound : Integer;
                      G              : TVector);
{ ----------------------------------------------------------------------
  Computes the gradient vector of a function of several variables by
  numerical differentiation
  ----------------------------------------------------------------------
  Input parameters  : Func    = function of several variables
                      X       = vector of variables
                      Lbound,
                      Ubound  = indices of first and last variables
  ----------------------------------------------------------------------
  Output parameter  : G       = gradient vector
  ---------------------------------------------------------------------- }

procedure NumHessGrad(Func           : TFuncNVar;
                      X              : TVector;
                      Lbound, Ubound : Integer;
                      G              : TVector;
                      H              : TMatrix);
{ ----------------------------------------------------------------------
  Computes gradient vector & hessian matrix by numerical differentiation
  ----------------------------------------------------------------------
  Input parameters  : as in NumGradient
  ----------------------------------------------------------------------
  Output parameters : G = gradient vector
                      H = hessian matrix
  ---------------------------------------------------------------------- }

function Marquardt(Func           : TFuncNVar;
                   HessGrad       : THessGrad;
                   X              : TVector;
                   Lbound, Ubound : Integer;
                   MaxIter        : Integer;
                   Tol            : Float;
                   var F_min      : Float;
                   H_inv          : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Minimization of a function of several variables by Marquardt's method
  ----------------------------------------------------------------------
  Input parameters  : Func     = objective function
                      HessGrad = procedure to compute gradient & hessian
                      X        = initial minimum coordinates
                      Lbound,
                      Ubound   = indices of first and last variables
                      MaxIter  = maximum number of iterations
                      Tol      = required precision
  ----------------------------------------------------------------------
  Output parameters : X     = refined minimum coordinates
                      F_min = function value at minimum
                      H_inv = inverse hessian matrix
  ----------------------------------------------------------------------
  Possible results  : OPT_OK
                      OPT_SING
                      OPT_BIG_LAMBDA
                      OPT_NON_CONV
  ---------------------------------------------------------------------- }

function BFGS(Func           : TFuncNVar;
              Gradient       : TGradient;
              X              : TVector;
              Lbound, Ubound : Integer;
              MaxIter        : Integer;
              Tol            : Float;
              var F_min      : Float;
              H_inv          : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Minimization of a function of several variables by the
  Broyden-Fletcher-Goldfarb-Shanno method
  ----------------------------------------------------------------------
  Parameters : Gradient = procedure to compute gradient vector
               Other parameters as in Marquardt
  ----------------------------------------------------------------------
  Possible results : OPT_OK
                     OPT_NON_CONV
  ---------------------------------------------------------------------- }

implementation

var
  Eps              : Float;      { Fractional increment for numer. derivation }
  X1               : TVector;    { Initial point for line minimization }
  DeltaX1          : TVector;    { Direction for line minimization }
  Xt               : TVector;    { Minimum found by line minimization }
  Lbound1, Ubound1 : Integer;    { Bounds of X1 and DeltaX1 }
  LinObjFunc       : TFuncNVar;  { Objective function for line minimization }
  LogFile          : Text;       { Stores the result of each minimization step }


  procedure MinBrack(Func : TFunc; var A, B, C, Fa, Fb, Fc : Float);
{ ----------------------------------------------------------------------
  Given two points (A, B) this procedure finds a triplet (A, B, C)
  such that:

  1) A < B < C
  2) A, B, C are within the golden ratio
  3) Func(B) < Func(A) and Func(B) < Func(C).

  The corresponding function values are returned in Fa, Fb, Fc
  ---------------------------------------------------------------------- }

  begin
    if A > B then
      Swap(A, B);
    Fa := Func(A);
    Fb := Func(B);
    if Fb > Fa then
      begin
        Swap(A, B);
        Swap(Fa, Fb);
      end;
    C := B + GOLD * (B - A);
    Fc := Func(C);
    while Fc < Fb do
      begin
        A := B;
        B := C;
        Fa := Fb;
        Fb := Fc;
        C := B + GOLD * (B - A);
        Fc := Func(C);
      end;
    if A > C then
      begin
        Swap(A, C);
        Swap(Fa, Fc);
      end;
  end;

  function GoldSearch(Func           : TFunc;
                      A, B           : Float;
                      MaxIter        : Integer;
                      Tol            : Float;
                      var Xmin, Ymin : Float) : Integer;
  var
    C, Fa, Fb, Fc, F1, F2, MinTol, X0, X1, X2, X3 : Float;
    Iter                                          : Integer;
  begin
    MinTol := Sqrt(MACHEP);
    if Tol < MinTol then Tol := MinTol;
    MinBrack(Func, A, B, C, Fa, Fb, Fc);
    X0 := A;
    X3 := C;
    if (C - B) > (B - A) then
      begin
        X1 := B;
        X2 := B + CGOLD * (C - B);
        F1 := Fb;
        F2 := Func(X2);
      end
    else
      begin
        X1 := B - CGOLD * (B - A);
        X2 := B;
        F1 := Func(X1);
        F2 := Fb;
      end;
    Iter := 0;
    while (Iter <= MaxIter) and (Abs(X3 - X0) > Tol * (Abs(X1) + Abs(X2))) do
      if F2 < F1 then
        begin
          X0 := X1;
          X1 := X2;
          F1 := F2;
          X2 := X1 + CGOLD * (X3 - X1);
          F2 := Func(X2);
          Inc(Iter);
        end
      else
        begin
          X3 := X2;
          X2 := X1;
          F2 := F1;
          X1 := X2 - CGOLD * (X2 - X0);