' ******************************************************************
' Minimization of a function of several variables along a line
' ******************************************************************

' ------------------------------------------------------------------
' 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 LinMin (Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
            X() AS DOUBLE, DeltaX() AS DOUBLE, BYREF R AS DOUBLE,_ 
            MaxIter AS INTEGER, Tol AS DOUBLE, BYREF Fmin AS DOUBLE)
' ------------------------------------------------------------------
' Minimizes function Func from point X in the direction specified by
' DeltaX, using golden search as the minimization algorithm
' ------------------------------------------------------------------
' Input parameters  : X       = starting point
'                     DeltaX  = search direction
'                     R       = initial step, in fraction of |DeltaX|
'                     MaxIter = maximum number of iterations
'                     Tol     = required precision
' ----------------------------------------------------------------------
' Output parameters : X    = refined minimum coordinates
'                     Fmin = function value at minimum
'                     R    = step corresponding to the minimum
' ----------------------------------------------------------------------
' Possible results  : OptOk      = no error
'                     OptNonConv = non-convergence
' ----------------------------------------------------------------------

  CONST MinTol = 1.5D-8  ' Minimal tolerance (~ Sqrt(MachEp))
 
  DIM AS INTEGER I, L, N, Iter
  DIM AS DOUBLE  A, B, C, Fa, Fb, Fc, R0, R1, R2, R3, F1, F2, Norm
  
  L = LBOUND(X)
  N = UBOUND(X)

  DIM P(L TO N) AS DOUBLE
 
  IF R <= 0 THEN R = 1

  Norm = 0
  FOR I = L TO N
    Norm = Norm + DeltaX(I) * DeltaX(I)
  NEXT I
  Norm = SQR(Norm)

  ' Bracket the minimum (see procedure MinBrack in module GOLDSRCH)

  A = 0: B = R * Norm
 
  FOR I = L TO N
    P(I) = X(I)
  NEXT I
 
  Fa = Func(P())

  FOR I = L TO N
    P(I) = X(I) + B * DeltaX(I)
  NEXT I

  Fb = Func(P())

  IF Fb > Fa THEN
    SWAP A, B
    SWAP Fa, Fb
  END IF

  C = B + Gold * (B - A)

  FOR I = L TO N
    P(I) = X(I) + C * DeltaX(I)
  NEXT I

  Fc = Func(P())

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

    FOR I = L TO N
      P(I) = X(I) + C * DeltaX(I)
    NEXT I

    Fc = Func(P())
  LOOP

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

  ' Refine the minimum (see procedure GoldSearch in module GOLDSRCH)

  IF Tol < MinTol THEN Tol = MinTol
 
  R0 = A: R3 = C

  IF (C - B) > (B - A) THEN
    R1 = B
    R2 = B + CGold * (C - B)
    F1 = Fb
  
    FOR I = L TO N
      P(I) = X(I) + R2 * DeltaX(I)
    NEXT I
  
    F2 = Func(P())
  ELSE
    R1 = B - CGold * (B - A)
    R2 = B

    FOR I = L TO N
      P(I) = X(I) + R1 * DeltaX(I)
    NEXT I

    F1 = Func(P())
    F2 = Fb
  END IF

  Iter = 0

  DO WHILE (Iter <= MaxIter) AND (ABS(R3 - R0) > Tol * (ABS(R1) + ABS(R2)))
    IF F2 < F1 THEN
      R0 = R1
      R1 = R2
      F1 = F2
      R2 = R1 + CGold * (R3 - R1)

      FOR I = L TO N
        P(I) = X(I) + R2 * DeltaX(I)
      NEXT I
    
      F2 = Func(P())
    ELSE
      R3 = R2
      R2 = R1
      F2 = F1
      R1 = R2 - CGold * (R2 - R0)
    
      FOR I = L TO N
        P(I) = X(I) + R1 * DeltaX(I)
      NEXT I
    
      F1 = Func(P())
    END IF

    Iter = Iter + 1
  LOOP

  IF F1 < F2 THEN
    R = R1
    Fmin = F1
  ELSE
    R = R2
    Fmin = F2
  END IF

  FOR I = L TO N
    X(I) = X(I) + R * DeltaX(I)
  NEXT I

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

