pikaia Subroutine

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

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)
logical, intent(in), optional :: omp

if OpenMP is being used

Calls

proc~~pikaia~~CallsGraph proc~pikaia pikaia_class%pikaia proc~adjmut pikaia_class%adjmut proc~pikaia->proc~adjmut proc~cross pikaia_class%cross proc~pikaia->proc~cross proc~decode pikaia_class%decode proc~pikaia->proc~decode proc~encode pikaia_class%encode proc~pikaia->proc~encode proc~func_wrapper pikaia_class%func_wrapper proc~pikaia->proc~func_wrapper proc~genrep pikaia_class%genrep proc~pikaia->proc~genrep proc~mutate pikaia_class%mutate proc~pikaia->proc~mutate proc~newpop pikaia_class%newpop proc~pikaia->proc~newpop proc~report pikaia_class%report proc~pikaia->proc~report proc~rninit pikaia_class%rninit proc~pikaia->proc~rninit proc~rnkpop pikaia_class%rnkpop proc~pikaia->proc~rnkpop proc~select_parents pikaia_class%select_parents proc~pikaia->proc~select_parents proc~stdrep pikaia_class%stdrep proc~pikaia->proc~stdrep proc~urand pikaia_class%urand proc~pikaia->proc~urand proc~cross->proc~urand proc~mutate->proc~urand proc~newpop->proc~func_wrapper proc~newpop->proc~rnkpop none~initialize mt19937%initialize proc~rninit->none~initialize proc~rqsort rqsort proc~rnkpop->proc~rqsort proc~select_parents->proc~urand proc~stdrep->proc~urand proc~genrand64_real1 mt19937%genrand64_real1 proc~urand->proc~genrand64_real1 proc~init_by_array64 mt19937%init_by_array64 none~initialize->proc~init_by_array64 proc~init_genrand64 mt19937%init_genrand64 none~initialize->proc~init_genrand64 proc~init_genrand64_i4 mt19937%init_genrand64_i4 none~initialize->proc~init_genrand64_i4 proc~genrand64_int64 mt19937%genrand64_int64 proc~genrand64_real1->proc~genrand64_int64 proc~genrand64_int64->proc~init_genrand64 proc~init_by_array64->proc~init_genrand64 proc~init_genrand64_i4->none~initialize

Called by

proc~~pikaia~~CalledByGraph proc~pikaia pikaia_class%pikaia proc~solve_with_pikaia pikaia_class%solve_with_pikaia proc~solve_with_pikaia->proc~pikaia program~test_algen test_algen program~test_algen->proc~solve_with_pikaia

Source Code

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

    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)
    logical, intent(in), optional          :: omp !! if OpenMP is being used

    !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 = .true.
    if ( present( omp ) ) use_openmp = omp

    !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