' ******************************************************************
' Minimization of a function of several variables by Newton-Raphson
' 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     )
'
' Det(H) = 80000 * (X^2 - Y) + 400
'
' 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 and hessian of
' the function. It is recommended to use analytical derivatives
' whenever possible. Otherwise you can use the alternative code
' provided below.
' ******************************************************************

SUB HessGrad (X() AS DOUBLE, G() AS DOUBLE, H() AS DOUBLE)

  DIM AS DOUBLE A, B, C

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

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

  H(1, 1) = 1200 * C - 400 * X(2) + 2
  H(1, 2) = -400 * X(1)
  H(2, 1) = H(1, 2)
  H(2, 2) = 200
END SUB

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

'SUB HessGrad (X() AS DOUBLE, G() AS DOUBLE, H() AS DOUBLE)
'
' Define the relative increment to compute numerical derivatives
' (should not be lower than MachEp^1/3 which is about 6.06D-6
'
'  CONST Eta = 1D-6
'
'  NumHessGrad @Func, X(), Eta, G(), H()
'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 DOUBLE  Det                         ' Determinant of hessian
DIM AS INTEGER I, J                        ' Loop variables

Newton @Func, @HessGrad, X(), MaxIter, Tol, Fmin, G(), Hinv(), Det

SELECT CASE MathErr()
  CASE OptNonConv
    PRINT "Non-convergence!"
    STOP
  CASE OptSing
    PRINT "Singular Hessian matrix!"
    STOP
END SELECT

CLS

PRINT "Minimization of Rosenbrock's function (Newton-Raphson 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

PRINT : PRINT "Determinant of hessian:": PRINT

PRINT USING "#####.######"; Det

SLEEP
