' ******************************************************************
' Minimization of a function of several variables by the Broyden-
' Fletcher-Goldfarb-Shanno (BFGS) method.
'
' Example: Rosenbrock's function:
'
' F(X, Y) = 100 * (Y - X^2)^2 + (1 - X)^2
'
'               ( -400 * (Y - X^2) * X - 2 + 2 * X )
' Gradient: G = (                                  )
'               (  200 * Y - 200 * X^2             )
'
'               ( 1200 * X^2 - 400 * Y + 2     -400 * X )
' Hessian:  H = (                                       )
'               ( -400 * X                      200     )
'
' True minimum is at (1, 1), F = 0
'
' The inverse hessian at the minimum is:
'
'   ( 1/2     1       )
'   ( 1       401/200 )
'
' Ref: H. Rosenbrock, Comput. J., 1960, 3, 175
' ******************************************************************

#INCLUDE "fbmath.bi"

' ******************************************************************
' Define number of variables, number of iterations, and precision
' ******************************************************************

CONST Nvar    = 2     ' Number of variables
CONST MaxIter = 1000  ' Max number of iterations
CONST Tol     = 1D-6  ' Required precision (must be > Sqrt(MachEp))

' ******************************************************************
' Define starting point
' ******************************************************************

' Variables: X(1) = X, X(2) = Y
DIM AS DOUBLE X(1 TO Nvar) = {2, 2}

' ******************************************************************
' Define the function to be minimized
' ******************************************************************

FUNCTION Func(X() AS DOUBLE) AS DOUBLE

  DIM AS DOUBLE A, B

  A = X(2) - X(1) * X(1)
  B = 1 - X(1)

  Func = 100 * A * A + B * B
END FUNCTION

' ******************************************************************
' Define the subroutine which computes the gradient of the function.
' It is recommended to use analytical derivatives whenever possible.
' Otherwise you can use the alternative code provided below.
' ******************************************************************

SUB Gradient(X() AS DOUBLE, G() AS DOUBLE)

  DIM AS DOUBLE A, B

  A = X(2) - X(1) * X(1)
  B = 1 - X(1)

  G(1) = -400 * X(1) * A - 2 * B
  G(2) = 200 * A
END SUB

' ******************************************************************
' Alternative code to be used if the analytical derivatives are not
' available.
' ******************************************************************

' SUB Gradient (X() AS DOUBLE, G() AS DOUBLE)
'
'   Define the relative increment to compute numerical derivatives
'   (should not be lower than MachEp^1/2 which is about 1.5D-8)
'
'   CONST Eta = 1D-6
'
'   NumGrad @Func, X(), Eta, G()
' END SUB

' ******************************************************************
' Main program
' ******************************************************************

DIM AS DOUBLE  G(1 TO Nvar)                ' Gradient
DIM AS DOUBLE  Hinv(1 TO Nvar, 1 TO Nvar)  ' Inverse hessian matrix
DIM AS DOUBLE  Fmin                        ' Function value at min.
DIM AS INTEGER I, J                        ' Loop variables

BFGS @Func, @Gradient, X(), MaxIter, Tol, Fmin, G(), Hinv()

IF MathErr = OptNonConv THEN
  PRINT "Non-convergence!"
  STOP
END IF

CLS

PRINT "Minimization of Rosenbrock's function (BFGS method)"
PRINT "---------------------------------------------------"
PRINT

PRINT "Coordinates of minimum:": PRINT

FOR I = 1 TO Nvar
  PRINT USING "X(#) = #####.######"; I; X(I)
NEXT I

PRINT : PRINT "Function value:" : PRINT

PRINT "Fmin  = "; Fmin

PRINT : PRINT "Gradient:" : PRINT

FOR I = 1 TO Nvar
  PRINT USING "#####.######"; G(I)
NEXT I

PRINT : PRINT "Inverse hessian matrix:": PRINT

FOR I = 1 TO Nvar
  FOR J = 1 TO Nvar
    PRINT USING "#####.######"; Hinv(I, J);
  NEXT J
  PRINT
NEXT I

SLEEP
