' ******************************************************************
' Minimization of a function of one variable by Golden Search method
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE OptOk 0  
' No error

#DEFINE OptNonConv -1  
' Non-convergence

' ------------------------------------------------------------------
' Mathematical constants
' ------------------------------------------------------------------

#DEFINE Gold 1.618033988749895#   
' Golden Mean (1 + Sqrt(5))/2

#DEFINE CGold 0.3819660112501052#  
' 2 - Gold

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' **********************************************************************

SUB MinBrack (Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
              BYREF A AS DOUBLE, BYREF B AS DOUBLE, BYREF C AS DOUBLE,_
              BYREF Fa AS DOUBLE, BYREF Fb AS DOUBLE, BYREF Fc AS DOUBLE)
' ----------------------------------------------------------------------
' Given two points (A, B) this procedure finds a triplet (A, B, C)
' such that:
'
' 1) A < B < C
' 2) A, B, C are within the golden ratio
' 3) Func(B) < Func(A) and Func(B) < Func(C).
'
' The corresponding function values are returned in Fa, Fb, Fc
' ----------------------------------------------------------------------

  IF A > B THEN SWAP A, B

  Fa = Func(A)
  Fb = Func(B)
  IF Fb > Fa THEN
    SWAP A, B
    SWAP Fa, Fb
  END IF

  C = B + Gold * (B - A)
  Fc = Func(C)
  DO WHILE Fc < Fb
    A = B
    B = C
    Fa = Fb
    Fb = Fc
    C = B + Gold * (B - A)
    Fc = Func(C)
  LOOP

  IF A > C THEN
    SWAP A, C
    SWAP Fa, Fc
  END IF
END SUB

SUB GoldSearch(Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
               A AS DOUBLE, B AS DOUBLE, MaxIter AS INTEGER,_
               Tol AS DOUBLE, BYREF Xmin AS DOUBLE, BYREF Ymin AS DOUBLE)

  CONST MinTol = 1.5D-8  ' Minimal tolerance ~ sqrt(MachEp)

  DIM AS DOUBLE  C, Fa, Fb, Fc, X0, X1, X2, X3, F1, F2
  DIM AS INTEGER Iter 

  IF Tol < MinTol THEN Tol = MinTol

  MinBrack Func, A, B, C, Fa, Fb, Fc

  X0 = A
  X3 = C
  IF (C - B) > (B - A) THEN
    X1 = B
    X2 = B + CGold * (C - B)
    F1 = Fb
    F2 = Func(X2)
  ELSE
    X1 = B - CGold * (B - A)
    X2 = B
    F1 = Func(X1)
    F2 = Fb
  END IF

  Iter = 0

  DO WHILE (Iter <= MaxIter) AND (ABS(X3 - X0) > Tol * (ABS(X1) + ABS(X2)))
    IF F2 < F1 THEN
      X0 = X1
      X1 = X2
      F1 = F2
      X2 = X1 + CGold * (X3 - X1)
      F2 = Func(X2)
    ELSE
      X3 = X2
      X2 = X1
      F2 = F1
      X1 = X2 - CGold * (X2 - X0)
      F1 = Func(X1)
    END IF
    Iter = Iter + 1
  LOOP

  IF F1 < F2 THEN
    Xmin = X1
    Ymin = F1
  ELSE
    Xmin = X2
    Ymin = F2
  END IF

  IF Iter > MaxIter THEN
    ErrCode = OptNonConv
  ELSE
    ErrCode = OptOk
  END IF
END SUB

