' ******************************************************************
' Function minimization by the simplex method
' ******************************************************************

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

#DEFINE OptOk 0  
' No error

#DEFINE OptNonConv -1  
' Non-convergence

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

#DEFINE MaxNum 1.797693134862315D+308   
' Max. floating point number: 2^1024

' ------------------------------------------------------------------
' Other constants
' ------------------------------------------------------------------

#DEFINE True  -1
#DEFINE False  0
  
' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

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

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

SUB UpdateSimplex (L AS INTEGER, N AS INTEGER, Hi AS INTEGER,_
                   Ystar AS DOUBLE, Pstar() AS DOUBLE,_
                   Pt() AS DOUBLE, F() AS DOUBLE)
' ------------------------------------------------------------------
' Update the current simplex
' ------------------------------------------------------------------
   
  DIM AS INTEGER I
   
  FOR I = L TO N
    Pt(I, Hi) = Pstar(I)
  NEXT I
 
  F(Hi) = Ystar
END SUB

SUB Simplex(      Func    AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                  X()     AS DOUBLE,_ 
                  MaxIter AS INTEGER,_ 
                  Tol     AS DOUBLE,_ 
            BYREF Fmin    AS DOUBLE)
' ------------------------------------------------------------------
' Minimizes function Func by the simplex method
' ------------------------------------------------------------------
' Input parameters  : Func    = function to be minimized
'                     X       = initial minimum coordinates
'                     MaxIter = maximum number of iterations
'                     Tol     = required precision
' ------------------------------------------------------------------
' Output parameters : X    = refined minimum coordinates
'                     Fmin = function value at minimum
' ------------------------------------------------------------------
' The function MathErr() returns one of the following codes:
'   OptOk      = no error
'   OptNonConv = non-convergence
' ------------------------------------------------------------------

  ' Quit if no iteration required
  IF MaxIter < 1 THEN
    Fmin = Func(X())
    ErrCode = OptOk
    EXIT SUB
  END IF

  DIM AS INTEGER I, J, L, N, M, Lo, Hi, Iter, Flag
  DIM AS DOUBLE  Corr, Cmax, F0, Ystar, Y2star, Z
        
  L = LBOUND(X)  ' Array bounds
  N = UBOUND(X)
  M = N + 1
                                            
  DIM AS DOUBLE X0(L TO N)                     ' Parameters
  DIM AS DOUBLE Pt(L TO N, L TO M)             ' Vertices
  DIM AS DOUBLE Fct(L TO M)                    ' Function values
  DIM AS DOUBLE Pbar(L TO N)                   ' Centroid coordinates
  DIM AS DOUBLE Pstar(L TO N), P2star(L TO N)  ' New vertices
      
  ' Construct initial simplex
  FOR J = L TO M
    FOR I = L TO N
      IF I = J - 1 THEN
        IF X0(I) <> 0 THEN X0(I) = 1.5 * X(I) ELSE X0(I) = 1
      ELSE
        X0(I) = X(I)
      END IF
      Pt(I, J) = X0(I)
    NEXT I
    Fct(J) = Func(X0())
  NEXT J
    
  Iter = 0
  F0 = MaxNum
  FOR I = L TO N
    X0(I) = X(I)    ' Initial parameters
  NEXT I
    
  DO
    Lo = 1          ' Index of vertex where f is minimal
    Hi = 1          ' Index of vertex where f is maximal
    FOR I = L TO M
      IF Fct(I) < Fct(Lo) THEN
        Lo = I
      ELSEIF Fct(I) > Fct(Hi) THEN
        Hi = I
      END IF
    NEXT I
                                
    IF Fct(Lo) < F0 THEN  ' New minimum found
      Iter = Iter + 1
      F0 = Fct(Lo)
    END IF
    
    ' Find centroid of points other than Pt(Hi)
    FOR I = L TO N
      Z = 0
      FOR J = L TO M
        IF J <> Hi THEN Z = Z + Pt(I, J)
      NEXT J
      Pbar(I) = Z / (N - L + 1)
    NEXT I
    
    ' Reflect worst vertex through centroid
    FOR I = L TO N
      Pstar(I) = 2 * Pbar(I) - Pt(I, Hi)
    NEXT I
    
    Ystar = Func(Pstar())
    
    ' If reflection successful, try extension
    IF Ystar < Fct(Lo) THEN
      FOR I = L TO N
        P2star(I) = 3 * Pstar(I) - 2 * Pbar(I)
      NEXT I
    
      Y2star = Func(P2star())
    
      ' Retain extension or contraction
      IF Y2star < Fct(Lo) THEN
        UpdateSimplex L, N, Hi, Y2star, P2star(), Pt(), Fct()
      ELSE
        UpdateSimplex L, N, Hi, Ystar, Pstar(), Pt(), Fct()
      END IF
    ELSE
      I = 0
      Flag = False
      DO
        I = I + 1
        IF I <> Hi AND Fct(I) > Ystar THEN Flag = True
      LOOP UNTIL Flag OR I = M
    
      IF Flag THEN
        UpdateSimplex L, N, Hi, Ystar, Pstar(), Pt(), Fct()
      ELSE
        ' Contraction on the reflection side of the centroid
        IF Ystar <= Fct(Hi) THEN
          UpdateSimplex L, N, Hi, Ystar, Pstar(), Pt(), Fct()
        END IF

        ' Contraction on the opposite side of the centroid
        FOR I = L TO N
          P2star(I) = .5 * (Pt(I, Hi) + Pbar(I))
        NEXT I
    
        Y2star = Func(P2star())
    
        IF Y2star <= Fct(Hi) THEN
          UpdateSimplex L, N, Hi, Y2star, P2star(), Pt(), Fct()
        ELSE
          ' Contract whole simplex
          FOR J = L TO M
            FOR I = L TO N
              Pt(I, J) = .5 * (Pt(I, J) + Pt(I, Lo))
            NEXT I
          NEXT J
        END IF
      END IF
    END IF
    
    FOR I = L TO N
      X0(I) = Pt(I, Lo)  ' New parameter
    NEXT I
       
    ' Test convergence
    Cmax = 0
    FOR I = L TO N
      Corr = ABS(Pt(I, Hi) - Pt(I, Lo))
      IF Corr > Cmax THEN Cmax = Corr
    NEXT I

    ' Update iteration count
    Iter = Iter + 1
    IF Iter > MaxIter THEN
      ErrCode = OptNonConv
      EXIT SUB
    END IF
  LOOP UNTIL Cmax < Tol
    
  Fmin = Fct(Lo)
  FOR I = L TO N
    X(I) = X0(I)
  NEXT I
    
  ErrCode = OptOk
END SUB



