module genet_mod

use simulated_anneal

implicit none

contains

!***********************************************************************
      subroutine pikaia(ff,n,ctrl,x,f,status)
!=======================================================================
!     Optimization (maximization) of user-supplied "fitness" function
!     ff  over n-dimensional parameter space  x  using a basic genetic
!     algorithm method.
!
!     Paul Charbonneau & Barry Knapp
!     High Altitude Observatory
!     National Center for Atmospheric Research
!     Boulder CO 80307-3000
!     <paulchar@hao.ucar.edu>
!     <knapp@hao.ucar.edu>
!
!     Version 1.2   [ 2002 April 3 ]
!
!     Genetic algorithms are heuristic search techniques that
!     incorporate in a computational setting, the biological notion
!     of evolution by means of natural selection.  This subroutine
!     implements the three basic operations of selection, crossover,
!     and mutation, operating on "genotypes" encoded as strings.
!
!     Version 1.2 differs from version 1.0 (December 1995) in that
!     it includes (1) two-point crossover, (2) creep mutation, and
!     (3) dynamical adjustment of the mutation rate based on metric
!     distance in parameter space.
!
!     References:
!
!        Charbonneau, Paul. "An introduction to gemetic algorithms for
!           numerical optimization", NCAR Technical Note TN-450+IA
!           (April 2002)
!
!        Charbonneau, Paul. "Release Notes for PIKAIA 1.2",
!           NCAR Technical Note TN-451+STR (April 2002)
!
!        Charbonneau, Paul, and Knapp, Barry. "A User's Guide
!           to PIKAIA 1.0" NCAR Technical Note TN-418+IA
!           (December 1995)
!
!        Goldberg, David E.  Genetic Algorithms in Search, Optimization,
!           & Machine Learning.  Addison-Wesley, 1989.
!
!        Davis, Lawrence, ed.  Handbook of Genetic Algorithms.
!           Van Nostrand Reinhold, 1991.
!
!=======================================================================
!     USES: ff, urand, setctl, report, rnkpop, select, encode, decode,
!           cross, mutate, genrep, stdrep, newpop, adjmut
      implicit none
 
!     Input:
      integer   n

INTERFACE
   real function ff(n, x) result(f)
!
   IMPLICIT NONE
   INTEGER, INTENT(IN)    :: n
   REAL , INTENT(IN)      :: x(:)
! Local variables
   integer                :: i,j,mini,maxi,nh,nodes,inode,stat
   real ,allocatable      :: obf(:)
   real ,allocatable      :: bb(:,:)
   real ,allocatable      :: xx(:,:),spike(:,:)
   integer, external      :: omp_get_thread_num,omp_get_num_threads
   END function ff
END INTERFACE

!      real      ff
!      external  ff

!
!      o Integer  n  is the parameter space dimension, i.e., the number
!        of adjustable parameters. 
!
!      o Function  ff  is a user-supplied scalar function of n vari-
!        ables, which must have the calling sequence f = ff(n,x), where
!        x is a real parameter array of length n.  This function must
!        be written so as to bound all parameters to the interval [0,1];
!        that is, the user must determine a priori bounds for the para-
!        meter space, and ff must use these bounds to perform the appro-
!        priate scalings to recover true parameter values in the
!        a priori ranges.
!
!        By convention, ff should return higher values for more optimal
!        parameter values (i.e., individuals which are more "fit").
!        For example, in fitting a function through data points, ff
!        could return the inverse of chi**2.
!
!        In most cases initialization code will have to be written
!        (either in a driver or in a separate subroutine) which loads
!        in data values and communicates with ff via one or more labeled
!        common blocks.  An example exercise driver and fitness function
!        are provided in the accompanying file, xpkaia.f.
!
!
!      Input/Output:
       real ctrl(12)
!
!      o Array  ctrl  is an array of control flags and parameters, to
!        control the genetic behavior of the algorithm, and also printed
!        output.  A default value will be used for any control variable
!        which is supplied with a value less than zero.  On exit, ctrl
!        contains the actual values used as control variables.  The
!        elements of ctrl and their defaults are:
!
!           ctrl( 1) - number of individuals in a population (default
!                      is 100)
!           ctrl( 2) - number of generations over which solution is
!                      to evolve (default is 500)
!           ctrl( 3) - number of significant digits (i.e., number of
!                      genes) retained in chromosomal encoding (default
!                      is 6)  (Note: This number is limited by the
!                      machine floating point precision.  Most 32-bit
!                      floating point representations have only 6 full
!                      digits of precision.  To achieve greater preci-
!                      sion this routine could be converted to double
!                      precision, but note that this would also require
!                      a double precision random number generator, which
!                      likely would not have more than 9 digits of
!                      precision if it used 4-byte integers internally.)
!           ctrl( 4) - crossover probability; must be  <= 1.0 (default
!                      is 0.85). If crossover takes place, either one
!                      or two splicing points are used, with equal
!                      probabilities
!           ctrl( 5) - mutation mode; 1/2/3/4/5 (default is 2)
!                      1=one-point mutation, fixed rate
!                      2=one-point, adjustable rate based on fitness 
!                      3=one-point, adjustable rate based on distance
!                      4=one-point+creep, fixed rate
!                      5=one-point+creep, adjustable rate based on fitness
!                      6=one-point+creep, adjustable rate based on distance
!           ctrl( 6) - initial mutation rate; should be small (default
!                      is 0.005) (Note: the mutation rate is the proba-
!                      bility that any one gene locus will mutate in
!                      any one generation.)
!           ctrl( 7) - minimum mutation rate; must be >= 0.0 (default
!                      is 0.0005)
!           ctrl( 8) - maximum mutation rate; must be <= 1.0 (default
!                      is 0.25)
!           ctrl( 9) - relative fitness differential; range from 0
!                      (none) to 1 (maximum).  (default is 1.)
!           ctrl(10) - reproduction plan; 1/2/3=Full generational
!                      replacement/Steady-state-replace-random/Steady-
!                      state-replace-worst (default is 3)
!           ctrl(11) - elitism flag; 0/1=off/on (default is 0)
!                      (Applies only to reproduction plans 1 and 2)
!           ctrl(12) - printed output 0/1/2=None/Minimal/Verbose
!                      (default is 0)
!
!
!     Output:
      real      x(n), f
      integer   status
!
!      o Array  x(1:n)  is the "fittest" (optimal) solution found,
!         i.e., the solution which maximizes fitness function ff
!
!      o Scalar  f  is the value of the fitness function at x
!
!      o Integer  status  is an indicator of the success or failure
!         of the call to pikaia (0=success; non-zero=failure)
!
!
!     Constants
      integer   NMAX, PMAX, DMAX
!      parameter (NMAX = ntot, PMAX = 128, DMAX = 6)
!
!      o NMAX is the maximum number of adjustable parameters
!        (n <= NMAX)
!
!      o PMAX is the maximum population (ctrl(1) <= PMAX)
!
!      o DMAX is the maximum number of Genes (digits) per Chromosome
!        segement (parameter) (ctrl(3) <= DMAX)
!
!
!     Local variables
      integer        np, nd, ngen, imut, irep, ielite, ivrb, k, ip, ig,&
                     ip1, ip2, new, newtot
      real           pcross, pmut, pmutmn, pmutmx, fdif
!
      real,allocatable           ::  ph(:,:), oldph(:,:), newph(:,:), fitns(:)
!
      integer,allocatable        :: gn1(:), gn2(:)
      integer,allocatable        :: ifit(:), jfit(:)

!
!     User-supplied uniform random number generator
!      real           urand
!      external       urand
!
!     Function urand should not take any arguments.  If the user wishes
!     to be able to initialize urand, so that the same sequence of
!     random numbers can be repeated, this capability could be imple-
!     mented with a separate subroutine, and called from the users
!     driver program.  An example urand function (and initialization
!     subroutine) which uses the function ran0 (the "minimal standard"
!     random number generator of Park and Miller [Comm. ACM 31, 1192-
!     1201, Oct 1988; Comm. ACM 36 No. 7, 105-110, July 1993]) is
!     provided.
!
!
!     Set control variables from input and defaults

      NMAX = ntot+1
      PMAX = 128
      DMAX = 6

      allocate(ph(NMAX,2), oldph(NMAX,PMAX), newph(NMAX,PMAX), fitns(PMAX))
      allocate(gn1(NMAX*DMAX), gn2(NMAX*DMAX),ifit(PMAX), jfit(PMAX))

      call setctl (ctrl,n,np,ngen,nd,pcross,pmutmn,pmutmx,pmut,imut,fdif,irep,ielite,ivrb,status)

!!!!!!!!!!!!!!!!
! DANGER       !
!              !
      np=pmax  !
!!!!!!!!!!!!!!!!
      if (status .ne. 0) then
         write(*,*) ' Control vector (ctrl) argument(s) invalid'
         return
      endif

!     Make sure locally-dimensioned arrays are big enough
      if (n.gt.NMAX .or. np.gt.PMAX .or. nd.gt.DMAX) then
         write(*,*)&
            ' Number of parameters, population, or genes too large'
         status = -1
         return
      endif
 
!     Compute initial (random but bounded) phenotypes
      do ip=1,np
         do k=1,n
            oldph(k,ip)=urand()
!  write(0,*) ip,k

         end do
! Changed from ff(n,oldph(1,ip)) to ff(n,oldph(:,ip))
! ff needs a vector as input
         fitns(ip) = ff(n,oldph(:,ip))

      end do

!     Rank initial population by fitness order
      call rnkpop(np,fitns,ifit,jfit)
write(0,*) "passei 02"
 
!     Main Generation Loop
      do 10 ig=1,ngen
 
!        Main Population Loop
         newtot=0
         do 20 ip=1,np/2
 
!           1. pick two parents
            call select(np,jfit,fdif,ip1)
   21       call select(np,jfit,fdif,ip2)
write(0,*) "passei 03"

            if (ip1.eq.ip2) goto 21
 
!           2. encode parent phenotypes
            call encode(n,nd,oldph(1,ip1),gn1)
            call encode(n,nd,oldph(1,ip2),gn2)
 write(0,*) "passei 04"

!           3. breed
            call cross(n,nd,pcross,gn1,gn2)
            call mutate(n,nd,pmut,gn1,imut)
            call mutate(n,nd,pmut,gn2,imut)
 write(0,*) "passei 05"

!           4. decode offspring genotypes
            call decode(n,nd,gn1,ph(1,1))
            call decode(n,nd,gn2,ph(1,2))
 write(0,*) "passei 06"

!           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
 write(0,*) "passei 07"

!        End of Main Population Loop
   20    continue
 
!        if running full generational replacement: swap populations
         if (irep.eq.1)&
            call newpop(ff,ielite,NMAX,n,np,oldph,newph,&
               ifit,jfit,fitns,newtot)
 write(0,*) "passei 08"

!        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)
!
         if (ivrb.gt.0) call report&
            (ivrb,NMAX,n,np,nd,oldph,fitns,ifit,pmut,ig,newtot)
 
!     End of Main Generation Loop
   10 continue
write(0,*) "passei 09"
!
!     Return best phenotype and its fitness
      do 30 k=1,n
         x(k) = oldph(k,ifit(np))
   30 continue
      f = fitns(ifit(np))
!
      end subroutine
!*********************************************************************
real function objfcn3(n, x) result(f)
!
IMPLICIT NONE
INTEGER, INTENT(IN)    :: n
REAL , INTENT(IN)      :: x(:)
! Local variables
integer                :: i,j,mini,maxi,nh,nodes,inode,stat
real ,allocatable      :: obf(:)
real ,allocatable      :: bb(:,:)
real ,allocatable      :: xx(:,:),spike(:,:)
integer, external      :: omp_get_thread_num,omp_get_num_threads

!$OMP PARALLEL
nodes = omp_get_num_threads()
!$OMP END PARALLEL

allocate(obf(nodes),bb(nw,nodes),xx(2*nw-1,nodes),spike(2*nw-1,nodes))
obf=0.0
f = 0.0
nh=11
!$OMP PARALLEL DO PRIVATE (inode,&
                               j,&
                            mini,&
                            maxi,&
                               i) schedule(dynamic,1)  num_threads(nodes) 
                               
do j = 1, nx*nz
   inode = omp_get_thread_num()+1
   bb(:,inode)=x((j-1)*nw+1:j*nw)
   call tcai1_lop2 ( dble(xx(:,inode)), dble(x((j-1)*nw+1:j*nw)),dble(bb(:,inode)))
   xx(:,inode)=xx(:,inode)/maxval(abs(xx(:,inode)))
   spike(:,inode)=0.0
   where (xx(:,inode) .eq. 1.)
      spike(:,inode)=1.
   end where
   obf(inode) = obf(inode) + sum((xx(:,inode) - spike(:,inode))**2)
!   mini=(max(1,int(j/nz)+1-nh)-1)*nz+1
!write(*, '(2i6)') j,inode
   mini=(max(1,int(j/nz)+1)-1)*nz+1
   maxi=(min(nx,int(j/nz)+1+nh)-1)*nz
   do i = mini,maxi
      if (i.ne.j) then
         call tcai1_lop2 ( dble(xx(:,inode)), dble(x((i-1)*nw+1:i*nw)),dble(bb(:,inode)))
         xx(:,inode)=xx(:,inode)/maxval(abs(xx(:,inode)))
         obf(inode) = obf(inode) + sum(xx(:,inode)**2)
      end if
   end do
end do
!$OMP END PARALLEL DO

f=-sum(obf)

write(*, '("  The obj.fcn: f = ",g25.8)') f
RETURN
END function objfcn3

!*********************************************************************
      real function urand() result(urand1)
!=====================================================================
!     Return the next pseudo-random deviate from a sequence which is
!     uniformly distributed in the interval [0,1]
!
!     Uses the function ran0, the "minimal standard" random number
!     generator of Park and Miller (Comm. ACM 31, 1192-1201, Oct 1988;
!     Comm. ACM 36 No. 7, 105-110, July 1993).
!=====================================================================
      implicit none
!
!     Input - none
!
!     Output
!      real     urand1
!
!     Local
      integer  iseed
!      real     ran0
!      external ran0
!
!     Common block to make iseed visible to rninit (and to save
!     it between calls)
      common /rnseed/ iseed
!
      urand1 = ran0( iseed )
      return
      end function urand
!*********************************************************************
      subroutine rninit( seed )
!=====================================================================
!     Initialize random number generator urand with given seed
!=====================================================================
      implicit none
!
!     Input
      integer seed
!
!     Output - none
!
!     Local
      integer iseed
!
!     Common block to communicate with urand
      common /rnseed/ iseed
!
!     Set the seed value
      iseed = seed
      if(iseed.le.0) iseed=123456
      return
      end subroutine 
!*********************************************************************
      real function ran0( seed ) result(ran1)
!=====================================================================
!     "Minimal standard" pseudo-random number generator of Park and
!     Miller.  Returns a uniform random deviate r s.t. 0 < r < 1.0.
!     Set seed to any non-zero integer value to initialize a sequence,
!     then do not change seed between calls for successive deviates
!     in the sequence.
!
!     References:
!        Park, S. and Miller, K., "Random Number Generators: Good Ones
!           are Hard to Find", Comm. ACM 31, 1192-1201 (Oct. 1988)
!        Park, S. and Miller, K., in "Remarks on Choosing and Imple-
!           menting Random Number Generators", Comm. ACM 36 No. 7,
!           105-110 (July 1993)
!=====================================================================
! *** Declaration section ***
!
      implicit none
!
!     Input/Output:
      integer seed
!
!     Output:
!      real ran1
!
!     Constants:
      integer A,M,Q,R
      parameter (A=48271,M=2147483647,Q=44488,R=3399)
      real SCALE,EPS,RNMX
      parameter (SCALE=1./M,EPS=1.2e-7,RNMX=1.-EPS)
!
!     Local:
      integer j
!
! *** Executable section ***
!
      j = seed/Q
      seed = A*(seed-j*Q)-R*j
      if (seed .lt. 0) seed = seed+M
      ran1 = min(seed*SCALE,RNMX)
      return
      end function ran0
!**********************************************************************
      subroutine rqsort(n,a,p)
!======================================================================
!     Return integer array p which indexes array a in increasing order.
!     Array a is not disturbed.  The Quicksort algorithm is used.
!
!     B. G. Knapp, 86/12/23
!
!     Reference: N. Wirth, Algorithms and Data Structures,
!     Prentice-Hall, 1986
!======================================================================
      implicit none

!     Input:
      integer   n
      real      a(n)

!     Output:
      integer   p(n)

!     Constants
      integer   LGN, Q
      parameter (LGN=32, Q=11)
!        (LGN = log base 2 of maximum n;
!         Q = smallest subfile to use quicksort on)

!     Local:
      real      x
      integer   stackl(LGN),stackr(LGN),s,t,l,m,r,i,j

!     Initialize the stack
      stackl(1)=1
      stackr(1)=n
      s=1

!     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

!           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
!           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
!           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
!           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 subroutine

!********************************************************************
      subroutine setctl(ctrl,n,np,ngen,nd,pcross,pmutmn,pmutmx,pmut,imut,&
          fdif,irep,ielite,ivrb,status)
!===================================================================
!     Set control variables and flags from input and defaults
!===================================================================
      implicit none
!
!     Input
      integer  n
!
!     Input/Output
      real     ctrl(12)
!
!     Output
      integer  np, ngen, nd, imut, irep, ielite, ivrb, status
      real     pcross, pmutmn, pmutmx, pmut, fdif
!
!     Local
      integer  i
      real     DFAULT(12)
      save     DFAULT
      data     DFAULT /100,500,5,.85,2,.005,.0005,.25,1,1,1,0/
!
      do 1 i=1,12
         if (ctrl(i).lt.0.) 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)
      status = 0
!
!     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
 
!     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)
         status = 5
      endif
   10 format(' ERROR: illegal value for imut (ctrl(5))')
 
      if (fdif.gt.1.) then
         write(*,11)
         status = 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)
         status = 10
      endif
   12 format(' ERROR: illegal value for irep (ctrl(10))')
 
      if (pcross.gt.1.0 .or. pcross.lt.0.) then
         write(*,13)
         status = 4
      endif
   13 format(' ERROR: illegal value for pcross (ctrl(4))')
 
      if (ielite.ne.0 .and. ielite.ne.1) then
         write(*,14)
         status = 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.33 .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 subroutine
!********************************************************************
      subroutine report(ivrb,ndim,n,np,nd,oldph,fitns,ifit,pmut,ig,nnew)
!
!     Write generation report to standard output
!
      implicit none
 
!     Input:
      integer np,ifit(np),ivrb,ndim,n,nd,ig,nnew
      real oldph(ndim,np),fitns(np),pmut
!
!     Output: none
!
!     Local
      real bestft,pmutpv
      save bestft,pmutpv
      integer ndpwr,k
      logical rpt
      data bestft,pmutpv /0,0/
!
      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
 
!        Power of 10 to make integer genotypes for display
         ndpwr = nint(10.**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)')&
               nint(ndpwr*oldph(k,ifit(np  ))),&
               nint(ndpwr*oldph(k,ifit(np-1))),&
               nint(ndpwr*oldph(k,ifit(np/2)))
   15    continue
 
      endif
      end subroutine

!**********************************************************************
!                         GENETICS MODULE
!**********************************************************************
!
!     ENCODE:    encodes phenotype into genotype
!                called by: PIKAIA
!
!     DECODE:    decodes genotype into phenotype
!                called by: PIKAIA
!
!     CROSS:     Breeds two offspring from two parents
!                called by: PIKAIA
!
!     MUTATE:    Introduces random mutation in a genotype
!                called by: PIKAIA
!
!     ADJMUT:    Implements variable mutation rate
!                called by: PIKAIA
!
!**********************************************************************
      subroutine encode(n,nd,ph,gn)
!======================================================================
!     encode phenotype parameters into integer genotype
!     ph(k) are x,y coordinates [ 0 < x,y < 1 ]
!======================================================================
!
      implicit none
!
!     Inputs:
      integer   n, nd
      real      ph(n)
!
!     Output:
      integer   gn(n*nd)
!
!     Local:
      integer   ip, i, j, ii
      real      z
!
      z=10.**nd
      ii=0
      do 1 i=1,n
         ip=int(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 subroutine
 
!**********************************************************************
      subroutine decode(n,nd,gn,ph)
!======================================================================
!     decode genotype into phenotype parameters
!     ph(k) are x,y coordinates [ 0 < x,y < 1 ]
!======================================================================
!
      implicit none
!
!     Inputs:
      integer   n, nd, gn(n*nd)
!
!     Output:
      real      ph(n)
!
!     Local:
      integer   ip, i, j, ii
      real      z
!
      z=10.**(-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 subroutine
 
!**********************************************************************
      subroutine cross(n,nd,pcross,gn1,gn2)
!======================================================================
!     breeds two parent chromosomes into two offspring chromosomes
!     breeding occurs through crossover. If the crossover probability
!     test yields true (crossover taking place), either one-point or
!     two-point crossover is used, with equal probabilities.
!
!     Compatibility with version 1.0: To enforce 100% use of one-point
!     crossover, un-comment appropriate line in source code below
!======================================================================
!
      implicit none
!
!     Inputs:
      integer        n, nd
      real           pcross
!
!     Input/Output:
      integer        gn1(n*nd), gn2(n*nd)
!
!     Local:
      integer        i, ispl, ispl2, itmp, t
!
!     Function
!      real           urand
!      external       urand
 
 
!     Use crossover probability to decide whether a crossover occurs
      if (urand().lt.pcross) then

!        Compute first crossover point
         ispl=int(urand()*n*nd)+1
 
!        Now choose between one-point and two-point crossover 
         if (urand().lt.0.5) then
            ispl2=n*nd
         else
            ispl2=int(urand()*n*nd)+1
!           Un-comment following line to enforce one-point crossover
!           ispl2=n*nd
            if (ispl2.lt.ispl) then
               itmp=ispl2
               ispl2=ispl
               ispl=itmp
            endif
         endif 
 
!        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 subroutine
 
!**********************************************************************
      subroutine mutate(n,nd,pmut,gn,imut)
!======================================================================
!     Mutations occur at rate pmut at all gene loci
!        imut=1    Uniform mutation, constant rate
!        imut=2    Uniform mutation, variable rate based on fitness
!        imut=3    Uniform mutation, variable rate based on distance
!        imut=4    Uniform or creep mutation, constant rate
!        imut=5    Uniform or creep mutation, variable rate based on
!                  fitness
!        imut=6    Uniform or creep mutation, variable rate based on
!                  distance
!======================================================================
!
      implicit none
!
!     Input:
      integer        n, nd, imut
      real           pmut
!
!     Input/Output:
      integer        gn(n*nd)
!
!     Local:
      integer        i,j,k,l,ist,inc,loc,kk

!
!     Function:
!      real           urand
!      external       urand
!
!     Decide which type of mutation is to occur
      if(imut.ge.4.and.urand().le.0.5)then

!     CREEP MUTATION OPERATOR
!     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
!     Construct integer
                  loc=(i-1)*nd+j
                  inc=nint ( urand() )*2-1
                  ist=(i-1)*nd+1
                  gn(loc)=gn(loc)+inc
!                 write(*,*) ist,loc,inc
!     This is where we carry over the one (up to two digits)
!     first take care of decrement below 0 case
		  if(inc.lt.0 .and. gn(loc).lt.0)then
		     if(j.eq.1)then
                        gn(loc)=0
                     else
                        do 3 k=loc,ist+1,-1
                           gn(k)=9
                           gn(k-1)=gn(k-1)-1
                        if( gn(k-1).ge.0 )goto 4
    3                   continue
!    we popped under 0.00000 lower bound; fix it up
                        if( gn(ist).lt.0.)then
                           do 5 l=ist,loc
                              gn(l)=0
    5                      continue
                        endif
    4                   continue
                     endif
                  endif
                  if(inc.gt.0 .and. gn(loc).gt.9)then
                     if(j.eq.1)then
                        gn(loc)=9
                     else
                        do 6 k=loc,ist+1,-1
                           gn(k)=0
                           gn(k-1)=gn(k-1)+1
                           if( gn(k-1).le.9 )goto 7
    6                   continue
!                       we popped over 9.99999 upper bound; fix it up
                       if( gn(ist).gt.9 )then
                           do 8 l=ist,loc
                              gn(l)=9
    8                      continue
                        endif
    7                   continue
                     endif
                  endif
               endif
    2       continue
    1    continue

      else

!     UNIFORM MUTATION OPERATOR
!     Subject each locus to random mutation at the rate pmut
         do 10 i=1,n*nd
            if (urand().lt.pmut) then
               gn(i)=int(urand()*10.)
            endif
   10    continue
      endif

      return
      end subroutine
 
!**********************************************************************
      subroutine adjmut(ndim,n,np,oldph,fitns,ifit,pmutmn,pmutmx,pmut,imut)
!======================================================================
!     dynamical adjustment of mutation rate;
!        imut=2 or imut=5 : adjustment based on fitness differential
!                           between best and median individuals
!        imut=3 or imut=6 : adjustment based on metric distance
!                           between best and median individuals
!======================================================================
!
      implicit none
!
!     Input:
      integer        n, ndim, np, ifit(np), imut
      real           oldph(ndim,np), fitns(np), pmutmn, pmutmx
!
!     Input/Output:
      real           pmut
!
!     Local:
      integer        i
      real           rdif, rdiflo, rdifhi, delta
      parameter      (rdiflo=0.05, rdifhi=0.25, delta=1.5)

      if(imut.eq.2.or.imut.eq.5)then
!     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
!     Adjustment based on normalized metric distance
         rdif=0.
         do 1 i=1,n
            rdif=rdif+( oldph(i,ifit(np))-oldph(i,ifit(np/2)) )**2
    1    continue
         rdif=sqrt( rdif ) / float(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 subroutine


!**********************************************************************
!                       REPRODUCTION MODULE
!**********************************************************************
!
!     SELECT:   Parent selection by roulette wheel algorithm
!               called by: PIKAIA
!
!     RNKPOP:   Ranks initial population
!               called by: PIKAIA, NEWPOP
!
!     GENREP:   Inserts offspring into population, for full
!               generational replacement
!               called by: PIKAIA
!
!     STDREP:   Inserts offspring into population, for steady-state
!               reproduction
!               called by: PIKAIA
!               calls:     FF
!
!     NEWPOP:   Replaces old generation with new generation
!               called by: PIKAIA
!               calls:     FF, RNKPOP
!
!**********************************************************************
      subroutine select(np,jfit,fdif,idad)
!======================================================================
!     Selects a parent from the population, using roulette wheel
!     algorithm with the relative fitnesses of the phenotypes as
!     the "hit" probabilities [see Davis 1991, chap. 1].
!======================================================================
!     USES: urand
      implicit none
!
!     Input:
      integer        np, jfit(np)
      real           fdif
!
!     Output:
      integer        idad
!
!     Local:
      integer        np1, i
      real           dice, rtfit
!
!     Function:
!      real           urand
!      external       urand
!
!
      np1 = np+1
      dice = urand()*np*np1
      rtfit = 0.
      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
!     Assert: loop will never exit by falling through
 
    2 return
      end subroutine
 
!**********************************************************************
      subroutine rnkpop(n,arrin,indx,rank)
!======================================================================
!     Calls external sort routine to produce key index and rank order
!     of input array arrin (which is not altered).
!======================================================================
!     USES: rqsort
      implicit none
!
!     Input
      integer    n
      real       arrin(n)
!
!     Output
      integer    indx(n),rank(n)
!
!     Local
      integer    i
!
!     External sort subroutine
!      external rqsort
!
!
!     Compute the key index
      call rqsort(n,arrin,indx)
!
!     ...and the rank order
      do 1 i=1,n
         rank(indx(i)) = n-i+1
    1 continue
      return
      end subroutine
 
!***********************************************************************
      subroutine genrep(ndim,n,np,ip,ph,newph)
!=======================================================================
!     full generational replacement: accumulate offspring into new
!     population array
!=======================================================================
!
      implicit none
 
!     Input:
      integer        ndim, n, np, ip
      real           ph(ndim,2)
!
!     Output:
      real           newph(ndim,np)
!
!     Local:
      integer        i1, i2, k
!
!
!     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 subroutine
 
!**********************************************************************
      subroutine stdrep(ff,ndim,n,np,irep,ielite,ph,oldph,fitns,ifit,jfit,nnew)
!======================================================================
!     steady-state reproduction: insert offspring pair into population
!     only if they are fit enough (replace-random if irep=2 or
!     replace-worst if irep=3).
!======================================================================
!     USES: ff, urand
      implicit none
!
!     Input:
      integer        ndim, n, np, irep, ielite
      real           ff, ph(ndim,2)
      external       ff
!
!     Input/Output:
      real           oldph(ndim,np), fitns(np)
      integer        ifit(np), jfit(np)
!
!     Output:
      integer        nnew
 
!     Local:
      integer        i, j, k, i1, if1
      real           fit
!
!     External function
!      real           urand
!      external       urand
!
!
      nnew = 0
      do 1 j=1,2
 
!        1. compute offspring fitness (with callers fitness function)
         fit=ff(n,ph(1,j))
 
!        2. if fit enough, insert in population
         do 20 i=np,1,-1
            if (fit.gt.fitns(ifit(i))) then
 
!              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
 
!              offspring is fit enough for insertion, and is unique
 
!              (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=int(urand()*np)+1
               else
                  i1=int(urand()*(np-1))+1
               endif
               if1 = ifit(i1)
               fitns(if1)=fit
               do 21 k=1,n
                  oldph(k,if1)=ph(k,j)
   21          continue
 
!              (ii) shift and update ranking arrays
               if (i.lt.i1) then
 
!                 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
 
!                 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 subroutine
 
!**********************************************************************
      subroutine newpop(ff,ielite,ndim,n,np,oldph,newph,ifit,jfit,fitns,nnew)
!======================================================================
!     replaces old population by new; recomputes fitnesses & ranks
!======================================================================
!     USES: ff, rnkpop
      implicit none
!
!     Input:
      integer        ndim, np, n, ielite
      real           ff
      external       ff
!
!     Input/Output:
      real           oldph(ndim,np), newph(ndim,np)
!
!     Output:
      integer        ifit(np), jfit(np), nnew
      real           fitns(np)
!
!     Local:
      integer        i, k
!
!
      nnew = np
 
!     if using elitism, introduce in new population fittest of old
!     population (if greater than fitness of the individual it is
!     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
 
!     replace population
      do 2 i=1,np
         do 3 k=1,n
            oldph(k,i)=newph(k,i)
    3    continue
 
!        get fitness using callers fitness function
         fitns(i)=ff(n,oldph(1,i))
    2 continue
 
!     compute new population fitness rank order
      call rnkpop(np,fitns,ifit,jfit)
 
      return
      end subroutine
 
end module
