' ******************************************************************
' Numerical Hessian and Gradient
'
' 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     )
'
' Ref: H. Rosenbrock, Comput. J., 1960, 3, 175
' ******************************************************************

#INCLUDE "fbmath.bi"

CONST Nvar = 2     ' Number of variables
CONST Eta  = 1D-6  ' Relative increment (must be > 6.06D-6)

DIM AS INTEGER I, J
DIM AS DOUBLE  F

DIM AS DOUBLE X(1 TO Nvar) = {1, 2}

DIM AS DOUBLE G(1 TO Nvar)
DIM AS DOUBLE H(1 TO Nvar, 1 TO Nvar)

FUNCTION Func(X() AS DOUBLE) AS DOUBLE
' Rosenbrock's function

  DIM A AS DOUBLE, B AS DOUBLE

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

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

F = Func(X())

NumHessGrad @Func, X(), Eta, G(), H()

CLS
PRINT "Function :"
PRINT
PRINT USING "#####.######"; F

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

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

SLEEP
