F1 := Func(X1);
          Inc(Iter);
        end;
    if F1 < F2 then
      begin
        Xmin := X1;
        Ymin := F1;
      end
    else
      begin
        Xmin := X2;
        Ymin := F2;
      end;
    if Iter > MaxIter then
      GoldSearch := OPT_NON_CONV
    else
      GoldSearch := OPT_OK;
  end;

  procedure CreateLogFile;
  begin
    Assign(LogFile, LogFileName);
    Rewrite(LogFile);
  end;

  function Simplex(Func           : TFuncNVar;
                   X              : TVector;
                   Lbound, Ubound : Integer;
                   MaxIter        : Integer;
                   Tol            : Float;
                   var F_min      : Float) : Integer;
  const
    STEP = 1.50;  { Step used to construct the initial simplex }
  var
    P             : TMatrix;  { Simplex coordinates }
    F             : TVector;  { Function values }
    Pbar          : TVector;  { Centroid coordinates }
    Pstar, P2star : TVector;  { New vertices }
    Ystar, Y2star : Float;    { New function values }
    F0            : Float;    { Function value at minimum }
    N             : Integer;  { Number of parameters }
    M             : Integer;  { Index of last vertex }
    L, H          : Integer;  { Vertices with lowest & highest F values }
    I, J          : Integer;  { Loop variables }
    Iter          : Integer;  { Iteration count }
    Corr, MaxCorr : Float;    { Corrections }
    Sum           : Float;
    Flag          : Boolean;

    procedure UpdateSimplex(Y : Float; Q : TVector);
    { Update "worst" vertex and function value }
    begin
      F[H] := Y;
      CopyVector(P[H], Q, Lbound, Ubound);
    end;

  begin
    if WriteLogFile then
      begin
        CreateLogFile;
        WriteLn(LogFile, 'Simplex');
        WriteLn(LogFile, 'Iter         F');
      end;

    N := Ubound - Lbound + 1;
    M := Succ(Ubound);

    DimMatrix(P, M, Ubound);
    DimVector(F, M);
    DimVector(Pbar, Ubound);
    DimVector(Pstar, Ubound);
    DimVector(P2star, Ubound);

    Iter := 1;
    F0 := MAXNUM;

    { Construct initial simplex }
    for I := Lbound to M do
      CopyVector(P[I], X, Lbound, Ubound);
    for I := Lbound to Ubound do
      P[I,I] := P[I,I] * STEP;

    { Evaluate function at each vertex }
    for I := Lbound to M do
      F[I] := Func(P[I]);

    repeat
      { Find vertices (L,H) having the lowest and highest
        function values, i.e. "best" and "worst" vertices }
      L := Lbound;
      H := Lbound;
      for I := Succ(Lbound) to M do
        if F[I] < F[L] then
          L := I
        else if F[I] > F[H] then
          H := I;
      if F[L] < F0 then
        F0 := F[L];

      if WriteLogFile then
        WriteLn(LogFile, Iter:4, '   ', F0:12);

      { Find centroid of points other than P(H) }
      for J := Lbound to Ubound do
        begin
          Sum := 0.0;
          for I := Lbound to M do
            if I <> H then Sum := Sum + P[I,J];
          Pbar[J] := Sum / N;
        end;

      { Reflect worst vertex through centroid }
      for J := Lbound to Ubound do
        Pstar[J] := 2.0 * Pbar[J] - P[H,J];
      Ystar := Func(Pstar);

      { If reflection successful, try extension }
      if Ystar < F[L] then
        begin
          for J := Lbound to Ubound do
            P2star[J] := 3.0 * Pstar[J] - 2.0 * Pbar[J];
          Y2star := Func(P2star);

          { Retain extension or contraction }
          if Y2star < F[L] then
            UpdateSimplex(Y2star, P2star)
          else
            UpdateSimplex(Ystar, Pstar);
        end
      else
        begin
          I := Lbound;
          Flag := False;
          repeat
            if (I <> H) and (F[I] > Ystar) then Flag := True;
            Inc(I);
          until Flag or (I > M);
          if Flag then
            UpdateSimplex(Ystar, Pstar)
          else
            begin
              { Contraction on the reflection side of the centroid }
              if Ystar <= F[H] then
                UpdateSimplex(Ystar, Pstar);

              { Contraction on the opposite side of the centroid }
              for J := Lbound to Ubound do
                P2star[J] := 0.5 * (P[H,J] + Pbar[J]);
              Y2star := Func(P2star);
              if Y2star <= F[H] then
                UpdateSimplex(Y2star, P2star)
              else
                { Contract whole simplex }
                for I := Lbound to M do
                  for J := Lbound to Ubound do
                    P[I,J] := 0.5 * (P[I,J] + P[L,J]);
            end;
        end;

      { Test convergence }
      MaxCorr := 0.0;
      for J := Lbound to Ubound do
        begin
          Corr := Abs(P[H,J] - P[L,J]);
          if Corr > MaxCorr then MaxCorr := Corr;
        end;
      Inc(Iter);
    until (MaxCorr < Tol) or (Iter > MaxIter);

    CopyVector(X, P[L], Lbound, Ubound);
    F_min := F[L];

    if WriteLogFile then
      Close(LogFile);

    if Iter > MaxIter then
      Simplex := OPT_NON_CONV
    else
      Simplex := OPT_OK;
  end;

  function F1dim(R : Float) : Float;
{ ----------------------------------------------------------------------
  Function used by LinMin to find the minimum of the objective function
  LinObjFunc in the direction specified by the global variables X1 and
  DeltaX1. R is the step in this direction.
  ---------------------------------------------------------------------- }
  var
    I : Integer;
  begin
    for I := Lbound1 to Ubound1 do
      Xt[I] := X1[I] + R * DeltaX1[I];
    F1dim := LinObjFunc(Xt);
  end;

  function LinMin(Func           : TFuncNVar;
                  X, DeltaX      : TVector;
                  Lbound, Ubound : Integer;
                  MaxIter        : Integer;
                  Tol            : Float;
                  var F_min      : Float) : Integer;
  var
    I, ErrCode : Integer;
    R : Float;
  begin
    { Redimension global vectors }
    DimVector(X1, Ubound);
    DimVector(DeltaX1, Ubound);
    DimVector(Xt, Ubound);

    Lbound1 := Lbound;
    Ubound1 := Ubound;

    { Initialize global variables }
    LinObjFunc := Func;
    for I := Lbound to Ubound do
      begin
        X1[I] := X[I];
        DeltaX1[I] := DeltaX[I]
      end;

    { Perform golden search }
    ErrCode := GoldSearch(F1dim, 0.0, 1.0, MaxIter, Tol, R, F_min);

    { Update variables }
    if ErrCode = OPT_OK then
      for I := Lbound to Ubound do
        X[I] := X[I] + R * DeltaX[I];

    LinMin := ErrCode;
  end;

  procedure NumGradient(Func           : TFuncNVar;
                        X              : TVector;
                        Lbound, Ubound : Integer;
                        G              : TVector);
  var
    Temp, Delta, Fplus, Fminus : Float;
    I                          : Integer;
  begin
    for I := Lbound to Ubound do
      begin
        Temp := X[I];
        if Temp <> 0.0 then Delta := Eps * Abs(Temp) else Delta := Eps;
        X[I] := Temp - Delta;
        Fminus := Func(X);
        X[I] := Temp + Delta;
        Fplus := Func(X);
        G[I] := (Fplus - Fminus) / (2.0 * Delta);
        X[I] := Temp;
      end;
  end;

  procedure NumHessGrad(Func           : TFuncNVar;
                        X              : TVector;
                        Lbound, Ubound : Integer;
                        G              : TVector;
                        H              : TMatrix);
  var
    Delta, Xminus, Xplus, Fminus, Fplus : TVector;
    Temp1, Temp2, F, F2plus             : Float;
    I, J                                : Integer;
  begin
    DimVector(Delta, Ubound);   { Increments   }
    DimVector(Xminus, Ubound);  { X - Delta    }
    DimVector(Xplus, Ubound);   { X + Delta    }
    DimVector(Fminus, Ubound);  { F(X - Delta) }
    DimVector(Fplus, Ubound);   { F(X + Delta) }

    F := Func(X);

    for I := Lbound to Ubound do
      begin
        if X[I] <> 0.0 then
          Delta[I] := Eps * Abs(X[I])
        else
          Delta[I] := Eps;
        Xplus[I] := X[I] + Delta[I];
        Xminus[I] := X[I] - Delta[I];
      end;

    for I := Lbound to Ubound do
      begin
        Temp1 := X[I];
        X[I] := Xminus[I];
        Fminus[I] := Func(X);
        X[I] := Xplus[I];
        Fplus[I] := Func(X);
        X[I] := Temp1;
      end;

    for I := Lbound to Ubound do
      begin
        G[I] := (Fplus[I] - Fminus[I]) / (2.0 * Delta[I]);
        H[I,I] := (Fplus[I] + Fminus[I] - 2.0 * F) / Sqr(Delta[I]);
      end;

    for I := Lbound to Pred(Ubound) do
      begin
        Temp1 := X[I];
        X[I] := Xplus[I];
        for J := Succ(I) to Ubound do
          begin
            Temp2 := X[J];
            X[J] := Xplus[J];
            F2plus := Func(X);
            H[I,J] := (F2plus - Fplus[I] - Fplus[J] + F) / (Delta[I] * Delta[J]);
            H[J,I] := H[I,J];
            X[J] := Temp2;