! ch13_differential_evolution.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! Bounded minimization by differential evolution: current-to-best/1/bin.
! Edit the objective, bounds, and search settings in the driver.
! No derivatives are required. Console output only.
! Repeat with other seeds; global optimality is not guaranteed.
!
! Compile (Fortran 95):
! gfortran -std=f95 -Wall -Wextra -O2 ch13_differential_evolution.f95 -o ch13_differential_evolution
! Run: ./ch13_differential_evolution (Windows: ch13_differential_evolution.exe)

module differential_evolution_methods
    implicit none
    private
    integer, parameter :: dp=selected_real_kind(15,307)
    public :: dp, initialize_random, differential_evolution
contains
    subroutine initialize_random(seed)
        integer, intent(in) :: seed
        integer, allocatable :: seeds(:)
        integer :: n,i
        call random_seed(size=n)
        allocate(seeds(n))
        do i=1,n
            seeds(i)=1+modulo(modulo(seed,1000000)+37*modulo(i,1000000),1000000)
        end do
        call random_seed(put=seeds)
        deallocate(seeds)
        ! Repeatable on the same compiler/runtime; sequences depend on compiler.
    end subroutine initialize_random
    subroutine differential_evolution(objective,lower,upper,N,K,F,CR,pbest,fbest,evaluations)
        interface
            function objective(p) result(value)
                integer, parameter :: dp=selected_real_kind(15,307)
                real(dp), intent(in) :: p(:)
                real(dp) :: value
            end function objective
        end interface
        real(dp), intent(in) :: lower(:),upper(:),F,CR
        integer, intent(in) :: N,K
        real(dp), intent(out) :: pbest(:),fbest
        integer, intent(out) :: evaluations
        real(dp) :: pop(size(lower),N),next_pop(size(lower),N),values(N),next_values(N)
        real(dp) :: mutant(size(lower)),trial(size(lower)),u,value,width(size(lower))
        integer :: dim,i,j,g,r1,r2,jrand,idx(1)

        dim=size(lower)
        if (dim<1 .or. size(upper)/=dim .or. size(pbest)/=dim) stop 'Mismatched vector dimensions.'
        if (.not. all(finite_value(lower)) .or. .not. all(finite_value(upper))) stop 'Bounds must be finite.'
        width=upper-lower
        if (.not. all(finite_value(width)) .or. any(width<=0.0_dp)) stop 'Invalid bound widths.'
        if (N<4 .or. K<1) stop 'Require at least four members and a positive generation count.'
        if (.not. finite_value(F) .or. F<=0.0_dp .or. F>2.0_dp) stop 'Require 0 < F <= 2.'
        if (.not. finite_value(CR) .or. CR<0.0_dp .or. CR>1.0_dp) stop 'Require 0 <= CR <= 1.'

        do i=1,N
            call random_number(trial)
            pop(:,i)=lower+width*trial
            values(i)=checked_value(objective,pop(:,i))
        end do
        evaluations=N
        do g=1,K
            idx=minloc(values)
            pbest=pop(:,idx(1))
            next_pop=pop
            next_values=values
            ! Use the frozen population and best member throughout a generation.
            do i=1,N
                do
                    call random_number(u)
                    r1=1+int(N*u)
                    if (r1/=i) exit
                end do
                do
                    call random_number(u)
                    r2=1+int(N*u)
                    if (r2/=i .and. r2/=r1) exit
                end do
                mutant=pop(:,i)+F*(pbest-pop(:,i))+F*(pop(:,r1)-pop(:,r2))
                trial=pop(:,i)
                call random_number(u)
                jrand=1+int(dim*u)
                do j=1,dim
                    call random_number(u)
                    ! At least one coordinate comes from the mutant.
                    if (u<CR .or. j==jrand) trial(j)=mutant(j)
                    ! Resample invalid coordinates within the specified bounds.
                    if (.not. finite_value(trial(j)) .or. trial(j)<lower(j) .or. trial(j)>upper(j)) then
                        call random_number(u)
                        trial(j)=lower(j)+width(j)*u
                    end if
                end do
                value=checked_value(objective,trial)
                evaluations=evaluations+1
                if (value<=values(i)) then
                    next_pop(:,i)=trial
                    next_values(i)=value
                end if
            end do
            pop=next_pop
            values=next_values
        end do
        idx=minloc(values)
        pbest=pop(:,idx(1))
        fbest=values(idx(1))
    end subroutine differential_evolution

    function checked_value(objective,p) result(value)
        interface
            function objective(p) result(value)
                integer, parameter :: dp=selected_real_kind(15,307)
                real(dp), intent(in) :: p(:)
                real(dp) :: value
            end function objective
        end interface
        real(dp), intent(in) :: p(:)
        real(dp) :: value
        value=objective(p)
        if (.not. finite_value(value)) stop 'Objective returned a nonfinite value.'
    end function checked_value

    elemental logical function finite_value(x)
        real(dp), intent(in) :: x
        finite_value=abs(x)<=huge(x)
    end function finite_value
end module differential_evolution_methods

module evolution_example
    use differential_evolution_methods, only: dp
    implicit none
contains
    function objective(p) result(value)
        real(dp), intent(in) :: p(:)
        real(dp) :: value,pi
        ! Shifted Rastrigin function: exact minimum is 0 at every p(j)=0.5.
        pi=acos(-1.0_dp)
        value=sum((p-0.5_dp)**2+10.0_dp*(1.0_dp-cos(2.0_dp*pi*(p-0.5_dp))))
    end function objective
end module evolution_example

program ch13_differential_evolution
    use differential_evolution_methods, only: dp, initialize_random, differential_evolution
    use evolution_example, only: objective
    implicit none
    ! Editable settings: one lower and upper bound for each variable.
    integer, parameter :: dim=2,population_size=50,max_generations=400,seed=3
    real(dp), parameter :: lower(dim)=(/-4.5_dp,-4.5_dp/),upper(dim)=(/5.5_dp,5.5_dp/)
    real(dp), parameter :: F=0.9_dp,CR=0.9_dp
    real(dp) :: pbest(dim),fbest
    integer :: evaluations,j
    ! integer :: clock_seed  ! Uncomment for clock-based seeding.

    call initialize_random(seed)
    ! Alternative: comment out the line above, uncomment the clock_seed
    ! declaration above, and uncomment both lines below.
    ! call system_clock(count=clock_seed)
    ! call initialize_random(clock_seed)

    call differential_evolution(objective,lower,upper,population_size,max_generations,F,CR, &
        pbest,fbest,evaluations)
    write(*,'(/,a)') 'Differential evolution'
    do j=1,dim
        write(*,'(a,i0,a,es23.15)') 'Best p(',j,') = ',pbest(j)
    end do
    write(*,'(a,es23.15)') 'Best objective value: ',fbest
    write(*,'(a,i0)') 'Generations:          ',max_generations
    write(*,'(a,i0)') 'Function calls:       ',evaluations
    write(*,'(a)') 'Status: Search budget completed.'
    write(*,'(a)') 'Repeat with other seeds; global optimality is not guaranteed.'
end program ch13_differential_evolution