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