TEXT   101

unit optim pas

Guest on 23rd April 2022 01:21:29 AM

  1. { **********************************************************************
  2.   *                            Unit OPTIM.PAS                          *
  3.   *                              Version 1.8                           *
  4.   **********************************************************************
  5.   This unit implements the following methods for function minimization:
  6.  
  7.     * Brent's method for a function of 1 variable
  8.     * Marquardt's, BFGS, Simplex and Simulated Annealing methods
  9.       for a function of several variables
  10.   **********************************************************************
  11.   References:
  12.   1) 'Numerical Recipes' by Press et al.
  13.   2) D. W. MARQUARDT, J. Soc. Indust. Appl. Math., 1963, 11, 431-441
  14.   3) J. A. NELDER & R. MEAD, Comput. J., 1964, 7, 308-313
  15.   4) R. O'NEILL, Appl. Statist., 1971, 20, 338-345
  16.   ********************************************************************** }
  17.  
  18. unit Optim;
  19.  
  20. interface
  21.  
  22. uses
  23.   FMath, Matrices, Random, Stat;
  24.  
  25. { **********************************************************************
  26.   Error codes
  27.   ********************************************************************** }
  28.  
  29. const
  30.   BIG_LAMBDA = - 2;  { Too high Marquardt's parameter }
  31.   NON_CONV   = - 3;  { Non-convergence }
  32.  
  33. { **********************************************************************
  34.   Functional types
  35.   ********************************************************************** }
  36.  
  37. type
  38.   { Function of one variable }
  39.   TFunc = function(X : Float) : Float;
  40.  
  41.   { Function of several variables }
  42.   TFuncNVar = function(X : PVector) : Float;
  43.  
  44.   { Procedure to compute gradient vector }
  45.   TGradient = procedure(Func           : TFuncNVar;
  46.                         X              : PVector;
  47.                         Lbound, Ubound : Integer;
  48.                         G              : PVector);
  49.  
  50.   { Procedure to compute gradient vector and hessian matrix }
  51.   THessGrad = procedure(Func           : TFuncNVar;
  52.                         X              : PVector;
  53.                         Lbound, Ubound : Integer;
  54.                         G              : PVector;
  55.                         H              : PMatrix);
  56.  
  57. { **********************************************************************
  58.   Log file
  59.   ********************************************************************** }
  60.  
  61. const
  62.   WriteLogFile : Boolean = False;        { Write iteration info to log file }
  63.   LogFileName  : String  = 'OPTIM.LOG';  { Name of log file }
  64.  
  65. { **********************************************************************
  66.   Parameters for simulated annealing
  67.   ********************************************************************** }
  68.  
  69. const
  70.   SA_Nt        : Integer = 50;   { Nb of iterations at each temperature }
  71.   SA_Cv        : Float   = 0.5;  { Coef. of variation for random numbers }
  72.   SA_MinSigma  : Float   = 1.0;  { Minimal S.D. for random numbers }
  73.   SA_SigmaFact : Float   = 0.1;  { Overall S.D. reduction factor }
  74.   SA_NCycles   : Integer = 1;    { Number of SA cycles }
  75.  
  76. { **********************************************************************
  77.   Minimization routines
  78.   ********************************************************************** }
  79.  
  80. function Brent(Func             : TFunc;
  81.                X1, X2           : Float;
  82.                MaxIter          : Integer;
  83.                Tol              : Float;
  84.                var X_min, F_min : Float) : Integer;
  85. { ----------------------------------------------------------------------
  86.   Minimization of a function of 1 variable by Brent's method
  87.   ----------------------------------------------------------------------
  88.   Input parameters  : Func    = function to be minimized
  89.                       X1, X2  = two values near the minimum
  90.                       MaxIter = maximum number of iterations
  91.                       Tol     = required precision
  92.   ----------------------------------------------------------------------
  93.   Output parameters : X_min = refined minimum coordinate
  94.                       F_min = function value at minimum
  95.   ----------------------------------------------------------------------
  96.   Possible results  : MAT_OK
  97.                       NON_CONV
  98.   ---------------------------------------------------------------------- }
  99.  
  100. procedure NumGradient(Func           : TFuncNVar;
  101.                       X              : PVector;
  102.                       Lbound, Ubound : Integer;
  103.                       G              : PVector);
  104. { ----------------------------------------------------------------------
  105.   Computes the gradient vector of a function of several variables by
  106.   numerical differentiation
  107.   ----------------------------------------------------------------------
  108.   Input parameters  : Func    = function of several variables
  109.                       X       = vector of variables
  110.                       Lbound,
  111.                       Ubound  = indices of first and last variables
  112.   ----------------------------------------------------------------------
  113.   Output parameter  : G       = gradient vector
  114.   ---------------------------------------------------------------------- }
  115.  
  116. procedure NumHessGrad(Func           : TFuncNVar;
  117.                       X              : PVector;
  118.                       Lbound, Ubound : Integer;
  119.                       G              : PVector;
  120.                       H              : PMatrix);
  121. { ----------------------------------------------------------------------
  122.   Computes gradient vector & hessian matrix by numerical differentiation
  123.   ----------------------------------------------------------------------
  124.   Input parameters  : as in NumGradient
  125.   ----------------------------------------------------------------------
  126.   Output parameters : G = gradient vector
  127.                       H = hessian matrix
  128.   ---------------------------------------------------------------------- }
  129.  
  130. function LinMin(Func           : TFuncNVar;
  131.                 X, DeltaX      : PVector;
  132.                 Lbound, Ubound : Integer;
  133.                 MaxIter        : Integer;
  134.                 Tol            : Float;
  135.                 var F_min      : Float) : Integer;
  136. { ----------------------------------------------------------------------
  137.   Minimizes function Func from point X in the direction specified by
  138.   DeltaX
  139.   ----------------------------------------------------------------------
  140.   Input parameters  : Func    = objective function
  141.                       X       = initial minimum coordinates
  142.                       DeltaX  = direction in which minimum is searched
  143.                       Lbound,
  144.                       Ubound  = indices of first and last variables
  145.                       MaxIter = maximum number of iterations
  146.                       Tol     = required precision
  147.   ----------------------------------------------------------------------
  148.   Output parameters : X     = refined minimum coordinates
  149.                       F_min = function value at minimum
  150.   ----------------------------------------------------------------------
  151.   Possible results  : MAT_OK
  152.                       NON_CONV
  153.   ---------------------------------------------------------------------- }
  154.  
  155. function Marquardt(Func           : TFuncNVar;
  156.                    HessGrad       : THessGrad;
  157.                    X              : PVector;
  158.                    Lbound, Ubound : Integer;
  159.                    MaxIter        : Integer;
  160.                    Tol            : Float;
  161.                    var F_min      : Float;
  162.                    H_inv          : PMatrix) : Integer;
  163. { ----------------------------------------------------------------------
  164.   Minimization of a function of several variables by Marquardt's method
  165.   ----------------------------------------------------------------------
  166.   Input parameters  : HessGrad = procedure to compute
  167.                                  gradient and hessian
  168.                       Other parameters as in LinMin
  169.   ----------------------------------------------------------------------
  170.   Output parameters : X     = refined minimum coordinates
  171.                       F_min = function value at minimum
  172.                       H_inv = inverse hessian matrix
  173.   ----------------------------------------------------------------------
  174.   Possible results  : MAT_OK
  175.                       MAT_SINGUL
  176.                       BIG_LAMBDA
  177.                       NON_CONV
  178.   ---------------------------------------------------------------------- }
  179.  
  180. function BFGS(Func           : TFuncNVar;
  181.               Gradient       : TGradient;
  182.               X              : PVector;
  183.               Lbound, Ubound : Integer;
  184.               MaxIter        : Integer;
  185.               Tol            : Float;
  186.               var F_min      : Float;
  187.               H_inv          : PMatrix) : Integer;
  188. { ----------------------------------------------------------------------
  189.   Minimization of a function of several variables by the
  190.   Broyden-Fletcher-Goldfarb-Shanno method
  191.   ----------------------------------------------------------------------
  192.   Parameters : Gradient = procedure to compute gradient vector
  193.                Other parameters as in Marquardt
  194.   ----------------------------------------------------------------------
  195.   Possible results : MAT_OK
  196.                      NON_CONV
  197.   ---------------------------------------------------------------------- }
  198.  
  199. function Simplex(Func           : TFuncNVar;
  200.                  X              : PVector;
  201.                  Lbound, Ubound : Integer;
  202.                  MaxIter        : Integer;
  203.                  Tol            : Float;
  204.                  var F_min      : Float) : Integer;
  205. { ----------------------------------------------------------------------
  206.   Minimization of a function of several variables by the simplex method
  207.   of Nelder and Mead
  208.   ----------------------------------------------------------------------
  209.   Parameters : as in Marquardt
  210.   ----------------------------------------------------------------------
  211.   Possible results : MAT_OK
  212.                      NON_CONV
  213.   ---------------------------------------------------------------------- }
  214.  
  215.   function SimAnn(Func           : TFuncNVar;
  216.                   X              : PVector;
  217.                   Lbound, Ubound : Integer;
  218.                   MaxIter        : Integer;
  219.                   T_Ratio        : Float;
  220.                   var F_min      : Float) : Integer;
  221. { ----------------------------------------------------------------------
  222.   Minimization of a function of several variables by simulated annealing
  223.   ----------------------------------------------------------------------
  224.   Input parameters : MaxIter = number of annealing steps
  225.                      T_Ratio = ratio of final to initial temperature
  226.                      Other parameters as in marquardt
  227.   ----------------------------------------------------------------------
  228.   Output parameter : F_min   = function value
  229.   ----------------------------------------------------------------------
  230.   The function always returns 0
  231.   ---------------------------------------------------------------------- }
  232.  
  233. implementation
  234.  
  235. var
  236.   Eps              : Float;      { Fractional increment for numer. derivation }
  237.   X1               : PVector;    { Initial point for line minimization }
  238.   DeltaX1          : PVector;    { Direction for line minimization }
  239.   Lbound1, Ubound1 : Integer;    { Bounds of X1 and DeltaX1 }
  240.   LinObjFunc       : TFuncNVar;  { Objective function for line minimization }
  241.   LogFile          : Text;       { Log file }
  242.  
  243.   procedure MinBrak(Func : TFunc; var A, B, C : Float);
  244. { --------------------------------------------------------------------
  245.   Procedure used by Brent's method to initially bracket the minimum.
  246.   Given two points A and B, this routine finds
  247.   a triplet of points (A < B < C) such that :
  248.   Func(B) < Func(A) and Func(B) < Func(C)
  249.   -------------------------------------------------------------------- }
  250.   label 1;
  251.   const
  252.     GLIMIT = 100.0;
  253.     TINY   = 1.0E-20;
  254.   var
  255.     Fa, Fb, Fc, Ulim, U, R, Q, D, Fu : Float;
  256.   begin
  257.     Fa := Func(A);
  258.     Fb := Func(B);
  259.     if Fb > Fa then
  260.       begin
  261.         FSwap(A, B);
  262.         FSwap(Fa, Fb);
  263.       end;
  264.     C := B + GOLD * (B - A);
  265.     Fc := Func(C);
  266. 1:  if Fb >= Fc then
  267.       begin
  268.         R := (B - A) * (Fb - Fc);
  269.         Q := (B - C) * (Fb - Fa);
  270.         D := Q - R;
  271.         if Abs(D) < TINY then D := Sgn(D) * TINY;
  272.         U := B - ((B - C) * Q - (B - A) * R) / (2.0 * D);
  273.         Ulim := B + GLIMIT * (C - B);
  274.         if (B - U) * (U - C) > 0.0 then
  275.           begin
  276.             Fu := Func(U);
  277.             if Fu < Fc then
  278.               begin
  279.                 A := B;
  280.                 Fa := Fb;
  281.                 B := U;
  282.                 Fb := Fu;
  283.                 goto 1
  284.               end
  285.             else if Fu > Fb then
  286.               begin
  287.                 C := U;
  288.                 Fc := Fu;
  289.                 goto 1
  290.               end;
  291.             U := C + GOLD * (C - B);
  292.             Fu := Func(U)
  293.           end
  294.         else if (C - U) * (U - Ulim) > 0.0 then
  295.           begin
  296.             Fu := Func(U);
  297.             if Fu < Fc then
  298.               begin
  299.                 B := C;
  300.                 C := U;
  301.                 U := C + GOLD * (C - B);
  302.                 Fb := Fc;
  303.                 Fc := Fu;
  304.                 Fu := Func(U)
  305.               end
  306.           end
  307.         else if (U - Ulim) * (Ulim - C) >= 0.0 then
  308.           begin
  309.             U := Ulim;
  310.             Fu := Func(U)
  311.           end
  312.         else
  313.           begin
  314.             U := C + GOLD * (C - B);
  315.             Fu := Func(U)
  316.           end;
  317.         A := B;
  318.         B := C;
  319.         C := U;
  320.         Fa := Fb;
  321.         Fb := Fc;
  322.         Fc := Fu;
  323.         goto 1
  324.       end
  325.   end;
  326.  
  327.   function Brent(Func             : TFunc;
  328.                  X1, X2           : Float;
  329.                  MaxIter          : Integer;
  330.                  Tol              : Float;
  331.                  var X_min, F_min : Float) : Integer;
  332.   label 1, 2, 3;
  333.   const
  334.     TINY = 1.0E-10;
  335.   var
  336.     A, B, D, E, Etemp, Fu, Fv, Fw, Fx       : Float;
  337.     P, Q, R, Tol1, Tol2, U, V, W, X, X3, Xm : Float;
  338.     Iter                                    : Integer;
  339.   begin
  340.     MinBrak(Func, X1, X2, X3);
  341.     if X1 < X3 then A := X1 else A := X3;
  342.     if X1 > X3 then B := X1 else B := X3;
  343.     V := X2;
  344.     W := V;
  345.     X := V;
  346.     D := 0.0;
  347.     E := 0.0;
  348.     Fx := Func(X);
  349.     Fv := Fx;
  350.     Fw := Fx;
  351.     for Iter := 1 to MaxIter do
  352.       begin
  353.         Xm := 0.5 * (A + B);
  354.         Tol1 := Tol * Abs(X) + TINY;
  355.         Tol2 := 2.0 * Tol1;
  356.         if Abs(X - Xm) <= (Tol2 - 0.5 * (B - A)) then goto 3;
  357.         if Abs(E) > Tol1 then
  358.           begin
  359.             R := (X - W) * (Fx - Fv);
  360.             Q := (X - V) * (Fx - Fw);
  361.             P := (X - V) * Q - (X - W) * R;
  362.             Q := 2.0 * (Q - R);
  363.             if Q > 0.0 then P := - P;
  364.             Q := Abs(Q);
  365.             Etemp := E;
  366.             E := D;
  367.             if (Abs(P) >= Abs(0.5 * Q * Etemp)) or (P <= Q * (A - X))
  368.             or (P >= Q * (B - X)) then goto 1;
  369.             D := P / Q;
  370.             U := X + D;
  371.             if ((U - A) < Tol2) or ((B - U) < Tol2) then
  372.               D := Sgn(Xm - X) * Abs(Tol1);
  373.             goto 2
  374.           end;
  375. 1:      if X >= Xm then E := A - X else E := B - X;
  376.         D := CGOLD * E;
  377. 2:      if Abs(D) >= Tol1 then U := X + D else U := X + Sgn(D) * Abs(Tol1);
  378.         Fu := Func(U);
  379.         if Fu <= Fx then
  380.           begin
  381.             if U >= X then A := X else B := X;
  382.             V := W;
  383.             Fv := Fw;
  384.             W := X;
  385.             Fw := Fx;
  386.             X := U;
  387.             Fx := Fu
  388.           end
  389.         else
  390.           begin
  391.             if U < X then A := U else B := U;
  392.             if (Fu <= Fw) or (W = X) then
  393.               begin
  394.                 V := W;
  395.                 Fv := Fw;
  396.                 W := U;
  397.                 Fw := Fu
  398.               end
  399.             else if (Fu <= Fv) or (V = X) or (V = 2) then
  400.               begin
  401.                 V := U;
  402.                 Fv := Fu
  403.               end
  404.           end
  405.       end;
  406.     Brent := NON_CONV;
  407. 3:  X_min := X;
  408.     F_min := Fx;
  409.     Brent := MAT_OK;
  410.   end;
  411.  
  412.   {$F+}
  413.   procedure NumGradient(Func           : TFuncNVar;
  414.                         X              : PVector;
  415.                         Lbound, Ubound : Integer;
  416.                         G              : PVector);
  417.   var
  418.     Temp, Delta, Fplus, Fminus : Float;
  419.     I                          : Integer;
  420.   begin
  421.     for I := Lbound to Ubound do
  422.       begin
  423.         Temp := X^[I];
  424.         if Temp <> 0.0 then Delta := Eps * Abs(Temp) else Delta := Eps;
  425.         X^[I] := Temp - Delta;
  426.         Fminus := Func(X);
  427.         X^[I] := Temp + Delta;
  428.         Fplus := Func(X);
  429.         G^[I] := (Fplus - Fminus) / (2.0 * Delta);
  430.         X^[I] := Temp;
  431.       end;
  432.   end;
  433.   {$F-}
  434.  
  435.   {$F+}
  436.   procedure NumHessGrad(Func           : TFuncNVar;
  437.                         X              : PVector;
  438.                         Lbound, Ubound : Integer;
  439.                         G              : PVector;
  440.                         H              : PMatrix);
  441.   var
  442.     Delta, Xminus, Xplus, Fminus, Fplus : PVector;
  443.     Temp1, Temp2, F, F2plus             : Float;
  444.     I, J                                : Integer;
  445.   begin
  446.     DimVector(Delta, Ubound);   { Increments   }
  447.     DimVector(Xminus, Ubound);  { X - Delta    }
  448.     DimVector(Xplus, Ubound);   { X + Delta    }
  449.     DimVector(Fminus, Ubound);  { F(X - Delta) }
  450.     DimVector(Fplus, Ubound);   { F(X + Delta) }
  451.  
  452.     F := Func(X);
  453.  
  454.     for I := Lbound to Ubound do
  455.       begin
  456.         if X^[I] <> 0.0 then
  457.           Delta^[I] := Eps * Abs(X^[I])
  458.         else
  459.           Delta^[I] := Eps;
  460.         Xplus^[I] := X^[I] + Delta^[I];
  461.         Xminus^[I] := X^[I] - Delta^[I];
  462.       end;
  463.  
  464.     for I := Lbound to Ubound do
  465.       begin
  466.         Temp1 := X^[I];
  467.         X^[I] := Xminus^[I];
  468.         Fminus^[I] := Func(X);
  469.         X^[I] := Xplus^[I];
  470.         Fplus^[I] := Func(X);
  471.         X^[I] := Temp1;
  472.       end;
  473.  
  474.     for I := Lbound to Ubound do
  475.       begin
  476.         G^[I] := (Fplus^[I] - Fminus^[I]) / (2.0 * Delta^[I]);
  477.         H^[I]^[I] := (Fplus^[I] + Fminus^[I] - 2.0 * F) / Sqr(Delta^[I]);
  478.       end;
  479.  
  480.     for I := Lbound to Pred(Ubound) do
  481.       begin
  482.         Temp1 := X^[I];
  483.         X^[I] := Xplus^[I];
  484.         for J := Succ(I) to Ubound do
  485.           begin
  486.             Temp2 := X^[J];
  487.             X^[J] := Xplus^[J];
  488.             F2plus := Func(X);
  489.             H^[I]^[J] := (F2plus - Fplus^[I] - Fplus^[J] + F) /
  490.               (Delta^[I] * Delta^[J]);
  491.             H^[J]^[I] := H^[I]^[J];
  492.             X^[J] := Temp2;
  493.           end;
  494.         X^[I] := Temp1;
  495.       end;
  496.  
  497.     DelVector(Delta, Ubound);
  498.     DelVector(Xminus, Ubound);
  499.     DelVector(Xplus, Ubound);
  500.     DelVector(Fminus, Ubound);
  501.     DelVector(Fplus, Ubound);
  502.   end;
  503.   {$F-}
  504.  
  505.   {$F+}
  506.   function F1dim(R : Float) : Float;
  507. { ----------------------------------------------------------------------
  508.   Function used by LinMin to find the minimum of the objective function
  509.   LinObjFunc in the direction specified by the global variables X1 and
  510.   DeltaX1. R is the step in this direction.
  511.   ---------------------------------------------------------------------- }
  512.   const
  513.     Xt : PVector = nil;
  514.   var
  515.     I : Integer;
  516.   begin
  517.     if Xt = nil then
  518.       DimVector(Xt, Ubound1);
  519.     for I := Lbound1 to Ubound1 do
  520.       Xt^[I] := X1^[I] + R * DeltaX1^[I];
  521.     F1dim := LinObjFunc(Xt);
  522.   end;
  523.   {$F-}
  524.  
  525.   function LinMin(Func           : TFuncNVar;
  526.                   X, DeltaX      : PVector;
  527.                   Lbound, Ubound : Integer;
  528.                   MaxIter        : Integer;
  529.                   Tol            : Float;
  530.                   var F_min      : Float) : Integer;
  531.   var
  532.     I, ErrCode : Integer;
  533.     R          : Float;
  534.   begin
  535.     { Redimension global vectors }
  536.     DelVector(X1, Ubound1);
  537.     DelVector(DeltaX1, Ubound1);
  538.     DimVector(X1, Ubound);
  539.     DimVector(DeltaX1, Ubound);
  540.  
  541.     Lbound1 := Lbound;
  542.     Ubound1 := Ubound;
  543.  
  544.     { Initialize global variables }
  545.     LinObjFunc := Func;
  546.     for I := Lbound to Ubound do
  547.       begin
  548.         X1^[I] := X^[I];
  549.         DeltaX1^[I] := DeltaX^[I]
  550.       end;
  551.  
  552.     { Call Brent's method }
  553.     {$IFDEF FPK}
  554.       ErrCode := Brent(@F1dim, 0.0, 1.0, MaxIter, Tol, R, F_min);
  555.     {$ELSE}
  556.       ErrCode := Brent(F1dim, 0.0, 1.0, MaxIter, Tol, R, F_min);
  557.     {$ENDIF}
  558.  
  559.     { Update variables }
  560.     if ErrCode = MAT_OK then
  561.       for I := Lbound to Ubound do
  562.         begin
  563.           DeltaX^[I] := R * DeltaX^[I];
  564.           X^[I] := X^[I] + DeltaX^[I];
  565.         end;
  566.     LinMin := ErrCode;
  567.   end;
  568.  
  569.   function ParamConv(OldX, X        : PVector;
  570.                      Lbound, Ubound : Integer;
  571.                      Tol            : Float) : Boolean;
  572. { ----------------------------------------------------------------------
  573.   Check for convergence on parameters
  574.   ---------------------------------------------------------------------- }
  575.   var
  576.     I    : Integer;
  577.     Ref  : Float;
  578.     Conv : Boolean;
  579.   begin
  580.     I := Lbound;
  581.     Conv := True;
  582.     repeat
  583.       if OldX^[I] <> 0.0 then Ref := Tol * Abs(OldX^[I]) else Ref := Tol;
  584.       Conv := Conv and (Abs(X^[I] - OldX^[I]) < Ref);
  585.       Inc(I);
  586.     until (Conv = False) or (I > Ubound);
  587.     ParamConv := Conv;
  588.   end;
  589.  
  590.   procedure CreateLogFile;
  591.   begin
  592.     Assign(LogFile, LogFileName);
  593.     Rewrite(LogFile);
  594.   end;
  595.  
  596.   function Marquardt(Func           : TFuncNVar;
  597.                      HessGrad       : THessGrad;
  598.                      X              : PVector;
  599.                      Lbound, Ubound : Integer;
  600.                      MaxIter        : Integer;
  601.                      Tol            : Float;
  602.                      var F_min      : Float;
  603.                      H_inv          : PMatrix) : Integer;
  604.   const
  605.     LAMBDA0   = 1.0E-2;   { Initial lambda value }
  606.     LAMBDAMAX = 1.0E+3;   { Highest lambda value }
  607.     FTOL      = 1.0E-10;  { Tolerance on function decrease }
  608.   var
  609.     Lambda,
  610.     Lambda1   : Float;    { Marquardt's lambda }
  611.     I         : Integer;  { Loop variable }
  612.     OldX      : PVector;  { Old parameters }
  613.     G         : PVector;  { Gradient vector }
  614.     H         : PMatrix;  { Hessian matrix }
  615.     A         : PMatrix;  { Modified Hessian matrix }
  616.     DeltaX    : PVector;  { New search direction }
  617.     F1        : Float;    { New minimum }
  618.     Lambda_Ok : Boolean;  { Successful Lambda decrease }
  619.     Conv      : Boolean;  { Convergence reached }
  620.     Done      : Boolean;  { Iterations done }
  621.     Iter      : Integer;  { Iteration count }
  622.     ErrCode   : Integer;  { Error code }
  623.   begin
  624.     if WriteLogFile then
  625.       begin
  626.         CreateLogFile;
  627.         WriteLn(LogFile, 'Marquardt');
  628.         WriteLn(LogFile, 'Iter         F            Lambda');
  629.       end;
  630.  
  631.     Lambda := LAMBDA0;
  632.     ErrCode := MAT_OK;
  633.  
  634.     DimVector(OldX, Ubound);
  635.     DimVector(G, Ubound);
  636.     DimMatrix(H, Ubound, Ubound);
  637.     DimMatrix(A, Ubound, Ubound);
  638.     DimVector(DeltaX, Ubound);
  639.  
  640.     F_min := Func(X);    { Initial function value }
  641.     LinObjFunc := Func;  { Function for line minimization }
  642.  
  643.     Iter := 1;
  644.     Conv := False;
  645.     Done := False;
  646.  
  647.     repeat
  648.       if WriteLogFile then
  649.         WriteLn(LogFile, Iter:4, '   ', F_min:12, '   ', Lambda:12);
  650.  
  651.       { Save current parameters }
  652.       CopyVector(OldX, X, Lbound, Ubound);
  653.  
  654.       { Compute Gradient and Hessian }
  655.       HessGrad(Func, X, Lbound, Ubound, G, H);
  656.       CopyMatrix(A, H, Lbound, Lbound, Ubound, Ubound);
  657.  
  658.       { Change sign of gradient }
  659.       for I := Lbound to Ubound do
  660.         G^[I] := - G^[I];
  661.  
  662.       if Conv then  { Newton-Raphson iteration }
  663.         begin
  664.           ErrCode := GaussJordan(A, G, Lbound, Ubound, H_inv, DeltaX);
  665.           if ErrCode = MAT_OK then
  666.             for I := Lbound to Ubound do
  667.               X^[I] := OldX^[I] + DeltaX^[I];
  668.           Done := True;
  669.         end
  670.       else          { Marquardt iteration }
  671.         begin
  672.           repeat
  673.             { Multiply each diagonal term of H by (1 + Lambda) }
  674.             Lambda1 := 1.0 + Lambda;
  675.             for I := Lbound to Ubound do
  676.               A^[I]^[I] := Lambda1 * H^[I]^[I];
  677.  
  678.             ErrCode := GaussJordan(A, G, Lbound, Ubound, H_inv, DeltaX);
  679.  
  680.             if ErrCode = MAT_OK then
  681.               begin
  682.                 { Initialize parameters }
  683.                 CopyVector(X, OldX, Lbound, Ubound);
  684.  
  685.                 { Minimize in the direction specified by DeltaX }
  686.                 ErrCode := LinMin(Func, X, DeltaX,
  687.                                   Lbound, Ubound, 100, 0.01, F1);
  688.  
  689.                 { Check that the function has decreased. Otherwise
  690.                   increase Lambda, without exceeding LAMBDAMAX }
  691.                 Lambda_Ok := (F1 - F_min) < F_min * FTOL;
  692.                 if not Lambda_Ok then Lambda := 10.0 * Lambda;
  693.                 if Lambda > LAMBDAMAX then ErrCode := BIG_LAMBDA;
  694.               end;
  695.           until Lambda_Ok or (ErrCode <> MAT_OK);
  696.  
  697.           { Check for convergence }
  698.           Conv := ParamConv(OldX, X, Lbound, Ubound, Tol);
  699.  
  700.           { Prepare next iteration }
  701.           Lambda := 0.1 * Lambda;
  702.           F_min := F1;
  703.         end;
  704.  
  705.       Inc(Iter);
  706.       if Iter > MaxIter then ErrCode := NON_CONV;
  707.     until Done or (ErrCode <> MAT_OK);
  708.  
  709.     DelVector(OldX, Ubound);
  710.     DelVector(G, Ubound);
  711.     DelMatrix(H, Ubound, Ubound);
  712.     DelMatrix(A, Ubound, Ubound);
  713.     DelVector(DeltaX, Ubound);
  714.  
  715.     if WriteLogFile then
  716.       Close(LogFile);
  717.  
  718.     Marquardt := ErrCode;
  719.   end;
  720.  
  721.   function BFGS(Func           : TFuncNVar;
  722.                 Gradient       : TGradient;
  723.                 X              : PVector;
  724.                 Lbound, Ubound : Integer;
  725.                 MaxIter        : Integer;
  726.                 Tol            : Float;
  727.                 var F_min      : Float;
  728.                 H_inv          : PMatrix) : Integer;
  729.   var
  730.     I, J, Iter, ErrCode            : Integer;
  731.     F1, Fae, Fad, Fac              : Float;
  732.     OldX, DeltaX, G, OldG, dG, HdG : PVector;
  733.     Conv                           : Boolean;
  734.   begin
  735.     if WriteLogFile then
  736.       begin
  737.         CreateLogFile;
  738.         WriteLn(LogFile, 'BFGS');
  739.         WriteLn(LogFile, 'Iter         F');
  740.       end;
  741.  
  742.     DimVector(OldX, Ubound);
  743.     DimVector(DeltaX, Ubound);
  744.     DimVector(G, Ubound);
  745.     DimVector(OldG, Ubound);
  746.     DimVector(dG, Ubound);
  747.     DimVector(HdG, Ubound);
  748.  
  749.     Iter := 1;
  750.     F1 := Func(X);                         { Initial function value }
  751.     Gradient(Func, X, Lbound, Ubound, G);  { Initial gradient }
  752.     LinObjFunc := Func;                    { Function for line minimization }
  753.  
  754.     for I := Lbound to Ubound do
  755.       begin
  756.         { Initialize inverse hessian to unit matrix }
  757.         for J := Lbound to Ubound do
  758.           if I = J then H_inv^[I]^[J] := 1.0 else H_inv^[I]^[J] := 0.0;
  759.         { Initial direction = opposite gradient }
  760.         DeltaX^[I] := - G^[I];
  761.       end;
  762.  
  763.     repeat
  764.       if WriteLogFile then
  765.         WriteLn(LogFile, Iter:4, '   ', F1:12);
  766.  
  767.       { Save current parameters }
  768.       CopyVector(OldX, X, Lbound, Ubound);
  769.  
  770.       { Minimize along the direction specified by DeltaX }
  771.       ErrCode := LinMin(Func, X, DeltaX, Lbound, Ubound, 100, 0.01, F_min);
  772.  
  773.       { Check for convergence on parameters }
  774.       Conv := ParamConv(OldX, X, Lbound, Ubound, Tol);
  775.  
  776.       if not Conv then
  777.         begin
  778.           { Save old function and gradient and compute new ones }
  779.           F1 := F_min;
  780.           F_min := Func(X);
  781.           CopyVector(OldG, G, Lbound, Ubound);
  782.           Gradient(Func, X, Lbound, Ubound, G);
  783.  
  784.           { Compute difference of gradients }
  785.           for I := Lbound to Ubound do
  786.             dG^[I] := G^[I] - OldG^[I];
  787.  
  788.           { Multiply by inverse hessian }
  789.           for I := Lbound to Ubound do
  790.             begin
  791.               HdG^[I] := 0.0;
  792.               for J := Lbound to Ubound do
  793.                 HdG^[I] := HdG^[I] + H_inv^[I]^[J] * dG^[J];
  794.             end;
  795.  
  796.           { Scalar products in denominator of BFGS formula }
  797.           Fac := 0.0; Fae := 0.0;
  798.           for I := Lbound to Ubound do
  799.             begin
  800.               Fac := Fac + dG^[I] * DeltaX^[I];
  801.               Fae := Fae + dG^[I] * HdG^[I];
  802.             end;
  803.  
  804.           { Apply BFGS correction }
  805.           Fac := 1.0 / Fac;
  806.           Fad := 1.0 / Fae;
  807.           for I := Lbound to Ubound do
  808.             dG^[I] := Fac * DeltaX^[I] - Fad * HdG^[I];
  809.           for I := Lbound to Ubound do
  810.             for J := Lbound to Ubound do
  811.               H_inv^[I]^[J] := H_inv^[I]^[J] + Fac * DeltaX^[I] * DeltaX^[J]
  812.                                              - Fad * HdG^[I] * HdG^[J]
  813.                                              + Fae * dG^[I] * dG^[J];
  814.  
  815.           { Compute new search direction }
  816.           for I := Lbound to Ubound do
  817.             begin
  818.               DeltaX^[I] := 0.0;
  819.               for J := Lbound to Ubound do
  820.                 DeltaX^[I] := DeltaX^[I] - H_inv^[I]^[J] * G^[J];
  821.             end
  822.         end;
  823.  
  824.       Inc(Iter);
  825.     until Conv or (Iter > MaxIter);
  826.  
  827.     DelVector(OldX, Ubound);
  828.     DelVector(DeltaX, Ubound);
  829.     DelVector(G, Ubound);
  830.     DelVector(OldG, Ubound);
  831.     DelVector(dG, Ubound);
  832.     DelVector(HdG, Ubound);
  833.  
  834.     if WriteLogFile then
  835.       Close(LogFile);
  836.  
  837.     if Iter > MaxIter then
  838.       BFGS := NON_CONV
  839.     else
  840.       BFGS := MAT_OK;
  841.   end;
  842.  
  843.   function Simplex(Func           : TFuncNVar;
  844.                    X              : PVector;
  845.                    Lbound, Ubound : Integer;
  846.                    MaxIter        : Integer;
  847.                    Tol            : Float;
  848.                    var F_min      : Float) : Integer;
  849.   const
  850.     STEP = 1.50;  { Step used to construct the initial simplex }
  851.   var
  852.     P             : PMatrix;  { Simplex coordinates }
  853.     F             : PVector;  { Function values }
  854.     Pbar          : PVector;  { Centroid coordinates }
  855.     Pstar, P2star : PVector;  { New vertices }
  856.     Ystar, Y2star : Float;    { New function values }
  857.     F0            : Float;    { Function value at minimum }
  858.     N             : Integer;  { Number of parameters }
  859.     M             : Integer;  { Index of last vertex }
  860.     L, H          : Integer;  { Vertices with lowest & highest F values }
  861.     I, J          : Integer;  { Loop variables }
  862.     Iter          : Integer;  { Iteration count }
  863.     Corr, MaxCorr : Float;    { Corrections }
  864.     Sum           : Float;
  865.     Flag          : Boolean;
  866.  
  867.     procedure UpdateSimplex(Y : Float; Q : PVector);
  868.     { Update "worst" vertex and function value }
  869.     begin
  870.       F^[H] := Y;
  871.       CopyVector(P^[H], Q, Lbound, Ubound);
  872.     end;
  873.  
  874.   begin
  875.     if WriteLogFile then
  876.       begin
  877.         CreateLogFile;
  878.         WriteLn(LogFile, 'Simplex');
  879.         WriteLn(LogFile, 'Iter         F');
  880.       end;
  881.  
  882.     N := Ubound - Lbound + 1;
  883.     M := Succ(Ubound);
  884.  
  885.     DimMatrix(P, M, Ubound);
  886.     DimVector(F, M);
  887.     DimVector(Pbar, Ubound);
  888.     DimVector(Pstar, Ubound);
  889.     DimVector(P2star, Ubound);
  890.  
  891.     Iter := 1;
  892.     F0 := MAXNUM;
  893.  
  894.     { Construct initial simplex }
  895.     for I := Lbound to M do
  896.       CopyVector(P^[I], X, Lbound, Ubound);
  897.     for I := Lbound to Ubound do
  898.       P^[I]^[I] := P^[I]^[I] * STEP;
  899.  
  900.     { Evaluate function at each vertex }
  901.     for I := Lbound to M do
  902.       F^[I] := Func(P^[I]);
  903.  
  904.     repeat
  905.       { Find vertices (L,H) having the lowest and highest
  906.         function values, i.e. "best" and "worst" vertices }
  907.       L := Lbound;
  908.       H := Lbound;
  909.       for I := Succ(Lbound) to M do
  910.         if F^[I] < F^[L] then
  911.           L := I
  912.         else if F^[I] > F^[H] then
  913.           H := I;
  914.       if F^[L] < F0 then
  915.         F0 := F^[L];
  916.  
  917.       if WriteLogFile then
  918.         WriteLn(LogFile, Iter:4, '   ', F0:12);
  919.  
  920.       { Find centroid of points other than P(H) }
  921.       for J := Lbound to Ubound do
  922.         begin
  923.           Sum := 0.0;
  924.           for I := Lbound to M do
  925.             if I <> H then Sum := Sum + P^[I]^[J];
  926.           Pbar^[J] := Sum / N;
  927.         end;
  928.  
  929.       { Reflect worst vertex through centroid }
  930.       for J := Lbound to Ubound do
  931.         Pstar^[J] := 2.0 * Pbar^[J] - P^[H]^[J];
  932.       Ystar := Func(Pstar);
  933.  
  934.       { If reflection successful, try extension }
  935.       if Ystar < F^[L] then
  936.         begin
  937.           for J := Lbound to Ubound do
  938.             P2star^[J] := 3.0 * Pstar^[J] - 2.0 * Pbar^[J];
  939.           Y2star := Func(P2star);
  940.  
  941.           { Retain extension or contraction }
  942.           if Y2star < F^[L] then
  943.             UpdateSimplex(Y2star, P2star)
  944.           else
  945.             UpdateSimplex(Ystar, Pstar);
  946.         end
  947.       else
  948.         begin
  949.           I := Lbound;
  950.           Flag := False;
  951.           repeat
  952.             if (I <> H) and (F^[I] > Ystar) then Flag := True;
  953.             Inc(I);
  954.           until Flag or (I > M);
  955.           if Flag then
  956.             UpdateSimplex(Ystar, Pstar)
  957.           else
  958.             begin
  959.               { Contraction on the reflection side of the centroid }
  960.               if Ystar <= F^[H] then
  961.                 UpdateSimplex(Ystar, Pstar);
  962.  
  963.               { Contraction on the opposite side of the centroid }
  964.               for J := Lbound to Ubound do
  965.                 P2star^[J] := 0.5 * (P^[H]^[J] + Pbar^[J]);
  966.               Y2star := Func(P2star);
  967.               if Y2star <= F^[H] then
  968.                 UpdateSimplex(Y2star, P2star)
  969.               else
  970.                 { Contract whole simplex }
  971.                 for I := Lbound to M do
  972.                   for J := Lbound to Ubound do
  973.                     P^[I]^[J] := 0.5 * (P^[I]^[J] + P^[L]^[J]);
  974.             end;
  975.         end;
  976.  
  977.       { Test convergence }
  978.       MaxCorr := 0.0;
  979.       for J := Lbound to Ubound do
  980.         begin
  981.           Corr := Abs(P^[H]^[J] - P^[L]^[J]);
  982.           if Corr > MaxCorr then MaxCorr := Corr;
  983.         end;
  984.       Inc(Iter);
  985.     until (MaxCorr < Tol) or (Iter > MaxIter);
  986.  
  987.     CopyVector(X, P^[L], Lbound, Ubound);
  988.     F_min := F^[L];
  989.  
  990.     DelMatrix(P, M, Ubound);
  991.     DelVector(F, M);
  992.     DelVector(Pbar, Ubound);
  993.     DelVector(Pstar, Ubound);
  994.     DelVector(P2star, Ubound);
  995.  
  996.     if WriteLogFile then
  997.       Close(LogFile);
  998.  
  999.     if Iter > MaxIter then
  1000.       Simplex := NON_CONV
  1001.     else
  1002.       Simplex := MAT_OK;
  1003.   end;
  1004.  
  1005.   procedure GenRandomVector(Cv, MinSigma   : Float;
  1006.                             X, X1          : PVector;
  1007.                             Lbound, Ubound : Integer);
  1008. { ----------------------------------------------------------------------
  1009.   Generates a random vector X1 from X
  1010.   ---------------------------------------------------------------------- }
  1011.   var
  1012.     I     : Integer;
  1013.     Sigma : Float;   { S.D. of gaussian random numbers }
  1014.   begin
  1015.     for I := Lbound to Ubound do
  1016.       begin
  1017.         Sigma := (Cv * Abs(X^[I]) + MinSigma);
  1018.         X1^[I] := RanGauss(X^[I], Sigma);
  1019.       end;
  1020.   end;
  1021.  
  1022.   function InitTemp(Func           : TFuncNVar;
  1023.                     X              : PVector;
  1024.                     Lbound, Ubound : Integer;
  1025.                     var F          : Float) : Float;
  1026. { ----------------------------------------------------------------------
  1027.   Computes the initial temperature so that the probability
  1028.   of accepting an increase of the function is about 0.5
  1029.   ---------------------------------------------------------------------- }
  1030.   const
  1031.     N_MAX = 50;  { Max number of function evaluations }
  1032.     N_INC = 10;  { Max number of function increases }
  1033.   var
  1034.     T     : Float;    { Temperature }
  1035.     X1    : PVector;  { New coordinates }
  1036.     F1    : Float;    { New function value }
  1037.     F_inc : PVector;  { Function increases }
  1038.     N     : Integer;  { Counts function evaluations }
  1039.     K     : Integer;  { Counts function increases }
  1040.   begin
  1041.     DimVector(X1, Ubound);
  1042.     DimVector(F_inc, N_INC);
  1043.  
  1044.     K := 0;
  1045.     N := 0;
  1046.     T := 0.0;
  1047.     F := Func(X);
  1048.  
  1049.     { Compute K (<= N_INC) increases of the function }
  1050.     repeat
  1051.       Inc(N);
  1052.       GenRandomVector(SA_Cv, SA_MinSigma, X, X1, Lbound, Ubound);
  1053.       F1 := Func(X1);
  1054.       if F1 > F then
  1055.         begin
  1056.           Inc(K);
  1057.           F_inc^[K] := F1 - F;
  1058.         end;
  1059.     until (K = N_INC) or (N = N_MAX);
  1060.  
  1061.     { The median M of these K increases has a probability of 1/2.
  1062.       From Boltzmann's formula: Exp(-M/T) = 1/2 ==> T = M / Ln(2) }
  1063.     T := Median(F_inc, 1, K) / LN2;
  1064.     if T = 0.0 then T := 1.0;
  1065.     InitTemp := T;
  1066.  
  1067.     DelVector(X1, Ubound);
  1068.     DelVector(F_inc, N_INC);
  1069.   end;
  1070.  
  1071.   function Accept(DeltaF, T        : Float;
  1072.                   var N_inc, N_acc : Integer) : Boolean;
  1073. { ----------------------------------------------------------------------
  1074.   Checks if a variation DeltaF of the function at temperature T is
  1075.   acceptable. Updates the counters N_inc (number of increases of the
  1076.   function) and N_acc (number of accepted increases).
  1077.   ---------------------------------------------------------------------- }
  1078.   begin
  1079.     if DeltaF < 0.0 then
  1080.       Accept := True
  1081.     else
  1082.       begin
  1083.         Inc(N_inc);
  1084.         if Exp(- DeltaF / T) > RanMar then
  1085.           begin
  1086.             Accept := True;
  1087.             Inc(N_acc);
  1088.           end
  1089.         else
  1090.           Accept := False;
  1091.       end;
  1092.   end;
  1093.  
  1094.   procedure SimAnnCycle(Func           : TFuncNVar;
  1095.                         X              : PVector;
  1096.                         Lbound, Ubound : Integer;
  1097.                         MaxIter        : Integer;
  1098.                         T_Ratio        : Float;
  1099.                         var LogFile    : Text;
  1100.                         var F_min      : Float);
  1101. { ----------------------------------------------------------------------
  1102.   Performs one cycle of simulated annealing
  1103.   ---------------------------------------------------------------------- }
  1104.   var
  1105.     T0, T    : Float;    { Temperature }
  1106.     Rt       : Float;    { Temperature reduction factor }
  1107.     Rs       : Float;    { S.D. reduction factor }
  1108.     MinSigma : Float;    { Minimal S.D. }
  1109.     Cv       : Float;    { Coef. of variation }
  1110.     F, F1    : Float;    { Function values }
  1111.     DeltaF   : Float;    { Variation of function }
  1112.     X1       : PVector;  { New coordinates }
  1113.     K        : Integer;  { Loop variable }
  1114.     Iter     : Integer;  { Iteration count }
  1115.     N_inc    : Integer;  { Nb of function increases }
  1116.     N_acc    : Integer;  { Nb of accepted increases }
  1117.   begin
  1118.     { Dimension array }
  1119.     DimVector(X1, Ubound);
  1120.  
  1121.     { Compute initial temperature and reduction factors }
  1122.     T0 := InitTemp(Func, X, Lbound, Ubound, F);
  1123.     Rt := Exp(Ln(T_Ratio) / MaxIter);
  1124.     Rs := Exp(Ln(SA_SigmaFact) / MaxIter);
  1125.  
  1126.     T := T0;
  1127.     Iter := 0;
  1128.     Cv := SA_Cv;
  1129.     MinSigma := SA_MinSigma;
  1130.  
  1131.     repeat
  1132.       N_inc := 0;
  1133.       N_acc := 0;
  1134.  
  1135.       { Perform SA_Nt evaluations at constant temperature }
  1136.       for K := 1 to SA_Nt do
  1137.         begin
  1138.           { Generate new point }
  1139.           GenRandomVector(Cv, MinSigma, X, X1, Lbound, Ubound);
  1140.  
  1141.           { Compute new function value }
  1142.           F1 := Func(X1);
  1143.           DeltaF := F1 - F;
  1144.  
  1145.           { Check for acceptance }
  1146.           if Accept(DeltaF, T, N_inc, N_acc) then
  1147.             begin
  1148.               CopyVector(X, X1, Lbound, Ubound);
  1149.               F := F1;
  1150.             end;
  1151.         end;
  1152.  
  1153.       if WriteLogFile then
  1154.         WriteLn(LogFile, Iter:4, '   ', T:12, '   ', F:12, N_inc:6, N_acc:6);
  1155.  
  1156.       { Update temperature, S.D. factors, and iteration count }
  1157.       T := T * Rt;
  1158.       Cv := Cv * Rs;
  1159.       MinSigma := MinSigma * Rs;
  1160.  
  1161.       Inc(Iter);
  1162.     until Iter > MaxIter;
  1163.  
  1164.     F_min := F;
  1165.     DelVector(X1, Ubound);
  1166.   end;
  1167.  
  1168.   function SimAnn(Func           : TFuncNVar;
  1169.                   X              : PVector;
  1170.                   Lbound, Ubound : Integer;
  1171.                   MaxIter        : Integer;
  1172.                   T_Ratio        : Float;
  1173.                   var F_min      : Float) : Integer;
  1174.   var
  1175.     K : Integer;
  1176.   begin
  1177.     if WriteLogFile then
  1178.       CreateLogFile;
  1179.  
  1180.     { Initialize the Marsaglia random number generator
  1181.       using the standard Pascal generator }
  1182.     Randomize;
  1183.     RMarIn(System.Random(10000), System.Random(10000));
  1184.  
  1185.     for K := 1 to SA_NCycles do
  1186.       begin
  1187.         if WriteLogFile then
  1188.           WriteLn(LogFile, 'Simulated annealing: Cycle ', K, #10,
  1189.                     'Iter         T              F        Inc   Acc');
  1190.         SimAnnCycle(Func, X, Lbound, Ubound, MaxIter, T_Ratio, LogFile, F_min);
  1191.       end;
  1192.  
  1193.     if WriteLogFile then
  1194.       Close(LogFile);
  1195.  
  1196.     SimAnn := 0;
  1197.   end;
  1198.  
  1199. begin
  1200.   X1 := nil;
  1201.   DeltaX1 := nil;
  1202.   Ubound1 := 1;
  1203.   Eps := Pow(MACHEP, 0.333);
  1204. end.

Raw Paste


Login or Register to edit or fork this paste. It's free.