1178 lines
37 KiB
FortranFixed
1178 lines
37 KiB
FortranFixed
c*********************************************************************
|
|
double precision function urand()
|
|
c=====================================================================
|
|
c Return the next pseudo-random deviate from a sequence which is
|
|
c uniformly distributed in the interval [0,1]
|
|
c
|
|
c Uses the function ran0, the "minimal standard" random number
|
|
c generator of Park and Miller (Comm. ACM 31, 1192-1201, Oct 1988;
|
|
c Comm. ACM 36 No. 7, 105-110, July 1993).
|
|
c=====================================================================
|
|
implicit none
|
|
c
|
|
c Input - none
|
|
c
|
|
c Output
|
|
|
|
c Local
|
|
double precision ran2
|
|
external ran2
|
|
c
|
|
urand = ran2()
|
|
return
|
|
end
|
|
c**********************************************************************
|
|
subroutine rqsort(n,a,p)
|
|
c======================================================================
|
|
c Return integer array p which indexes array a in increasing order.
|
|
c Array a is not disturbed. The Quicksort algorithm is used.
|
|
c
|
|
c B. G. Knapp, 86/12/23
|
|
c
|
|
c Reference: N. Wirth, Algorithms and Data Structures,
|
|
c Prentice-Hall, 1986
|
|
c======================================================================
|
|
implicit none
|
|
|
|
c Input:
|
|
integer n
|
|
double precision a(n)
|
|
|
|
c Output:
|
|
integer p(n)
|
|
|
|
c Constants
|
|
integer LGN, Q
|
|
parameter (LGN=32, Q=11)
|
|
c (LGN = log base 2 of maximum n;
|
|
c Q = smallest subfile to use quicksort on)
|
|
|
|
c Local:
|
|
double precision x
|
|
integer stackl(LGN),stackr(LGN),s,t,l,m,r,i,j
|
|
|
|
c Initialize the stack
|
|
stackl(1)=1
|
|
stackr(1)=n
|
|
s=1
|
|
|
|
c Initialize the pointer array
|
|
do 1 i=1,n
|
|
p(i)=i
|
|
1 continue
|
|
|
|
2 if (s.gt.0) then
|
|
l=stackl(s)
|
|
r=stackr(s)
|
|
s=s-1
|
|
|
|
3 if ((r-l).lt.Q) then
|
|
|
|
c Use straight insertion
|
|
do 6 i=l+1,r
|
|
t = p(i)
|
|
x = a(t)
|
|
do 4 j=i-1,l,-1
|
|
if (a(p(j)).le.x) goto 5
|
|
p(j+1) = p(j)
|
|
4 continue
|
|
j=l-1
|
|
5 p(j+1) = t
|
|
6 continue
|
|
else
|
|
|
|
c Use quicksort, with pivot as median of a(l), a(m), a(r)
|
|
m=(l+r)/2
|
|
t=p(m)
|
|
if (a(t).lt.a(p(l))) then
|
|
p(m)=p(l)
|
|
p(l)=t
|
|
t=p(m)
|
|
endif
|
|
if (a(t).gt.a(p(r))) then
|
|
p(m)=p(r)
|
|
p(r)=t
|
|
t=p(m)
|
|
if (a(t).lt.a(p(l))) then
|
|
p(m)=p(l)
|
|
p(l)=t
|
|
t=p(m)
|
|
endif
|
|
endif
|
|
|
|
c Partition
|
|
x=a(t)
|
|
i=l+1
|
|
j=r-1
|
|
7 if (i.le.j) then
|
|
8 if (a(p(i)).lt.x) then
|
|
i=i+1
|
|
goto 8
|
|
endif
|
|
9 if (x.lt.a(p(j))) then
|
|
j=j-1
|
|
goto 9
|
|
endif
|
|
if (i.le.j) then
|
|
t=p(i)
|
|
p(i)=p(j)
|
|
p(j)=t
|
|
i=i+1
|
|
j=j-1
|
|
endif
|
|
goto 7
|
|
endif
|
|
|
|
c Stack the larger subfile
|
|
s=s+1
|
|
if ((j-l).gt.(r-i)) then
|
|
stackl(s)=l
|
|
stackr(s)=j
|
|
l=i
|
|
else
|
|
stackl(s)=i
|
|
stackr(s)=r
|
|
r=j
|
|
endif
|
|
goto 3
|
|
endif
|
|
goto 2
|
|
endif
|
|
return
|
|
end
|
|
c***********************************************************************
|
|
subroutine pikaia(ff,n,ctrl,x,f,icondi)
|
|
c=======================================================================
|
|
c Optimization (maximization) of user-supplied "fitness" function
|
|
c ff over n-dimensional parameter space x using a basic genetic
|
|
c algorithm method.
|
|
c
|
|
c Paul Charbonneau & Barry Knapp
|
|
c High Altitude Observatory
|
|
c National Center for Atmospheric Research
|
|
c Boulder CO 80307-3000
|
|
c <paulchar@hao.ucar.edu>
|
|
c <knapp@hao.ucar.edu>
|
|
c
|
|
c Version 1.2 [ 2002 April 3 ]
|
|
c
|
|
c Genetic algorithms are heuristic search techniques that
|
|
c incorporate in a computational setting, the biological notion
|
|
c of evolution by means of natural selection. This subroutine
|
|
c implements the three basic operations of selection, crossover,
|
|
c and mutation, operating on "genotypes" encoded as strings.
|
|
c
|
|
c Version 1.2 differs from version 1.0 (December 1995) in that
|
|
c it includes (1) two-point crossover, (2) creep mutation, and
|
|
c (3) dynamical adjustment of the mutation rate based on metric
|
|
c distance in parameter space.
|
|
c
|
|
c References:
|
|
c
|
|
c Charbonneau, Paul. "An introduction to gemetic algorithms for
|
|
c numerical optimization", NCAR Technical Note TN-450+IA
|
|
c (April 2002)
|
|
c
|
|
c Charbonneau, Paul. "Release Notes for PIKAIA 1.2",
|
|
c NCAR Technical Note TN-451+STR (April 2002)
|
|
c
|
|
c Charbonneau, Paul, and Knapp, Barry. "A User's Guide
|
|
c to PIKAIA 1.0" NCAR Technical Note TN-418+IA
|
|
c (December 1995)
|
|
c
|
|
c Goldberg, David E. Genetic Algorithms in Search, Optimization,
|
|
c & Machine Learning. Addison-Wesley, 1989.
|
|
c
|
|
c Davis, Lawrence, ed. Handbook of Genetic Algorithms.
|
|
c Van Nostrand Reinhold, 1991.
|
|
c
|
|
c=======================================================================
|
|
c USES: ff, urand, setctl, report, rnkpop, select, encode, decode,
|
|
c cross, mutate, genrep, stdrep, newpop, adjmut
|
|
implicit none
|
|
|
|
c Input:
|
|
integer n
|
|
double precision ff
|
|
external ff
|
|
c
|
|
c o Integer n is the parameter space dimension, i.e., the number
|
|
c of adjustable parameters.
|
|
c
|
|
c o Function ff is a user-supplied scalar function of n vari-
|
|
c ables, which must have the calling sequence f = ff(n,x), where
|
|
c x is a double precision parameter array of length n. This function must
|
|
c be written so as to bound all parameters to the interval [0,1];
|
|
c that is, the user must determine a priori bounds for the para-
|
|
c meter space, and ff must use these bounds to perform the appro-
|
|
c priate scalings to recover true parameter values in the
|
|
c a priori ranges.
|
|
c
|
|
c By convention, ff should return higher values for more optimal
|
|
c parameter values (i.e., individuals which are more "fit").
|
|
c For example, in fitting a function through data points, ff
|
|
c could return the inverse of chi**2.
|
|
c
|
|
c In most cases initialization code will have to be written
|
|
c (either in a driver or in a separate subroutine) which loads
|
|
c in data values and communicates with ff via one or more labeled
|
|
c common blocks. An example exercise driver and fitness function
|
|
c are provided in the accompanying file, xpkaia.f.
|
|
c
|
|
c
|
|
c Input/Output:
|
|
double precision ctrl(12)
|
|
c
|
|
c o Array ctrl is an array of control flags and parameters, to
|
|
c control the genetic behavior of the algorithm, and also printed
|
|
c output. A default value will be used for any control variable
|
|
c which is supplied with a value less than zero. On exit, ctrl
|
|
c contains the actual values used as control variables. The
|
|
c elements of ctrl and their defaults are:
|
|
c
|
|
c ctrl( 1) - number of individuals in a population (default
|
|
c is 100)
|
|
c ctrl( 2) - number of generations over which solution is
|
|
c to evolve (default is 500)
|
|
c ctrl( 3) - number of significant digits (i.e., number of
|
|
c genes) retained in chromosomal encoding (default
|
|
c is 6) (Note: This number is limited by the
|
|
c machine floating point precision. Most 32-bit
|
|
c floating point representations have only 6 full
|
|
c digits of precision. To achieve greater preci-
|
|
c sion this routine could be converted to double
|
|
c precision, but note that this would also require
|
|
c a double precision random number generator, which
|
|
c likely would not have more than 9 digits of
|
|
c precision if it used 4-byte integers internally.)
|
|
c ctrl( 4) - crossover probability; must be <= 1.0 (default
|
|
c is 0.85). If crossover takes place, either one
|
|
c or two splicing points are used, with equal
|
|
c probabilities
|
|
c ctrl( 5) - mutation mode; 1/2/3/4/5 (default is 2)
|
|
c 1=one-point mutation, fixed rate
|
|
c 2=one-point, adjustable rate based on fitness
|
|
c 3=one-point, adjustable rate based on distance
|
|
c 4=one-point+creep, fixed rate
|
|
c 5=one-point+creep, adjustable rate based on fitness
|
|
c 6=one-point+creep, adjustable rate based on distance
|
|
c ctrl( 6) - initial mutation rate; should be small (default
|
|
c is 0.005) (Note: the mutation rate is the proba-
|
|
c bility that any one gene locus will mutate in
|
|
c any one generation.)
|
|
c ctrl( 7) - minimum mutation rate; must be >= 0.0 (default
|
|
c is 0.0005)
|
|
c ctrl( 8) - maximum mutation rate; must be <= 1.0 (default
|
|
c is 0.25)
|
|
c ctrl( 9) - relative fitness differential; range from 0
|
|
c (none) to 1 (maximum). (default is 1.)
|
|
c ctrl(10) - reproduction plan; 1/2/3=Full generational
|
|
c replacement/Steady-state-replace-random/Steady-
|
|
c state-replace-worst (default is 3)
|
|
c ctrl(11) - elitism flag; 0/1=off/on (default is 0)
|
|
c (Applies only to reproduction plans 1 and 2)
|
|
c ctrl(12) - printed output 0/1/2=None/Minimal/Verbose
|
|
c (default is 0)
|
|
c
|
|
c
|
|
c Output:
|
|
double precision x(n), f
|
|
integer icondi
|
|
c
|
|
c o Array x(1:n) is the "fittest" (optimal) solution found,
|
|
c i.e., the solution which maximizes fitness function ff
|
|
c
|
|
c o Scalar f is the value of the fitness function at x
|
|
c
|
|
c o Integer icondi is an indicator of the success or failure
|
|
c of the call to pikaia (0=success; non-zero=failure)
|
|
c
|
|
c
|
|
c Constants
|
|
integer NMAX, PMAX, DMAX
|
|
parameter (NMAX = 200, PMAX = 5001, DMAX = 9)
|
|
c
|
|
c o NMAX is the maximum number of adjustable parameters
|
|
c (n <= NMAX). original NMAX was 32
|
|
c
|
|
c o PMAX is the maximum population (ctrl(1) <= PMAX)
|
|
c
|
|
c o DMAX is the maximum number of Genes (digits) per Chromosome
|
|
c segement (parameter) (ctrl(3) <= DMAX)
|
|
c
|
|
c
|
|
c Local variables
|
|
integer np, nd, ngen, imut, irep, ielite, ivrb, k, ip, ig,
|
|
+ ip1, ip2, new, newtot
|
|
double precision pcross, pmut, pmutmn, pmutmx, fdif
|
|
c
|
|
double precision ph(NMAX,2), oldph(NMAX,PMAX), newph(NMAX,PMAX)
|
|
c
|
|
integer gn1(NMAX*DMAX), gn2(NMAX*DMAX)
|
|
integer ifit(PMAX), jfit(PMAX)
|
|
double precision fitns(PMAX)
|
|
c
|
|
c User-supplied uniform random number generator
|
|
double precision urand
|
|
external urand
|
|
c
|
|
c Function urand should not take any arguments. If the user wishes
|
|
c to be able to initialize urand, so that the same sequence of
|
|
c random numbers can be repeated, this capability could be imple-
|
|
c mented with a separate subroutine, and called from the user's
|
|
c driver program. An example urand function (and initialization
|
|
c subroutine) which uses the function ran0 (the "minimal standard"
|
|
c random number generator of Park and Miller [Comm. ACM 31, 1192-
|
|
c 1201, Oct 1988; Comm. ACM 36 No. 7, 105-110, July 1993]) is
|
|
c provided.
|
|
c
|
|
c
|
|
c Set control variables from input and defaults
|
|
call setctl
|
|
+ (ctrl,n,np,ngen,nd,pcross,pmutmn,pmutmx,pmut,imut,
|
|
+ fdif,irep,ielite,ivrb,icondi)
|
|
|
|
c Make sure locally-dimensioned arrays are big enough
|
|
if (n.gt.NMAX .or. np.gt.PMAX .or. nd.gt.DMAX) then
|
|
write(*,*)n,NMAX,np,PMAX,nd,DMAX
|
|
|
|
write(*,*)
|
|
+ ' Number of parameters, population, or genes too large'
|
|
icondi = -1
|
|
return
|
|
endif
|
|
|
|
c Compute initial (random but bounded) phenotypes
|
|
do 1 ip=2,np
|
|
do 2 k=1,n
|
|
oldph(k,ip)=urand()
|
|
2 continue
|
|
1 continue
|
|
|
|
do k=1,n
|
|
oldph(k,1)=x(k)
|
|
enddo
|
|
do ip=1,np
|
|
fitns(ip) = ff(n,oldph(1,ip))
|
|
enddo
|
|
|
|
c Rank initial population by fitness order
|
|
call rnkpop(np,fitns,ifit,jfit)
|
|
|
|
c Main Generation Loop
|
|
do 10 ig=1,ngen
|
|
|
|
c Main Population Loop
|
|
newtot=0
|
|
do 20 ip=1,np/2
|
|
|
|
c 1. pick two parents
|
|
call select(np,jfit,fdif,ip1)
|
|
21 call select(np,jfit,fdif,ip2)
|
|
if (ip1.eq.ip2) goto 21
|
|
|
|
c 2. encode parent phenotypes
|
|
call encode(n,nd,oldph(1,ip1),gn1)
|
|
call encode(n,nd,oldph(1,ip2),gn2)
|
|
|
|
c 3. breed
|
|
call cross(n,nd,pcross,gn1,gn2)
|
|
call mutate(n,nd,pmut,gn1,imut)
|
|
call mutate(n,nd,pmut,gn2,imut)
|
|
|
|
c 4. decode offspring genotypes
|
|
call decode(n,nd,gn1,ph(1,1))
|
|
call decode(n,nd,gn2,ph(1,2))
|
|
|
|
c 5. insert into population
|
|
if (irep.eq.1) then
|
|
call genrep(NMAX,n,np,ip,ph,newph)
|
|
else
|
|
call stdrep(ff,NMAX,n,np,irep,ielite,
|
|
+ ph,oldph,fitns,ifit,jfit,new)
|
|
newtot = newtot+new
|
|
endif
|
|
|
|
c End of Main Population Loop
|
|
20 continue
|
|
|
|
c if running full generational replacement: swap populations
|
|
if (irep.eq.1)
|
|
+ call newpop(ff,ielite,NMAX,n,np,oldph,newph,
|
|
+ ifit,jfit,fitns,newtot)
|
|
|
|
c adjust mutation rate?
|
|
if (imut.eq.2 .or. imut.eq.3 .or. imut.eq.5 .or. imut.eq.6)
|
|
+ call adjmut(NMAX,n,np,oldph,fitns,ifit,pmutmn,pmutmx,
|
|
+ pmut,imut)
|
|
c
|
|
if (ivrb.gt.0) call report
|
|
+ (ivrb,NMAX,n,np,nd,oldph,fitns,ifit,pmut,ig,newtot)
|
|
|
|
c End of Main Generation Loop
|
|
10 continue
|
|
c
|
|
c Return best phenotype and its fitness
|
|
do 30 k=1,n
|
|
x(k) = oldph(k,ifit(np))
|
|
30 continue
|
|
f = fitns(ifit(np))
|
|
c
|
|
end
|
|
c********************************************************************
|
|
subroutine setctl
|
|
+ (ctrl,n,np,ngen,nd,pcross,pmutmn,pmutmx,pmut,imut,
|
|
+ fdif,irep,ielite,ivrb,icondi)
|
|
c===================================================================
|
|
c Set control variables and flags from input and defaults
|
|
c===================================================================
|
|
implicit none
|
|
c
|
|
c Input
|
|
integer n
|
|
c
|
|
c Input/Output
|
|
double precision ctrl(12)
|
|
c
|
|
c Output
|
|
integer np, ngen, nd, imut, irep, ielite, ivrb, icondi
|
|
double precision pcross, pmutmn, pmutmx, pmut, fdif
|
|
c
|
|
c Local
|
|
integer i
|
|
double precision DFAULT(12)
|
|
save DFAULT
|
|
data DFAULT /100.0d0,500.0d0,5.0d0,0.85d0,2.0d0,0.005d0,
|
|
&0.0005d0,0.25d0,1.0d0,1.0d0,1.0d0,0.0d0/
|
|
c
|
|
do 1 i=1,12
|
|
if (ctrl(i).lt.0.0d0) ctrl(i)=DFAULT(i)
|
|
1 continue
|
|
|
|
np = ctrl(1)
|
|
ngen = ctrl(2)
|
|
nd = ctrl(3)
|
|
pcross = ctrl(4)
|
|
imut = ctrl(5)
|
|
pmut = ctrl(6)
|
|
pmutmn = ctrl(7)
|
|
pmutmx = ctrl(8)
|
|
fdif = ctrl(9)
|
|
irep = ctrl(10)
|
|
ielite = ctrl(11)
|
|
ivrb = ctrl(12)
|
|
icondi = 0
|
|
c
|
|
c Print a header
|
|
if (ivrb.gt.0) then
|
|
|
|
write(*,2) ngen,np,n,nd,pcross,pmut,pmutmn,pmutmx,fdif
|
|
2 format(/1x,60('*'),/,
|
|
+ ' *',13x,'PIKAIA Genetic Algorithm Report ',13x,'*',/,
|
|
+ 1x,60('*'),//,
|
|
+ ' Number of Generations evolving: ',i4,/,
|
|
+ ' Individuals per generation: ',i4,/,
|
|
+ ' Number of Chromosome segments: ',i4,/,
|
|
+ ' Length of Chromosome segments: ',i4,/,
|
|
+ ' Crossover probability: ',f9.4,/,
|
|
+ ' Initial mutation rate: ',f9.4,/,
|
|
+ ' Minimum mutation rate: ',f9.4,/,
|
|
+ ' Maximum mutation rate: ',f9.4,/,
|
|
+ ' Relative fitness differential: ',f9.4)
|
|
if (imut.eq.1) write(*,3) 'Uniform, Constant Rate'
|
|
if (imut.eq.2) write(*,3) 'Uniform, Variable Rate (F)'
|
|
if (imut.eq.3) write(*,3) 'Uniform, Variable Rate (D)'
|
|
if (imut.eq.4) write(*,3) 'Uniform+Creep, Constant Rate'
|
|
if (imut.eq.5) write(*,3) 'Uniform+Creep, Variable Rate (F)'
|
|
if (imut.eq.6) write(*,3) 'Uniform+Creep, Variable Rate (D)'
|
|
3 format(
|
|
+ ' Mutation Mode: ',A)
|
|
if (irep.eq.1) write(*,4) 'Full generational replacement'
|
|
if (irep.eq.2) write(*,4) 'Steady-state-replace-random'
|
|
if (irep.eq.3) write(*,4) 'Steady-state-replace-worst'
|
|
4 format(
|
|
+ ' Reproduction Plan: ',A)
|
|
endif
|
|
|
|
c Check some control values
|
|
if (imut.ne.1 .and. imut.ne.2 .and. imut.ne.3 .and. imut.ne.4
|
|
+ .and. imut.ne.5 .and. imut.ne.6) then
|
|
write(*,10)
|
|
icondi = 5
|
|
endif
|
|
10 format(' ERROR: illegal value for imut (ctrl(5))')
|
|
|
|
if (fdif.gt.1.) then
|
|
write(*,11)
|
|
icondi = 9
|
|
endif
|
|
11 format(' ERROR: illegal value for fdif (ctrl(9))')
|
|
|
|
if (irep.ne.1 .and. irep.ne.2 .and. irep.ne.3) then
|
|
write(*,12)
|
|
icondi = 10
|
|
endif
|
|
12 format(' ERROR: illegal value for irep (ctrl(10))')
|
|
|
|
if (pcross.gt.1.0d0.or.pcross.lt.0.0d0)then
|
|
write(*,13)
|
|
icondi = 4
|
|
endif
|
|
13 format(' ERROR: illegal value for pcross (ctrl(4))')
|
|
|
|
if (ielite.ne.0 .and. ielite.ne.1) then
|
|
write(*,14)
|
|
icondi = 11
|
|
endif
|
|
14 format(' ERROR: illegal value for ielite (ctrl(11))')
|
|
|
|
if (irep.eq.1 .and. imut.eq.1 .and. pmut.gt.0.5 .and.
|
|
+ ielite.eq.0) then
|
|
write(*,15)
|
|
endif
|
|
15 format(' WARNING: dangerously high value for pmut (ctrl(6));',
|
|
+ /' (Should enforce elitism with ctrl(11)=1.)')
|
|
|
|
if (irep.eq.1 .and. imut.eq.2 .and. pmutmx.gt.0.5 .and.
|
|
+ ielite.eq.0) then
|
|
write(*,16)
|
|
endif
|
|
16 format(' WARNING: dangerously high value for pmutmx (ctrl(8));',
|
|
+ /' (Should enforce elitism with ctrl(11)=1.)')
|
|
|
|
if (fdif.lt.0.33d0.and. irep.ne.3) then
|
|
write(*,17)
|
|
endif
|
|
17 format(' WARNING: dangerously low value of fdif (ctrl(9))')
|
|
|
|
if (mod(np,2).gt.0) then
|
|
np=np-1
|
|
write(*,18) np
|
|
endif
|
|
18 format(' WARNING: decreasing population size (ctrl(1)) to np=',i4)
|
|
|
|
return
|
|
end
|
|
c********************************************************************
|
|
subroutine report
|
|
+ (ivrb,ndim,n,np,nd,oldph,fitns,ifit,pmut,ig,nnew)
|
|
c
|
|
c Write generation report to standard output
|
|
c
|
|
implicit none
|
|
|
|
c Input:
|
|
integer np,ifit(np),ivrb,ndim,n,nd,ig,nnew
|
|
double precision oldph(ndim,np),fitns(np),pmut
|
|
c
|
|
c Output: none
|
|
c
|
|
c Local
|
|
double precision bestft,pmutpv
|
|
save bestft,pmutpv
|
|
integer ndpwr,k
|
|
logical rpt
|
|
data bestft,pmutpv /0.0d0,0.0d0/
|
|
c
|
|
rpt=.false.
|
|
|
|
if (pmut.ne.pmutpv) then
|
|
pmutpv=pmut
|
|
rpt=.true.
|
|
endif
|
|
|
|
if (fitns(ifit(np)).ne.bestft) then
|
|
bestft=fitns(ifit(np))
|
|
rpt=.true.
|
|
endif
|
|
|
|
if (rpt .or. ivrb.ge.2) then
|
|
|
|
c Power of 10 to make integer genotypes for display
|
|
ndpwr = idnint(10.d0**nd)
|
|
|
|
write(*,'(/i6,i6,f10.6,4f10.6)') ig,nnew,pmut,
|
|
+ fitns(ifit(np)), fitns(ifit(np-1)), fitns(ifit(np/2))
|
|
do 15 k=1,n
|
|
write(*,'(22x,3i10)')
|
|
+ idnint(ndpwr*oldph(k,ifit(np ))),
|
|
+ idnint(ndpwr*oldph(k,ifit(np-1))),
|
|
+ idnint(ndpwr*oldph(k,ifit(np/2)))
|
|
15 continue
|
|
|
|
endif
|
|
end
|
|
|
|
c**********************************************************************
|
|
c GENETICS MODULE
|
|
c**********************************************************************
|
|
c
|
|
c ENCODE: encodes phenotype into genotype
|
|
c called by: PIKAIA
|
|
c
|
|
c DECODE: decodes genotype into phenotype
|
|
c called by: PIKAIA
|
|
c
|
|
c CROSS: Breeds two offspring from two parents
|
|
c called by: PIKAIA
|
|
c
|
|
c MUTATE: Introduces random mutation in a genotype
|
|
c called by: PIKAIA
|
|
c
|
|
c ADJMUT: Implements variable mutation rate
|
|
c called by: PIKAIA
|
|
c
|
|
c**********************************************************************
|
|
subroutine encode(n,nd,ph,gn)
|
|
c======================================================================
|
|
c encode phenotype parameters into integer genotype
|
|
c ph(k) are x,y coordinates [ 0 < x,y < 1 ]
|
|
c======================================================================
|
|
c
|
|
implicit none
|
|
c
|
|
c Inputs:
|
|
integer n, nd
|
|
double precision ph(n)
|
|
c
|
|
c Output:
|
|
integer gn(n*nd)
|
|
c
|
|
c Local:
|
|
integer ip, i, j, ii
|
|
double precision z
|
|
c
|
|
z=10.0d0**nd
|
|
ii=0
|
|
do 1 i=1,n
|
|
ip=idint(ph(i)*z)
|
|
do 2 j=nd,1,-1
|
|
gn(ii+j)=mod(ip,10)
|
|
ip=ip/10
|
|
2 continue
|
|
ii=ii+nd
|
|
1 continue
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine decode(n,nd,gn,ph)
|
|
c======================================================================
|
|
c decode genotype into phenotype parameters
|
|
c ph(k) are x,y coordinates [ 0 < x,y < 1 ]
|
|
c======================================================================
|
|
c
|
|
implicit none
|
|
c
|
|
c Inputs:
|
|
integer n, nd, gn(n*nd)
|
|
c
|
|
c Output:
|
|
double precision ph(n)
|
|
c
|
|
c Local:
|
|
integer ip, i, j, ii
|
|
double precision z
|
|
c
|
|
z=10.0d0**(-nd)
|
|
ii=0
|
|
do 1 i=1,n
|
|
ip=0
|
|
do 2 j=1,nd
|
|
ip=10*ip+gn(ii+j)
|
|
2 continue
|
|
ph(i)=ip*z
|
|
ii=ii+nd
|
|
1 continue
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine cross(n,nd,pcross,gn1,gn2)
|
|
c======================================================================
|
|
c breeds two parent chromosomes into two offspring chromosomes
|
|
c breeding occurs through crossover. If the crossover probability
|
|
c test yields true (crossover taking place), either one-point or
|
|
c two-point crossover is used, with equal probabilities.
|
|
c
|
|
c Compatibility with version 1.0: To enforce 100% use of one-point
|
|
c crossover, un-comment appropriate line in source code below
|
|
c======================================================================
|
|
c
|
|
implicit none
|
|
c
|
|
c Inputs:
|
|
integer n, nd
|
|
double precision pcross
|
|
c
|
|
c Input/Output:
|
|
integer gn1(n*nd), gn2(n*nd)
|
|
c
|
|
c Local:
|
|
integer i, ispl, ispl2, itmp, t
|
|
c
|
|
c Function
|
|
double precision urand
|
|
external urand
|
|
|
|
|
|
c Use crossover probability to decide whether a crossover occurs
|
|
if (urand().lt.pcross) then
|
|
|
|
c Compute first crossover point
|
|
ispl=idint(urand()*n*nd)+1
|
|
|
|
c Now choose between one-point and two-point crossover
|
|
if (urand().lt.0.5d0) then
|
|
ispl2=n*nd
|
|
else
|
|
ispl2=idint(urand()*n*nd)+1
|
|
c Un-comment following line to enforce one-point crossover
|
|
c ispl2=n*nd
|
|
if (ispl2.lt.ispl) then
|
|
itmp=ispl2
|
|
ispl2=ispl
|
|
ispl=itmp
|
|
endif
|
|
endif
|
|
|
|
c Swap genes from ispl to ispl2
|
|
do 10 i=ispl,ispl2
|
|
t=gn2(i)
|
|
gn2(i)=gn1(i)
|
|
gn1(i)=t
|
|
10 continue
|
|
endif
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine mutate(n,nd,pmut,gn,imut)
|
|
c======================================================================
|
|
c Mutations occur at rate pmut at all gene loci
|
|
c imut=1 Uniform mutation, constant rate
|
|
c imut=2 Uniform mutation, variable rate based on fitness
|
|
c imut=3 Uniform mutation, variable rate based on distance
|
|
c imut=4 Uniform or creep mutation, constant rate
|
|
c imut=5 Uniform or creep mutation, variable rate based on
|
|
c fitness
|
|
c imut=6 Uniform or creep mutation, variable rate based on
|
|
c distance
|
|
c======================================================================
|
|
c
|
|
implicit none
|
|
c
|
|
c Input:
|
|
integer n, nd, imut
|
|
double precision pmut
|
|
c
|
|
c Input/Output:
|
|
integer gn(n*nd)
|
|
c
|
|
c Local:
|
|
integer i,j,k,l,ist,inc,iloc,kk
|
|
|
|
c
|
|
c Function:
|
|
double precision urand
|
|
external urand
|
|
c
|
|
c Decide which type of mutation is to occur
|
|
if(imut.ge.4.and.urand().le.0.5d0)then
|
|
|
|
c CREEP MUTATION OPERATOR
|
|
c Subject each locus to random +/- 1 increment at the rate pmut
|
|
do 1 i=1,n
|
|
do 2 j=1,nd
|
|
if (urand().lt.pmut) then
|
|
c Construct integer
|
|
iloc=(i-1)*nd+j
|
|
inc=idnint(urand())*2-1
|
|
ist=(i-1)*nd+1
|
|
gn(iloc)=gn(iloc)+inc
|
|
c write(*,*) ist,iloc,inc
|
|
c This is where we carry over the one (up to two digits)
|
|
c first take care of decrement below 0 case
|
|
if(inc.lt.0 .and. gn(iloc).lt.0)then
|
|
if(j.eq.1)then
|
|
gn(iloc)=0
|
|
else
|
|
do 3 k=iloc,ist+1,-1
|
|
gn(k)=9
|
|
gn(k-1)=gn(k-1)-1
|
|
if( gn(k-1).ge.0 )goto 4
|
|
3 continue
|
|
c we popped under 0.00000 lower bound; fix it up
|
|
if( gn(ist).lt.0)then
|
|
do 5 l=ist,iloc
|
|
gn(l)=0
|
|
5 continue
|
|
endif
|
|
4 continue
|
|
endif
|
|
endif
|
|
if(inc.gt.0 .and. gn(iloc).gt.9)then
|
|
if(j.eq.1)then
|
|
gn(iloc)=9
|
|
else
|
|
do 6 k=iloc,ist+1,-1
|
|
gn(k)=0
|
|
gn(k-1)=gn(k-1)+1
|
|
if( gn(k-1).le.9 )goto 7
|
|
6 continue
|
|
c we popped over 9.99999 upper bound; fix it up
|
|
if( gn(ist).gt.9 )then
|
|
do 8 l=ist,iloc
|
|
gn(l)=9
|
|
8 continue
|
|
endif
|
|
7 continue
|
|
endif
|
|
endif
|
|
endif
|
|
2 continue
|
|
1 continue
|
|
|
|
else
|
|
|
|
c UNIFORM MUTATION OPERATOR
|
|
c Subject each locus to random mutation at the rate pmut
|
|
do 10 i=1,n*nd
|
|
if (urand().lt.pmut) then
|
|
gn(i)=idint(urand()*10.0d0)
|
|
endif
|
|
10 continue
|
|
endif
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine adjmut(ndim,n,np,oldph,fitns,ifit,pmutmn,pmutmx,
|
|
+ pmut,imut)
|
|
c======================================================================
|
|
c dynamical adjustment of mutation rate;
|
|
c imut=2 or imut=5 : adjustment based on fitness differential
|
|
c between best and median individuals
|
|
c imut=3 or imut=6 : adjustment based on metric distance
|
|
c between best and median individuals
|
|
c======================================================================
|
|
c
|
|
implicit none
|
|
c
|
|
c Input:
|
|
integer n, ndim, np, ifit(np), imut
|
|
double precision oldph(ndim,np), fitns(np), pmutmn, pmutmx
|
|
c
|
|
c Input/Output:
|
|
double precision pmut
|
|
c
|
|
c Local:
|
|
integer i
|
|
double precision rdif, rdiflo, rdifhi, delta
|
|
parameter (rdiflo=0.05d0,rdifhi=0.25d0,delta=1.5d0)
|
|
|
|
if(imut.eq.2.or.imut.eq.5)then
|
|
c Adjustment based on fitness differential
|
|
rdif=abs(fitns(ifit(np))-fitns(ifit(np/2)))/
|
|
+ (fitns(ifit(np))+fitns(ifit(np/2)))
|
|
else if(imut.eq.3.or.imut.eq.6)then
|
|
c Adjustment based on normalized metric distance
|
|
rdif=0.0d0
|
|
do 1 i=1,n
|
|
rdif=rdif+( oldph(i,ifit(np))-oldph(i,ifit(np/2)) )**2
|
|
1 continue
|
|
rdif=dsqrt(rdif)/dble(n)
|
|
endif
|
|
|
|
if(rdif.le.rdiflo)then
|
|
pmut=min(pmutmx,pmut*delta)
|
|
else if(rdif.ge.rdifhi)then
|
|
pmut=max(pmutmn,pmut/delta)
|
|
endif
|
|
|
|
return
|
|
end
|
|
|
|
|
|
c**********************************************************************
|
|
c REPRODUCTION MODULE
|
|
c**********************************************************************
|
|
c
|
|
c SELECT: Parent selection by roulette wheel algorithm
|
|
c called by: PIKAIA
|
|
c
|
|
c RNKPOP: Ranks initial population
|
|
c called by: PIKAIA, NEWPOP
|
|
c
|
|
c GENREP: Inserts offspring into population, for full
|
|
c generational replacement
|
|
c called by: PIKAIA
|
|
c
|
|
c STDREP: Inserts offspring into population, for steady-state
|
|
c reproduction
|
|
c called by: PIKAIA
|
|
c calls: FF
|
|
c
|
|
c NEWPOP: Replaces old generation with new generation
|
|
c called by: PIKAIA
|
|
c calls: FF, RNKPOP
|
|
c
|
|
c**********************************************************************
|
|
subroutine select(np,jfit,fdif,idad)
|
|
c======================================================================
|
|
c Selects a parent from the population, using roulette wheel
|
|
c algorithm with the relative fitnesses of the phenotypes as
|
|
c the "hit" probabilities [see Davis 1991, chap. 1].
|
|
c======================================================================
|
|
c USES: urand
|
|
implicit none
|
|
c
|
|
c Input:
|
|
integer np, jfit(np)
|
|
double precision fdif
|
|
c
|
|
c Output:
|
|
integer idad
|
|
c
|
|
c Local:
|
|
integer np1, i
|
|
double precision dice, rtfit
|
|
c
|
|
c Function:
|
|
double precision urand
|
|
external urand
|
|
c
|
|
c
|
|
np1 = np+1
|
|
dice = urand()*np*np1
|
|
rtfit = 0.0d0
|
|
do 1 i=1,np
|
|
rtfit = rtfit+np1+fdif*(np1-2*jfit(i))
|
|
if (rtfit.ge.dice) then
|
|
idad=i
|
|
goto 2
|
|
endif
|
|
1 continue
|
|
c Assert: loop will never exit by falling through
|
|
|
|
2 return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine rnkpop(n,arrin,indx,rank)
|
|
c======================================================================
|
|
c Calls external sort routine to produce key index and rank order
|
|
c of input array arrin (which is not altered).
|
|
c======================================================================
|
|
c USES: rqsort
|
|
implicit none
|
|
c
|
|
c Input
|
|
integer n
|
|
double precision arrin(n)
|
|
c
|
|
c Output
|
|
integer indx(n),rank(n)
|
|
c
|
|
c Local
|
|
integer i
|
|
c
|
|
c External sort subroutine
|
|
external rqsort
|
|
c
|
|
c
|
|
c Compute the key index
|
|
call rqsort(n,arrin,indx)
|
|
c
|
|
c ...and the rank order
|
|
do 1 i=1,n
|
|
rank(indx(i)) = n-i+1
|
|
1 continue
|
|
return
|
|
end
|
|
|
|
c***********************************************************************
|
|
subroutine genrep(ndim,n,np,ip,ph,newph)
|
|
c=======================================================================
|
|
c full generational replacement: accumulate offspring into new
|
|
c population array
|
|
c=======================================================================
|
|
c
|
|
implicit none
|
|
|
|
c Input:
|
|
integer ndim, n, np, ip
|
|
double precision ph(ndim,2)
|
|
c
|
|
c Output:
|
|
double precision newph(ndim,np)
|
|
c
|
|
c Local:
|
|
integer i1, i2, k
|
|
c
|
|
c
|
|
c Insert one offspring pair into new population
|
|
i1=2*ip-1
|
|
i2=i1+1
|
|
do 1 k=1,n
|
|
newph(k,i1)=ph(k,1)
|
|
newph(k,i2)=ph(k,2)
|
|
1 continue
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine stdrep
|
|
+ (ff,ndim,n,np,irep,ielite,ph,oldph,fitns,ifit,jfit,nnew)
|
|
c======================================================================
|
|
c steady-state reproduction: insert offspring pair into population
|
|
c only if they are fit enough (replace-random if irep=2 or
|
|
c replace-worst if irep=3).
|
|
c======================================================================
|
|
c USES: ff, urand
|
|
implicit none
|
|
c
|
|
c Input:
|
|
integer ndim, n, np, irep, ielite
|
|
double precision ff, ph(ndim,2)
|
|
external ff
|
|
c
|
|
c Input/Output:
|
|
double precision oldph(ndim,np), fitns(np)
|
|
integer ifit(np), jfit(np)
|
|
c
|
|
c Output:
|
|
integer nnew
|
|
|
|
c Local:
|
|
integer i, j, k, i1, if1
|
|
double precision fit
|
|
c
|
|
c External function
|
|
double precision urand
|
|
external urand
|
|
c
|
|
c
|
|
nnew = 0
|
|
do 1 j=1,2
|
|
|
|
c 1. compute offspring fitness (with caller's fitness function)
|
|
fit=ff(n,ph(1,j))
|
|
|
|
c 2. if fit enough, insert in population
|
|
do 20 i=np,1,-1
|
|
if (fit.gt.fitns(ifit(i))) then
|
|
|
|
c make sure the phenotype is not already in the population
|
|
if (i.lt.np) then
|
|
do 5 k=1,n
|
|
if (oldph(k,ifit(i+1)).ne.ph(k,j)) goto 6
|
|
5 continue
|
|
goto 1
|
|
6 continue
|
|
endif
|
|
|
|
c offspring is fit enough for insertion, and is unique
|
|
|
|
c (i) insert phenotype at appropriate place in population
|
|
if (irep.eq.3) then
|
|
i1=1
|
|
else if (ielite.eq.0 .or. i.eq.np) then
|
|
i1=idint(urand()*np)+1
|
|
else
|
|
i1=idint(urand()*(np-1))+1
|
|
endif
|
|
if1 = ifit(i1)
|
|
fitns(if1)=fit
|
|
do 21 k=1,n
|
|
oldph(k,if1)=ph(k,j)
|
|
21 continue
|
|
|
|
c (ii) shift and update ranking arrays
|
|
if (i.lt.i1) then
|
|
|
|
c shift up
|
|
jfit(if1)=np-i
|
|
do 22 k=i1-1,i+1,-1
|
|
jfit(ifit(k))=jfit(ifit(k))-1
|
|
ifit(k+1)=ifit(k)
|
|
22 continue
|
|
ifit(i+1)=if1
|
|
else
|
|
|
|
c shift down
|
|
jfit(if1)=np-i+1
|
|
do 23 k=i1+1,i
|
|
jfit(ifit(k))=jfit(ifit(k))+1
|
|
ifit(k-1)=ifit(k)
|
|
23 continue
|
|
ifit(i)=if1
|
|
endif
|
|
nnew = nnew+1
|
|
goto 1
|
|
endif
|
|
20 continue
|
|
|
|
1 continue
|
|
|
|
return
|
|
end
|
|
|
|
c**********************************************************************
|
|
subroutine newpop
|
|
+ (ff,ielite,ndim,n,np,oldph,newph,ifit,jfit,fitns,nnew)
|
|
c======================================================================
|
|
c replaces old population by new; recomputes fitnesses & ranks
|
|
c======================================================================
|
|
c USES: ff, rnkpop
|
|
implicit none
|
|
c
|
|
c Input:
|
|
integer ndim, np, n, ielite
|
|
double precision ff
|
|
external ff
|
|
c
|
|
c Input/Output:
|
|
double precision oldph(ndim,np), newph(ndim,np)
|
|
c
|
|
c Output:
|
|
integer ifit(np), jfit(np), nnew
|
|
double precision fitns(np)
|
|
c
|
|
c Local:
|
|
integer i, k
|
|
c
|
|
c
|
|
nnew = np
|
|
|
|
c if using elitism, introduce in new population fittest of old
|
|
c population (if greater than fitness of the individual it is
|
|
c to replace)
|
|
if (ielite.eq.1 .and. ff(n,newph(1,1)).lt.fitns(ifit(np))) then
|
|
do 1 k=1,n
|
|
newph(k,1)=oldph(k,ifit(np))
|
|
1 continue
|
|
nnew = nnew-1
|
|
endif
|
|
|
|
c replace population
|
|
do 2 i=1,np
|
|
do 3 k=1,n
|
|
oldph(k,i)=newph(k,i)
|
|
3 continue
|
|
|
|
c get fitness using caller's fitness function
|
|
fitns(i)=ff(n,oldph(1,i))
|
|
2 continue
|
|
|
|
c compute new population fitness rank order
|
|
call rnkpop(np,fitns,ifit,jfit)
|
|
|
|
return
|
|
end
|