' ******************************************************************
' Optimization by Genetic Algorithm
' ******************************************************************
' Ref.:  E. Perrin, A. Mandrille, M. Oumoun, C. Fonteix & I. Marc
'        Optimisation globale par strategie d'evolution
'        Technique utilisant la genetique des individus diploides
'        Recherche operationnelle / Operations Research
'        1997, 31, 161-201
' Thanks to Magali Camut for her contribution
' ******************************************************************

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

#DEFINE MatOk 0
' No error

#DEFINE MatErrDim -3
' Non-compatible dimensions

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

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

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

DECLARE FUNCTION Max (X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Maximum of 2 reals

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

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

' Global variables
DIM SHARED AS INTEGER GA_NP = 200  ' Population size
DIM SHARED AS INTEGER GA_NG = 40   ' Max number of generations
DIM SHARED AS DOUBLE  GA_SR = 0.6  ' Survival rate
DIM SHARED AS DOUBLE  GA_MR = 0.1  ' Mutation rate
DIM SHARED AS DOUBLE  GA_HR = 0.5  ' Proportion of homozygotes

DIM SHARED AS INTEGER WriteLogFile = 0  ' Log file ?
DIM SHARED AS INTEGER FF                ' File handle

SUB InitGAParams(NP AS INTEGER = 200, _
                 NG AS INTEGER = 40,  _
                 SR AS DOUBLE  = 0.5, _
                 MR AS DOUBLE  = 0.1, _
                 HR AS DOUBLE  = 0.5)
' ------------------------------------------------------------------
' Initialize Genetic Algorithm parameters
' ------------------------------------------------------------------

  IF NP > 0 THEN GA_NP = NP
  IF NG > 0 THEN GA_NG = NG

  IF SR > 0 AND SR < 1 THEN GA_SR = SR
  IF MR > 0 AND MR < 1 THEN GA_MR = MR
  IF HR > 0 AND HR < 1 THEN GA_HR = HR
END SUB

SUB GA_CreateLogFile(LogFileName AS STRING = "genalg.txt")
' ------------------------------------------------------------------
' Initialize log file
' ------------------------------------------------------------------

  FF = FREEFILE
  OPEN LogFileName FOR OUTPUT AS #FF
  PRINT #FF, "Genetic Algorithm"
  PRINT #FF, Hdr
  WriteLogFile = -1
END SUB

SUB Mutate(I AS INTEGER, _
           C1() AS DOUBLE, C2() AS DOUBLE, D() AS DOUBLE, P() AS DOUBLE, _
           Xmin() AS DOUBLE, Range() AS DOUBLE)
' ------------------------------------------------------------------
' Mutate individual I
' ------------------------------------------------------------------

  DIM AS INTEGER J

  FOR J = LBOUND(Xmin) TO UBOUND(Xmin)
    C1(I,J) = Xmin(J) + RND * Range(J)
    C2(I,J) = Xmin(J) + RND * Range(J)
    D(I,J) = RND
    P(I,J) = D(I,J) * C1(I,J) + (1 - D(I,J)) * C2(I,J)
  NEXT J
END SUB

SUB Cross(I1 AS INTEGER, I2 AS INTEGER, I AS INTEGER, _
          C1() AS DOUBLE, C2() AS DOUBLE, D() AS DOUBLE, P() AS DOUBLE)
' ------------------------------------------------------------------
' Cross two individuals I1 and I2 --> new individual I
' ------------------------------------------------------------------

  DIM AS INTEGER J

  FOR J = LBOUND(P,2) TO UBOUND(P,2)
    IF RND < 0.5 THEN C1(I,J) = C1(I1,J) ELSE C1(I,J) = C1(I2,J)
    IF RND < 0.5 THEN C2(I,J) = C2(I1,J) ELSE C2(I,J) = C2(I2,J)
    D(I,J) = RND
    P(I,J) = D(I,J) * C1(I,J) + (1 - D(I,J)) * C2(I,J)
  NEXT J
END SUB

SUB Homozygote(I AS INTEGER, _
               C1() AS DOUBLE, C2() AS DOUBLE, P() AS DOUBLE)
' ------------------------------------------------------------------
' Make individual I homozygous
' ------------------------------------------------------------------

  DIM AS INTEGER J

  FOR J = LBOUND(P,2) TO UBOUND(P,2)
    C1(I,J) = P(I,J)
    C2(I,J) = P(I,J)
  NEXT J
END SUB

FUNCTION GA_Func(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                 I AS INTEGER, P() AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Computes objective function for individual I
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, J

  Lb = LBOUND(P, 2)
  Ub = UBOUND(P, 2)

  DIM AS DOUBLE X(Lb TO Ub)

  FOR J = Lb TO Ub
    X(J) = P(I,J)
  NEXT J

  GA_Func = Func(X())
END FUNCTION

SUB CompFunc(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE, X() AS DOUBLE, _
             C1() AS DOUBLE, C2() AS DOUBLE, D() AS DOUBLE, P() AS DOUBLE, _
             F() AS DOUBLE, BYREF Fmin AS DOUBLE)
' ------------------------------------------------------------------
' Computes function values
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J, K
  DIM AS DOUBLE  A

  STATIC Iter AS INTEGER = 0

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

  ' Compute function values
  FOR I = 1 TO GA_NP
    F(I) = GA_Func(Func, I, P())
  NEXT I

  ' Sort population according to function values
  FOR I = 1 TO GA_NP - 1
    K = I
    A = F(I)

    FOR J = I + 1 TO GA_NP
      IF F(J) < A THEN
        K = J
        A = F(J)
      END IF
    NEXT J

    SWAP F(I), F(K)

    FOR J = Lb TO Ub
      SWAP C1(I,J), C1(K,J)
      SWAP C2(I,J), C2(K,J)
      SWAP D(I,J), D(K,J)
      SWAP P(I,J), P(K,J)
    NEXT J
  NEXT I

  ' Update log file if necessary
  IF WriteLogFile THEN
    PRINT #FF, USING Mask; Iter; F(1)
  END IF

  ' Update minimum
  IF F(1) < Fmin THEN
    Fmin = F(1)
    FOR J = Lb TO Ub
      X(J) = P(1,J)
    NEXT J
  END IF

  Iter = Iter + 1
END SUB

SUB GenPop(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE, NS AS INTEGER, _
           C1() AS DOUBLE, C2() AS DOUBLE, D() AS DOUBLE, P() AS DOUBLE, _
           F() AS DOUBLE, Xmin() AS DOUBLE, Range() AS DOUBLE)
' ------------------------------------------------------------------
' Generates new population
' ------------------------------------------------------------------

  DIM AS INTEGER I, I1, I2
  DIM AS DOUBLE  F0

  FOR I = NS + 1 TO GA_NP
    I1 = INT(RND * NS) + 1

    DO
      I2 = INT(RND * NS) + 1
    LOOP UNTIL I2 <> I1

    F0 = Max(F(I1), F(I2))

    DO
      Cross I1, I2, I, C1(), C2(), D(), P()
    LOOP UNTIL GA_Func(Func, I, P()) <= F0
  NEXT I

  FOR I = 1 TO GA_NP
    IF RND < GA_MR THEN Mutate I, C1(), C2(), D(), P(), Xmin(), Range()
    IF RND < GA_HR THEN Homozygote I, C1(), C2(), P()
  NEXT I
END SUB

SUB GenAlg(      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 genetic algorithm
' ------------------------------------------------------------------
' 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 Lb, Ub, I, J, NS

  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

  ' Dimension arrays
  DIM AS DOUBLE C1(1 TO GA_NP, Lb TO Ub)
  DIM AS DOUBLE C2(1 TO GA_NP, Lb TO Ub)
  DIM AS DOUBLE D(1 TO GA_NP, Lb TO Ub)
  DIM AS DOUBLE P(1 TO GA_NP, Lb TO Ub)

  DIM AS DOUBLE F(1 TO GA_NP)

  DIM AS DOUBLE Range(Lb TO Ub)

  FOR I = Lb TO Ub
    Range(I) = Xmax(I) - Xmin(I)
  NEXT I

  NS = GA_NP * GA_SR  ' Number of survivors

  ' Initialize population and function value
  Fmin = 1.0E+300
  FOR I = 1 TO GA_NP
    Mutate I, C1(), C2(), D(), P(), Xmin(), Range()
  NEXT I

  CompFunc Func, X(), C1(), C2(), D(), P(), F(), Fmin

  FOR I = 1 TO GA_NG
    GenPop Func, NS, C1(), C2(), D(), P(), F(), Xmin(), Range()
    CompFunc Func, X(), C1(), C2(), D(), P(), F(), Fmin
  NEXT I

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