' ******************************************************************
' Optimization by Simulated Annealing
'
' Adapted from Fortran program SIMANN.FOR by Bill Goffe
' (http://www.netlib.org/simann)
' ******************************************************************

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

#DEFINE MatOk 0
' No error

#DEFINE MatErrDim -3
' Non-compatible dimensions

' ------------------------------------------------------------------
' Mathematical constant
' ------------------------------------------------------------------

#DEFINE InvLn2 1.442695040888963#
' 1/Ln(2)

' ------------------------------------------------------------------
' Machine-dependent constant
' ------------------------------------------------------------------

#DEFINE MaxLog 709.782712893384#
' Max. argument for Exp

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

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

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

DECLARE FUNCTION Median(X()    AS DOUBLE, _
                        Sorted AS INTEGER = -1) AS DOUBLE
' Median of sample X

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

' Strings used for printing results in log file
CONST Hdr  = " Iter          T                F           Inc     Acc"
CONST Mask = "#####    ##.#####^^^^^    ##.#####^^^^^    ####    ####"

' Global variables
DIM SHARED AS INTEGER SA_NT        = 5    ' Loops at constant temp.
DIM SHARED AS INTEGER SA_NS        = 15   ' Loops before step adjust
DIM SHARED AS DOUBLE  SA_RT        = 0.9  ' Temp. reduction factor
DIM SHARED AS INTEGER SA_NCycles   = 1    ' Number of cycles
DIM SHARED AS INTEGER WriteLogFile = 0    ' Log file ?
DIM SHARED AS INTEGER FF                  ' File handle

SUB InitSAParams(NT      AS INTEGER = 5,   _
                 NS      AS INTEGER = 15,  _
                 RT      AS DOUBLE  = 0.9, _
                 NCycles AS INTEGER = 1)
' ------------------------------------------------------------------
' Initialize Simulated Annealing parameters
' ------------------------------------------------------------------

  IF NT > 0 THEN SA_NT = NT
  IF NS > 0 THEN SA_NS = NS
  IF RT > 0 AND RT < 1 THEN SA_RT = RT
  IF NCycles > 0 THEN SA_NCycles = NCycles
END SUB

SUB SA_CreateLogFile(LogFileName AS STRING = "simann.txt")
' ------------------------------------------------------------------
' Initialize log file
' ------------------------------------------------------------------

  FF = FREEFILE
  OPEN LogFileName FOR OUTPUT AS #FF
  WriteLogFile = -1
END SUB

FUNCTION InitTemp(Func    AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                  X()     AS DOUBLE, _
                  Xmin()  AS DOUBLE, _
                  Range() AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Computes the initial temperature so that the probability
' of accepting an increase of the function is about 0.5
' ------------------------------------------------------------------

  CONST N_EVAL = 50      ' Number of function evaluations

  DIM AS DOUBLE  F, F1   ' Function values
  DIM AS INTEGER N_inc   ' Number of function increases
  DIM AS INTEGER Lb, Ub  ' Array bounds
  DIM AS INTEGER I       ' Index of function evaluation
  DIM AS INTEGER K       ' Index of parameter

  REDIM AS DOUBLE DeltaF(1 TO N_EVAL)  ' Function increases

  Lb = LBOUND(X)
  Ub = UBOUND(X)

  N_inc = 0
  F = Func(X())

  ' Compute N_EVAL function values, changing each parameter in turn
  K = Lb
  FOR I = 1 TO N_EVAL
    X(K) = Xmin(K) + RND * Range(K)
    F1 = Func(X())
    IF F1 > F THEN
      N_inc = N_inc + 1
      DeltaF(N_inc) = F1 - F
    END IF
    F = F1
    K = K + 1
    IF K > Ub THEN K = Lb
  NEXT I

  ' The median M of these N_inc increases has a probability of 1/2.
  ' From Boltzmann's formula: Exp(-M/T) = 1/2 ==> T = M / Ln(2)
  IF N_inc > 0 THEN
    REDIM PRESERVE DeltaF(1 TO N_inc)
    InitTemp = Median(DeltaF(), 0) * InvLn2
  ELSE
    InitTemp = 1.0
  END IF
END FUNCTION

FUNCTION Accept(      DeltaF AS DOUBLE, _
                      T      AS DOUBLE, _
                BYREF N_inc  AS INTEGER, _
                BYREF N_acc  AS INTEGER) AS INTEGER
' ------------------------------------------------------------------
' Checks if a variation DeltaF of the function at temperature T is
' acceptable. Updates the counters N_inc (number of increases of the
' function) and N_acc (number of accepted increases).
' ------------------------------------------------------------------

  DIM AS DOUBLE X

  IF DeltaF <= 0.0 THEN
    Accept = -1
    EXIT FUNCTION
  END IF

  N_inc = N_inc + 1
  X = DeltaF / T

  IF X > MaxLog THEN  ' Exp(-X) ~ 0
    Accept = 0
  ELSEIF EXP(- X) > RND THEN
    Accept = -1
    N_acc = N_acc + 1
  ELSE
    Accept = 0
  END IF
END FUNCTION

SUB SimAnnCycle(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                      X()    AS DOUBLE, _
                      Xmin() AS DOUBLE, _
                      Xmax() AS DOUBLE, _
                BYREF Fmin   AS DOUBLE)
' ----------------------------------------------------------------------
' Performs one cycle of simulated annealing
' ----------------------------------------------------------------------

  CONST SFact   = 2.0       ' Factor for step reduction
  CONST MinTemp = 1.0E-300  ' Minimal temperature
  CONST MinFunc = 1.0E-300  ' Minimal function value

  DIM AS INTEGER I, Iter, J, K, N_inc, N_acc, Lb, Ub
  DIM AS DOUBLE  F, F1, DeltaF, Ratio, T, OldX

  Lb = LBOUND(X)
  Ub = UBOUND(X)

  DIM AS DOUBLE  DeltaX(Lb TO Ub)
  DIM AS DOUBLE  Xopt(Lb TO Ub)
  DIM AS DOUBLE  Range(Lb TO Ub)
  DIM AS INTEGER Nacc(Lb TO Ub)

  ' Determine parameter range, step and optimum
  FOR K = Lb TO Ub
    Range(K) = Xmax(K) - Xmin(K)
    DeltaX(K) = 0.5 * Range(K)
    Xopt(K) = X(K)
  NEXT K

  ' Initialize function values
  F = Func(X())
  Fmin = F

  ' Initialize temperature and iteration count
  T = InitTemp(Func, X(), Xmin(), Range())
  Iter = 0

  DO
    N_inc = 0
    N_acc = 0

    ' Perform SA_NT evaluations at constant temperature
    FOR I = 1 TO SA_NT
      FOR J = 1 TO SA_NS
        FOR K = Lb TO Ub
          ' Save current parameter value
          OldX = X(K)

          ' Pick new value, keeping it within Range
          X(K) = X(K) + (2.0 * RND - 1.0) * DeltaX(K)
          IF X(K) < Xmin(K) OR X(K) > Xmax(K) THEN
            X(K) = Xmin(K) + RND * Range(K)
          END IF

          ' Compute new function value
          F1 = Func(X())
          DeltaF = F1 - F

          ' Check for acceptance
          IF Accept(DeltaF, T, N_inc, N_acc) THEN
            Nacc(K) = Nacc(K) + 1
            F = F1
          ELSE
            ' Restore parameter value
            X(K) = OldX
          END IF

          ' Update minimum if necessary
          IF F < Fmin THEN
            Xopt(K) = X(K)
            Fmin = F
          END IF
        NEXT K
      NEXT J

      ' Ajust step length to maintain an acceptance
      ' ratio of about 50% for each parameter
      FOR K = Lb TO Ub
        Ratio = Nacc(K) / SA_NS
        IF Ratio > 0.6 THEN
          ' Increase step length, keeping it within Range
          DeltaX(K) = DeltaX(K) * (1.0 + ((Ratio - 0.6) / 0.4) * SFact)
          IF DeltaX(K) > Range(K) THEN DeltaX(K) = Range(K)
        ELSEIF Ratio < 0.4 THEN
          ' Reduce step length
          DeltaX(K) = DeltaX(K) / (1.0 + ((0.4 - Ratio) / 0.4) * SFact)
        END IF

        ' Restore counter
        Nacc(K) = 0
      NEXT K
    NEXT I

    IF WriteLogFile THEN
      PRINT #FF, USING Mask; Iter; T; Fmin; N_inc; N_acc
    END IF

    ' Update temperature and iteration count
    T = T * SA_RT
    Iter = Iter + 1
  LOOP UNTIL N_acc = 0 OR T < MinTemp OR ABS(Fmin) < MinFunc

  FOR K = Lb TO Ub
    X(K) = Xopt(K)
  NEXT K
END SUB

SUB SimAnn(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                 X()    AS DOUBLE, _
                 Xmin() AS DOUBLE, _
                 Xmax() AS DOUBLE, _
           BYREF Fmin   AS DOUBLE)
' ------------------------------------------------------------------
' Minimization of a function of several variables
' by simulated annealing
' ------------------------------------------------------------------
' Input parameters : Func = objective function to be minimized
'                    X    = initial minimum coordinates
'                    Xmin = minimum value of X
'                    Xmax = maximum value of X
' ------------------------------------------------------------------
' Output parameters: X    = refined minimum coordinates
'                    Fmin = function value at minimum
' ------------------------------------------------------------------

  DIM AS INTEGER Cycle, Lb, Ub

  Lb = LBOUND(X)
  Ub = UBOUND(X)

  IF LBOUND(Xmin) <> Lb OR UBOUND(Xmin) <> Ub OR _
     LBOUND(Xmax) <> Lb OR UBOUND(Xmax) <> Ub THEN
       ErrCode = MatErrDim
       EXIT SUB
  END IF

  ErrCode = MatOk

  ' Initialize the random number generator
  RANDOMIZE TIMER

  FOR Cycle = 1 TO SA_NCycles
    IF WriteLogFile THEN
      PRINT #FF, ""
      PRINT #FF, "Simulated annealing: Cycle "; Cycle
      PRINT #FF, Hdr
    END IF
    SimAnnCycle(Func, X(), Xmin(), Xmax(), Fmin)
  NEXT Cycle

  IF WriteLogFile THEN
    CLOSE #FF
    WriteLogFile = 0
  END IF
END SUB
