pikaia Subroutine

private subroutine pikaia(me, x, f, status)

Optimization (maximization) of user-supplied "fitness" function over n-dimensional parameter space x using a basic genetic algorithm method.

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.

Authors

  • Paul Charbonneau & Barry Knapp (High Altitude Observatory, National Center for Atmospheric Research) Version 1.2 [ 2002 April 3 ]
  • Jacob Williams : 3/8/3015 : Refactoring and some new features.

References

  • Charbonneau, Paul. "An introduction to genetic 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.

Type Bound

pikaia_class

Arguments

Type IntentOptional Attributes Name
class(pikaia_class), intent(inout) :: me
real(kind=wp), intent(inout), dimension(:) :: x

Input - initial guess for solution vector. Output - the "fittest" (optimal) solution found, i.e., the solution which maximizes the fitness function.
real(kind=wp), intent(out) :: f

the (scalar) value of the fitness function at x
integer, intent(out) :: status

an indicator of the success or failure of the call to pikaia (0=success; non-zero=failure)

Calls

proc~~pikaia~~CallsGraph proc~pikaia pikaia_module::pikaia_class%pikaia proc~adjmut pikaia_module::pikaia_class%adjmut proc~pikaia->proc~adjmut proc~cross pikaia_module::pikaia_class%cross proc~pikaia->proc~cross proc~decode pikaia_module::pikaia_class%decode proc~pikaia->proc~decode proc~encode pikaia_module::pikaia_class%encode proc~pikaia->proc~encode proc~func_wrapper pikaia_module::pikaia_class%func_wrapper proc~pikaia->proc~func_wrapper proc~genrep pikaia_module::pikaia_class%genrep proc~pikaia->proc~genrep proc~mutate pikaia_module::pikaia_class%mutate proc~pikaia->proc~mutate proc~newpop pikaia_module::pikaia_class%newpop proc~pikaia->proc~newpop proc~report pikaia_module::pikaia_class%report proc~pikaia->proc~report proc~rninit pikaia_module::pikaia_class%rninit proc~pikaia->proc~rninit proc~rnkpop pikaia_module::pikaia_class%rnkpop proc~pikaia->proc~rnkpop proc~select_parents pikaia_module::pikaia_class%select_parents proc~pikaia->proc~select_parents proc~stdrep pikaia_module::pikaia_class%stdrep proc~pikaia->proc~stdrep proc~urand pikaia_module::pikaia_class%urand proc~pikaia->proc~urand proc~cross->proc~urand proc~mutate->proc~urand proc~newpop->proc~func_wrapper proc~newpop->proc~rnkpop initialize initialize proc~rninit->initialize proc~rqsort pikaia_module::rqsort proc~rnkpop->proc~rqsort proc~select_parents->proc~urand proc~stdrep->proc~urand genrand64_real1 genrand64_real1 proc~urand->genrand64_real1

Called by

proc~~pikaia~~CalledByGraph proc~pikaia pikaia_module::pikaia_class%pikaia proc~solve_with_pikaia pikaia_module::pikaia_class%solve_with_pikaia proc~solve_with_pikaia->proc~pikaia

Source Code

    subroutine pikaia(me,x,f,status)

    implicit none

    !subroutine arguments:
    class(pikaia_class),intent(inout)      :: me
    real(wp),dimension(:),intent(inout)    :: x      !! Input - initial guess for solution vector.
                                                     !! Output - the "fittest" (optimal) solution found,
                                                     !! i.e., the solution which maximizes the fitness function.
    real(wp),intent(out)                   :: f      !! the (scalar) value of the fitness function at x
    integer,intent(out)                    :: status !! an indicator of the success or failure
                                                     !! of the call to pikaia (0=success; non-zero=failure)

    !Local variables
    integer  :: k,ip,ig,ip1,ip2,new,newtot,istart,i_window,j
    real(wp) :: current_best_f, last_best_f, fguess
    logical  :: convergence
    logical  :: use_openmp !! if OpenMP is being used
    real(wp),dimension(me%n,2,me%np/2) :: ph
    real(wp),dimension(me%n,me%np)     :: oldph
    real(wp),dimension(me%n,me%np)     :: newph
    integer,dimension(me%n*me%nd)      :: gn1
    integer,dimension(me%n*me%nd)      :: gn2
    integer,dimension(me%np)           :: ifit
    integer,dimension(me%np)           :: jfit
    real(wp),dimension(me%np)          :: fitns
    real(wp),dimension(me%n)           :: xguess
    real(wp),dimension(2,me%np/2)      :: fits

    real(wp),parameter :: big = huge(1.0_wp)    !! a large number

    !initialize:
    call me%rninit()
    me%bestft   = -big
    me%pmutpv   = -big
    me%pmut     = me%pmuti  !set initial mutation rate (it can change)
    i_window    = 0
    last_best_f = -big
    convergence = .false.
    status      = 0

    ! if OpenMP is being used:
    use_openmp = .false.
!$  use_openmp = .true.

    !Handle the initial guess:
    if (me%initial_guess_frac==0.0_wp) then

        !initial guess not used (totally random population)

        istart = 1  !index to start random population members

    else

        !use the initial guess:

        xguess = x
        do k=1,me%n    !make sure they are all within the [0,1] bounds
            xguess(k) = max( 0.0_wp, min(1.0_wp,xguess(k)) )
        end do
        call me%ff(xguess,fguess)

        !how many elements in the population to set to xguess?:
        ! [at least 1, at most n]
        istart = max(1, min(me%np, int(me%np * me%initial_guess_frac)))

        do k=1,istart
            oldph(:,k) = xguess
            fitns(k)   = fguess
        end do

        istart = istart + 1  !index to start random population members

    end if

    !Compute initial (random but bounded) phenotypes
    do ip=istart,me%np
       do k=1,me%n
          oldph(k,ip) = me%urand()  !from [0,1]
       end do
    end do

    !$omp parallel do private(ip)
    do ip=istart,me%np
       call me%ff(oldph(:,ip),fitns(ip))
    end do
    !$omp end parallel do

    !Rank initial population by fitness order
    call me%rnkpop(fitns,ifit,jfit)

    !Main Generation Loop
    ! This is modified from the original for parallelization.
    ! Note that, now, in a generation, the population is not changed until
    ! all the new members are computed. So only the current members are used
    ! in this process.
    do ig=1,me%ngen

        !Main Population Loop
        newtot = 0
        do ip=1,me%np/2

            !1. pick two parents
            call me%select_parents(jfit,ip1,ip2)

            !2. encode parent phenotypes
            call me%encode(oldph(:,ip1),gn1)
            call me%encode(oldph(:,ip2),gn2)

            !3. breed
            call me%cross(gn1,gn2)
            call me%mutate(gn1)
            call me%mutate(gn2)

            !4. decode offspring genotypes
            call me%decode(gn1,ph(:,1,ip))
            call me%decode(gn2,ph(:,2,ip))

            !5. insert into population
            if (me%irep==1) then
                call me%genrep(ip,ph(:,:,ip),newph)
            else
                if (.not. use_openmp) then
                    ! compute all the fitnesses in the parallel
                    do j = 1, 2
                        ! compute offspring fitness (with caller's fitness function)
                        call me%ff(ph(:,j,ip),fits(j,ip))
                    end do
                    call me%stdrep(ph(:,:,ip),fits(:,ip),oldph,fitns,ifit,jfit,new)
                    newtot = newtot+new
               end if
            end if

        end do

        if (use_openmp) then
            !5. insert into population if not already done
            if (me%irep/=1) then
                ! compute all the fitnesses in the parallel
                !$omp parallel do private(ip)
                do ip=1,me%np/2
                    !$omp parallel do private(j)
                    do j = 1, 2
                        ! compute offspring fitness (with caller's fitness function)
                        call me%ff(ph(:,j,ip),fits(j,ip))
                    end do
                    !$omp end parallel do
                end do
                !$omp end parallel do
                newtot=0
                do ip=1,me%np/2
                    call me%stdrep(ph(:,:,ip),fits(:,ip),oldph,fitns,ifit,jfit,new)
                    newtot = newtot+new
                end do
            end if
        end if
        !End of Main Population Loop

        !if running full generational replacement: swap populations
        if (me%irep==1) call me%newpop(oldph,newph,ifit,jfit,fitns,newtot)

        !adjust mutation rate?
        if (any(me%imut==[2,3,5,6])) call adjmut(me,oldph,fitns,ifit)

        !report this iteration:
        if (me%ivrb>0) call me%report(oldph,fitns,ifit,ig,newtot)

        !report (unscaled) x:
        if (associated(me%iter_f)) &
            call me%iter_f(ig,me%xl+me%del*oldph(1:me%n,ifit(me%np)),fitns(ifit(me%np)))

        !JW additions: add a convergence criteria
        ! [stop if the last convergence_window iterations are all within the convergence_tol]
        current_best_f = fitns(ifit(me%np))    !current iteration best fitness
        if (abs(current_best_f-last_best_f)<=me%convergence_tol) then
            !this solution is within the tol from the previous one
            i_window = i_window + 1    !number of solutions within the convergence tolerance
        else
            i_window = 0    !a significantly better solution was found, reset window
        end if
        if (i_window>=me%convergence_window) then
            convergence = .true.
            exit !exit main loop -> convergence
        end if
        last_best_f = current_best_f    !to compare with next iteration

    end do    !End of Main Generation Loop

    !JW additions:
    if (me%ivrb>0) then
        if (convergence) then
            write(output_unit,'(A)') 'Solution Converged'
        else
            write(output_unit,'(A)') 'Iteration Limit Reached'
        end if
    end if

    !Return best phenotype and its fitness
    x = oldph(1:me%n,ifit(me%np))
    f = fitns(ifit(me%np))

    end subroutine pikaia