C ======================================================================
      SUBROUTINE CREATE(M,S,N,XL,XR,FVAL,WA)
C ======================================================================
*  Creates a population of M members.
*  Each member is an array with N components (i.e. a point in R^N)
*  XL(I), XR(I) are the lower and upper bounds for the I-th component. 
*  S(N,M) contains the population, i.e M,  N-vectors.
*  FVAL(J) contains the value of the objective function evaluated at the 
*  point corresponding to J-th column of S.
*  WA is a work array. 
** INPUT:  M,N,XL,XR
** OUTPUT: S,FVAL
      implicit double precision (a-h,o-z)
      dimension s(n,m),xl(n),xr(n),wa(n),fval(m)
      do 1 i=1,m
         do 2 j=1,n
           s(j,i) = (xr(j)-xl(j))*ranm()+xl(j)
           wa(j) = s(j,i)
2        continue
           fval(i) = funmin(wa,n)
1     continue
      end
C ======================================================================
      subroutine crossd(p1,p2,n,c1,c2) 
C ======================================================================
*  Implements the discrete crossover operation, between two parents.
*  P1, P2 are the parent members.
*  C1, C2 are the children  members.
*      
      implicit double precision (a-h,o-z)
      dimension p1(n),p2(n),c1(n),c2(n)
      logical l1
      do 1 i=1,n
         xi = ranm()
         l1 = xi .le. 5.d-1 
         if (l1) then
                 c1(i) = p1(i)
                 c2(i) = p2(i)
         else
                 c1(i) = p2(i)
                 c2(i) = p1(i)
         endif
1     continue
      end         
C ======================================================================
      subroutine crossi(p1,p2,n,c1,c2,xl,xr) 
C ======================================================================
*  Implements the intermediate crossover operation, between two parents.
*  P1, P2 are the parent members.
*  C1, C2 are the children  members.
*  Care is taken that the children are inside the box.
      implicit double precision (a-h,o-z)
      parameter (d=0.25d0)
      dimension p1(*),p2(*),c1(*),c2(*),xl(*),xr(*)
      da = 1+2*d
      a0 = -d
      do 1 i=1,n
         alfa1 = da*ranm()+a0
         alfa2 = da*ranm()+a0
         c1(i) = p1(i) + alfa1*(p2(i)-p1(i))
           if (c1(i) .lt. xl(i) ) c1(i) = xl(i)
           if (c1(i) .gt. xr(i) ) c1(i) = xr(i)
         c2(i) = p1(i) + alfa2*(p2(i)-p1(i))
           if (c2(i) .lt. xl(i) ) c2(i) = xl(i)
           if (c2(i) .gt. xr(i) ) c2(i) = xr(i)
1     continue
      end
C ======================================================================
      subroutine calpro(fval,m,pr) 
C ======================================================================
*  Calculates the selection probabilities.
** INPUT:   FVAL, M
** OUTPUT:  PR
      implicit double precision (a-h,o-z)
      parameter (eps = 1.d-7)
      dimension fval(m),pr(m)
      sum = 0
      fmax = fval(1)
         do 1 i = 1,m
            sum = sum + fval(i)
            if (fval(i) .gt. fmax) fmax = fval(i)
1        continue
         div = m*fmax - sum
         if (div .eq. 0.d0) STOP 'CALPRO'
         do 2 i =1,m
            pr(i) = (fmax - fval(i))/div
            pr(i) = (pr(i)+eps)/(1+m*eps)
2        continue
      end
C ======================================================================
      subroutine select(pr,m,j) 
C ======================================================================
*  Implements the selection operation.
*  Selects according to the probabilities ( PR(M) ) one point.
**  INPUT: PR, M
**  OUTPUT: J, the index of the selected point.

      implicit double precision (a-h,o-z)
      dimension pr(m)
      xi = ranm()
      seg = 0.d0
          do 1 i = 1,m
             seg = seg + pr(i)
             if (xi .le. seg) Then
                     j=i
                     return
             endif
1         continue
      end
C ======================================================================
      subroutine equate(POP,N,M,S)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension POP(n,m), S(n,m)
      do 1 i = 1,n
         do 2 j = 1,m
             POP(i,j) = S(i,j)
2        continue
1     continue
      end
C ======================================================================
      subroutine evaluate(POP,n,m,FVAL,wa)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension POP(n,m),FVAL(m),wa(n)
      do 1 i=1,m
         do 2 j=1,n
           wa(j) = POP(j,i)
2        continue
           fval(i) = funmin(wa,n)
1     continue
      end
C ======================================================================
      subroutine PICKG(POP,N,M,ICL,P)
C ======================================================================
*  Selects the ICL column of POP, and stores it to P
      implicit double precision (a-h,o-z)
      dimension POP(n,m),P(n) 
      do 1 i=1,n
         P(i) = POP(i,ICL)
1     continue
      end
C ======================================================================
      subroutine SETCOL(POP,N,M,ICL,P)
C ======================================================================
*  Assigns P, to the ICL column of POP.
      implicit double precision (a-h,o-z)
      dimension POP(n,m),P(n) 
      do 1 i=1,n
         POP(i,ICL)=P(i)
1     continue
      end
C ======================================================================
      subroutine MUTATE(XL,XR,PA1,n,IGENER,M,PMUTE,CH1)
C ======================================================================
      implicit double precision (a-h,o-z)
      dimension PA1(n),CH1(n),XL(n),XR(n)
      b = (1.d0-dble(IGENER)/dble(M))*PMUTE
      do 1 i=1,n
         xi = ranm()
         yi = ranm()
         if (xi .gt. 0.5d0) Then
            ch1(i) = pa1(i) + (xr(i)-pa1(i))*(1.d0-yi**b)
         else
            ch1(i) = pa1(i) + (pa1(i)-xl(i))*(1.d0-yi**b)
         endif
1     continue
      end
C ======================================================================
      subroutine ELITE(FVAL,M,INDX,K,SCR)
C ======================================================================
*  INPUT : M, FVAL(M),K
*  OUTPUT: INDX(M), the positions of the K smallest values in FVAL(M)
*  SCR(M) is a workspace array
      implicit double precision (a-h,o-z)
      dimension FVAL(M),INDX(M),SCR(M)
      fmax = fval(1)
      do 10 i=1,m
        scr(i) = fval(i)
        indx(i) = 0
        if (fmax .lt. fval(i))fmax = fval(i)
10    continue
      do 1 i=1,k
        call minel(scr,m,iel)
          indx(i) = iel
          scr(iel)= fmax
1     continue
      end
C ======================================================================
      subroutine MINEL(FVAL,M,IEL)
C ======================================================================
*  INPUT : M, FVAL(M)
*  OUTPUT: IEL, the position of the smallest value in FVAL(M)
      implicit double precision (a-h,o-z)
      dimension FVAL(M)
      fmin = fval(1)
      iel = 1
        do 1 i = 1,m
          if (fval(i) .le. fmin) Then
            fmin = fval(i)
            iel = i
          endif
1       continue
      end            
C ======================================================================
      PROGRAM GENOMINO
C ======================================================================
      implicit double precision (a-h,o-z)
      parameter (n= 7, m= 200, ng =  800)
      parameter (delpc =   5.d-2, selpc=  5.d-2,mutpc=  5.d-2)
      parameter (PMUTE=  5.d0)
      dimension POP(n,m),S(n,m),X(n),FVAL(m),PROB(m),WA(n)      
      dimension XL(n),XR(n),PA1(n),PA2(n),CH1(n),CH2(n)
      dimension INDX(m),SCR(m),IXAT(n),icode(4)
      data XL/(n)*-5.d0/
      data XR/(n)*5.d0/
      data IXAT/(n)*1/
      data icode/1,1,1,0/
      ndel = max(1.d0,delpc*m)
      nsel = max(1.d0,selpc*m)
      nmut = max(1.d0,mutpc*m)
      ncro = m-ndel-nsel-nmut
C  NCRO must be even. If not decrement it by one and increment NMUT
      if ((ncro/2)*2 .ne. ncro) Then
              ncro = ncro -1
              nmut = nmut +1
      endif        
      write(*,'(57(''=''))')
      write(*,*)'GENOMINO RUNS WITH THE FOLLOWING DATA:'
      write(*,*) 'D-Elitism number: ',ndel
      write(*,*) 'P-Elitism number: ',nsel
      write(*,*) 'Crossover number: ',ncro
      write(*,*) 'Mutation  number: ',nmut
      write(*,*) 'Mutation    Rate: ',pmute
      write(*,*) 'Population  Size: ',m
      write(*,*) 'Max. Generations: ',ng 
      write(*,'(57(''=''))')
C  Create initial generation
      CALL CREATE(M,POP,N,XL,XR,FVAL,WA)
      IGENER = 1
30    CONTINUE
      if (IGENER.eq.ng) go to 31 
      IGENER = IGENER + 1
*  Calculate probabilities from  FVAL
      CALL calpro(fval,m,PROB) 
*  Find the positions of the NDEL lowest values in FVAL, and store them in
*   array INDX(M). 
      CALL ELITE(FVAL,M,INDX,NDEL,SCR)
      do 21 i =1,ndel
         iel = indx(i)
         CALL PICKG(POP,N,M,IEL,WA)
*  Keep the best NDEL points in the next generation.
*  (Deterministic  Elitism)
*  S(n,m) is the next generation population.
         CALL SETCOL(S,N,M,I,WA)
21    continue         
      IFIL = ndel
*  Selection Operation. (Roulette Wheel)
*  Select for the next generation NSEL points, probabilistically.
*  Points corresponding to lower values of the objective, have higher 
*  probability to "survive". (Probabilistic Elitism).
      do 1 is = 1,nsel
         IFIL = IFIL + 1
         CALL select(prob,m,j) 
         CALL PICKG(POP,N,M,J,WA)
*  Start forming the next generation population S.
         CALL SETCOL(S,N,M,IFIL,WA)
1     continue
      do 2 ic = 1,ncro/2
         CALL select(PROB,M,J1)
*  To avoid identical parents, temporarily set PROB(J1)=0.
         temp = PROB(j1)
         PROB(J1)=0.
         CALL select(PROB,M,J2)
*  Reset PROB(J1) back to its true value.
         PROB(J1) = temp
*  Construct the parents PA1, PA2
         CALL PICKG(POP,N,M,J1,PA1)
         CALL PICKG(POP,N,M,J2,PA2)
*  Crossover Operation: Create the children CH1, CH2
*             CALL crossd(pa1,pa2,n,ch1,ch2) 
              CALL crossi(pa1,pa2,n,ch1,ch2,xl,xr) 
*  Continue filling the next generation population S.
         IFIL = IFIL + 1
         CALL SETCOL(S,N,M,IFIL,CH1)
         IFIL = IFIL + 1
         CALL SETCOL(S,N,M,IFIL,CH2)
2     CONTINUE
*  Mutation Operation
         do 3 j=1,nmut
           IFIL = IFIL + 1
           CALL select(PROB,M,J1)
           CALL PICKG(POP,N,M,J1,PA1)
           CALL MUTATE(XL,XR,PA1,n,IGENER,M,PMUTE,CH1)
           CALL SETCOL(S,N,M,IFIL,CH1)
3        continue
c  Start a new generation
      CALL equate(POP,N,M,S)
      CALL evaluate(POP,n,m,FVAL,wa)
            CALL MINEL(FVAL,M,IEL)
*           write(*,'(1x,5(g14.7,1x))') FVAL(iel),(S(j,iel),j=1,n) 
      go to 30
31    continue
      write(*,88)'GA:','CALLS','MINIMUM','POINT'
88    format(1x,a8,T13,a5,T19,a13,T32,A17)
      write(*,89) m*ng,FVAL(IEL),(S(j,iel),j=1,n)
            write(*,'(57(''-''))')
89    format(T10,i7,T20,D14.7,T40,d14.7,10(/t40,d14.7))
            NOF = 0
            CALL PICKG(POP,N,M,IEL,WA)
            CALL OPTIMA(N,NOF,wa,VAL,XL,XR,IXAT,ICODE,
     >          'in.dat','/dev/null',GRMSO,NF,NOG,NHO,NJO)
            write(*,88)'MERLIN:'
            write(*,89) nf,VAL,(wa(j),j=1,n)
            write(*,'(57(''=''))')
            write(*,*)'TOTAL FUNCTION CALLS: ',m*ng+nf
            write(*,*)'TOTAL GRADIENT CALLS: ',nog
      end
