Program evolve
integer :: chromNo=20, chromLength=32, generations, g, eliteNo !chromLength MUST be 32 or 64, depending on the bits used in rand etc.
real :: crossover=0.8, mutation=5E-3
real, allocatable :: fitness(:), eliteChrom(:,:), searchChrom(:,:)
integer, allocatable :: genePool(:,:), tempGene(:), tempGene2(:), nextGenePool(:,:)
pi = 4.D0*DATAN(1.D0)

!Create and initialize the gene pool:
call random_seed()
allocate(genePool(chromNo,chromLength))
allocate(tempGene(chromLength))
do i = 1, chromNo
        call random_number(rand)
        rand = (rand - 0.5)*10*pi !random number in (-5pi,5pi).
        call encode(rand,tempGene)
        do j = 1, chromLength
                genePool(i,j) = tempGene(j)
        end do
end do

!Now we've set up the gene pool. Lets run our genetic algo on it:
generations=500
!set how many genes are elite, around 10% of chromNo:
eliteNo=2
allocate(fitness(chromNo))     !stores the fitness
allocate(eliteChrom(eliteNo,2)) !stores the fitness and position of elite genes sorted by highest fitness first
allocate(searchChrom(chromNo,2)) !used for sorting, then we just copy the best 10% of this into eliteChrom
allocate(nextGenePool(chromNo,chromLength)) !the next generation
allocate(tempGene2(chromLength))  !we need two temp genes
!Big, outer loop that steps through the generations, g is the loop variable
do g = 1, generations
        !First we determine the fitness for all chromosomes
        do i = 1, chromNo
                do j = 1, chromLength
                         tempGene(j) = genePool(i,j)
                end do
                call decode(tempGene,val)
                !Now val is the value in (-5pi,5pi) represented by chromNo=i
                fitness(i) = exp(-val**2)*cos(val) !this is the fitness function
        end do
        !Now we find the top 10% genes, and simply copy them (elitism).
        !First we copy fitness(:) into searchChrom(:), and set the position-element
        do i = 1, chromNo
                searchChrom(i,1) = fitness(i)
                searchChrom(i,2) = i
        end do
        !Then we sort searchChrom(:), insertion sort, biggest elements last
        do i = 2, chromNo
                a=searchChrom(i,1)
                b=searchChrom(i,2)
                do j=i-1,1,-1
                        if (searchChrom(j,1)<=a) goto 10
                        searchChrom(j+1,1) = searchChrom(j,1)
                        searchChrom(j+1,2) = searchChrom(j,2)
                end do
                        j=0
10              searchChrom(j+1,1)=a
                searchChrom(j+1,2)=b
        end do
        !!!!!!!!!!!!!!!!!!!!!!
        !!THIS IS MORE MAGIC!! 
        !If you want just magic, comment out the following 4 lines. The accuracy
        !of the result will go to hell if you do so, I have no idea why.
        open(unit=666,file="/dev/null")
        write (666, *) searchChrom(:,1)
        write (666, *) searchChrom(:,2)
        close(666)
        !!!!END MORE MAGIC!!!
        !!!!!!!!!!!!!!!!!!!!!
        !Before we start doing anything fancy, we copy the best result into val,
        !and pass this to output at the end.
         call decode(genePool(searchChrom(chromNo,2),:),val)
        !searchChrom(:) is now sorted, and we copy the top 10% (remember biggest elements last)
        do i = 1, eliteNo
                eliteChrom(i,1) = searchChrom(chromNo-i+1,1)
                eliteChrom(i,2) = searchChrom(chromNo-1+1,2)
        end do

        !next, we copy the genes from eliteChrom into nextGenePool
        do i = 1, eliteNo
                intElite = eliteChrom(i,2)
                do j = 1, chromLength
                        nextGenePool(i,j) = genePool(searchChrom(chromNo-i+1,2),j)
                end do
        end do
        !Whew, that's elitism done.
        !Next is straight genetic programming of the remaining 90%.
        !We do {selection, crossover and mutation} enough times to fill up the remaining 90%
        !We were careful with choosing eliteNo, so the 2-step will work.
        do i = eliteNo+1, chromNo-1, 2
                !selecting two genes, with probability proportional to their fitness
                !we jump back here if mutation and/or crossover creates invalid chromosomes.
20              k=0 !comes from after setting of intMax
                l=0
                do while ((k .eq. 0) .AND. (l .eq. 0))
                        call random_number(chrom)
                        chrom = ceiling(chrom*chromNo) !choose a chromosome
                        call random_number(prob)     !probability of using this one
                        prob = prob*searchChrom(1,1) !scale prob with the highest fitness
                        if (fitness(chrom) .gt. prob) then
                                if (k .eq. 0 )then
                                        k = chrom
                                else
                                        l = chrom
                                end if
                        end if
                end do
                !now we know which two chromosomes to use, copy them to tempGene and tempGene2
                do j = 1, chromLength
                        tempGene(j) = genePool(k,j)
                        tempGene2(j) = genePool(l,j)
                end do
                !crossover
                call random_number(rand)
                if (rand .lt. crossover) then
                        call random_number(cross)
                        cross = ceiling(cross*chromLength)
                        do j = cross, chromLength
                                swap = tempGene(j)
                                tempGene(j) = tempGene2(j)
                                tempGene2(j) = swap
                        end do
                end if
                !mutation
                do j = 1, chromLength
                        !with small probability, flip bit j in tempGene(:)
                        call random_number(rand)
                        if (rand .lt. mutation) then
                                if (tempGene(j) .eq. 1) then
                                        tempGene(j) = 0
                                else
                                        tempGene(j) = 1
                                end if
                        end if
                        !do the same for tempGene2(:)
                        call random_number(rand)
                        if (rand .lt. mutation) then
                                if (tempGene2(j) .eq. 1) then
                                        tempGene2(j) = 0
                                else
                                        tempGene2(j) = 1
                                end if
                        end if
                end do
                !check if tempGene and tempGene2 are valid representations
                call intdecode(tempGene,int1)
                call intdecode(tempGene2,int2)
                intMax=2108287141 !maximum integer value, corresponding to 5*pi*2^27
                !if either is above the accepted value, we start selection, crossover and mutation from the beginning.
                if ((abs(int1) .gt. intMax).or.(abs(int2) .gt. intMax)) then
                        goto 20
                end if
                !now we have two new chromosomes that go into nextGenePool(:,:)
                do j = 1, chromLength
                        nextGenePool(i,j) = tempGene(j)
                        nextGenePool(i+1,j) = tempGene2(j)
                end do
        end do
        !Move nextGenePool into genePool, so we're ready for the next iteration
        genePool = nextGenePool
        !Save the best value in resmax:
        resmax = val
end do
write (*, *) "The global maximum value is at"
write (*, *) resmax
write (*, *) "and it is"
write (*, *) exp(-resmax**2)*cos(resmax)


deallocate(tempGene)
deallocate(tempGene2)
deallocate(genePool)
deallocate(nextGenePool)
deallocate(fitness)
deallocate(eliteChrom)
deallocate(searchChrom)

!call test()

End program evolve

!subroutine used to test our encode/decode subroutines.
subroutine test()
real rand, pi
integer tempGene(32)
pi = 4.D0*DATAN(1.D0)
call random_number(rand)
rand = (rand - 0.5)*10*pi !random number in (-5pi,5pi).
write (*, *) rand
call encode(rand,tempGene) !turns tempGene into our binary gene representation
call decode(tempGene,rand) !go back again
write (*, *) rand
write (*, *) tempGene
return
end

!subroutine that encodes our value in (-5pi,5pi) to a 32 element array of 1's and 0's
subroutine encode(num,gene)
integer gene(32), i, enc, bindec
real num, scaled
scaled = num*(2.0**27) !this stretches (-5pi,5pi) almost to the limits of a signed integer. Leaves ~2% of the int's unused.
enc = scaled
i = 32
do while (i .ge. 1)
        bindec = 2**(i-1)
        if ((enc - bindec) .ge. 0) then
                gene(i) = 1
                enc = enc - bindec
        else
                gene(i) = 0
        end if
        i = i - 1
end do
return
end

!subroutine that converts back again
subroutine decode(gene,num)
integer gene(32), i, dec, bindec
real num, scale
i = 32
do while (i .ge. 1)
        bindec=2**(i-1)
        dec = dec + bindec*gene(i)
        i = i - 1
end do
scale = dec
num = scale*7.450580596923828125E-9 !precomputed inverse of 2**27
return
end

!subroutine used when checking whether a mutated gene is allowed or forbidden
subroutine intdecode(gene,dec)
integer gene(32), i, dec, bindec
i = 32
do while (i .ge. 1)
        bindec=2**(i-1)
        dec = dec + bindec*gene(i)
        i = i - 1
end do
return
end