! ch13_genetic_algorithm.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! Bounded minimization by a real-valued genetic algorithm.
! Edit the objective, bounds, and search settings in the driver.
! Tournament selection compares objective values directly (smaller is better).
! No derivatives required. Repeat with other seeds and search budgets.
! The best point found is not proof of a global minimum. Console output only.
!
! Compile (Fortran 95):
! gfortran -std=f95 -Wall -Wextra -O2 ch13_genetic_algorithm.f95 -o ch13_genetic_algorithm
! Run: ./ch13_genetic_algorithm (Windows: ch13_genetic_algorithm.exe)

module genetic_algorithm_methods
    implicit none
    private
    integer, parameter :: dp=selected_real_kind(15,307)
    public :: dp, initialize_random, genetic_search
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 genetic_search(f,lower,upper,N,K,pc,pm,sigma,nt,pbest,fbest,evaluations)
        interface
            function f(p) result(value)
                integer, parameter :: dp=selected_real_kind(15,307)
                real(dp), intent(in) :: p(:)
                real(dp) :: value
            end function f
        end interface
        real(dp), intent(in) :: lower(:),upper(:),pc,pm,sigma(:)
        integer, intent(in) :: N,K,nt
        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) :: p1(size(lower)),p2(size(lower)),children(size(lower),2),child(size(lower)),u
        integer :: dim,i,j,g,count,idx(1)

        dim=size(lower)
        if (dim<1 .or. size(upper)/=dim .or. size(sigma)/=dim .or. size(pbest)/=dim) &
            stop 'Bounds, mutation widths, and result must have matching dimensions.'
        if (.not. all(finite_value(lower)) .or. .not. all(finite_value(upper))) stop 'Bounds must be finite.'
        if (.not. all(finite_value(upper-lower)) .or. any(upper<=lower)) stop 'Invalid bound widths.'
        if (.not. all(finite_value(sigma)) .or. any(sigma<=0.0_dp)) stop 'Mutation widths must be positive.'
        if (N<2 .or. K<1 .or. nt<1) stop 'Invalid population, generation, or tournament count.'
        if (.not. finite_value(pc) .or. pc<0.0_dp .or. pc>1.0_dp) stop 'Invalid crossover probability.'
        if (.not. finite_value(pm) .or. pm<0.0_dp .or. pm>1.0_dp) stop 'Invalid mutation probability.'

        do i=1,N
            call random_number(child)
            pop(:,i)=lower+(upper-lower)*child
            values(i)=checked_value(f,pop(:,i))
        end do
        evaluations=N
        do g=1,K
            ! Preserve one elite member, including its stored objective value.
            idx=minloc(values)
            next_pop(:,1)=pop(:,idx(1))
            next_values(1)=values(idx(1))
            count=1
            do while (count<N)
                call select_parent(pop,values,nt,p1)
                call select_parent(pop,values,nt,p2)
                children(:,1)=p1
                children(:,2)=p2
                call random_number(u)
                if (u<pc) then
                    ! Uniform crossover exchanges selected coordinates.
                    do j=1,dim
                        call random_number(u)
                        if (u<0.5_dp) then
                            children(j,1)=p2(j)
                            children(j,2)=p1(j)
                        end if
                    end do
                end if
                do i=1,2
                    if (count==N) exit
                    child=children(:,i)
                    call random_number(u)
                    if (u<pm) then
                        ! A selected child receives independent Gaussian kicks
                        ! in all coordinates; widths stay fixed across generations.
                        do j=1,dim
                            child(j)=child(j)+sigma(j)*normal_random()
                        end do
                    end if
                    if (.not. all(finite_value(child))) stop 'Mutation overflow: rescale variables.'
                    child=max(min(child,upper),lower)
                    count=count+1
                    next_pop(:,count)=child
                    next_values(count)=checked_value(f,child)
                    evaluations=evaluations+1
                end do
            end do
            pop=next_pop
            values=next_values
        end do
        idx=minloc(values)
        pbest=pop(:,idx(1))
        fbest=values(idx(1))
    end subroutine genetic_search

    subroutine select_parent(pop,values,nt,parent)
        real(dp), intent(in) :: pop(:,:),values(:)
        integer, intent(in) :: nt
        real(dp), intent(out) :: parent(:)
        integer :: trial,candidate,best,N
        real(dp) :: u
        N=size(values)
        call random_number(u)
        best=1+int(N*u)
        ! Sample tournament competitors with replacement.
        do trial=2,nt
            call random_number(u)
            candidate=1+int(N*u)
            if (values(candidate)<values(best)) best=candidate
        end do
        parent=pop(:,best)
    end subroutine select_parent

    function normal_random() result(z)
        real(dp) :: z,u1,u2
        ! Box-Muller transform; reject zero before taking a logarithm.
        do
            call random_number(u1)
            if (u1>0.0_dp) exit
        end do
        call random_number(u2)
        z=sqrt(-2.0_dp*log(u1))*cos(2.0_dp*acos(-1.0_dp)*u2)
    end function normal_random

    function checked_value(f,p) result(value)
        interface
            function f(p) result(value)
                integer, parameter :: dp=selected_real_kind(15,307)
                real(dp), intent(in) :: p(:)
                real(dp) :: value
            end function f
        end interface
        real(dp), intent(in) :: p(:)
        real(dp) :: value
        value=f(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 genetic_algorithm_methods

module genetic_example
    use genetic_algorithm_methods, only: dp
    implicit none
contains
    function objective(p) result(value)
        real(dp), intent(in) :: p(:)
        real(dp) :: value
        value=sin(3.0_dp*p(1))*cos(3.0_dp*p(2))+(p(1)-1.0_dp)**2+(p(2)+2.0_dp)**2
    end function objective
end module genetic_example

program ch13_genetic_algorithm
    use genetic_algorithm_methods, only: dp, initialize_random, genetic_search
    use genetic_example, only: objective
    implicit none
    ! Editable settings. Mutation probability applies to each child.
    integer, parameter :: dim=2,population_size=100,max_generations=200,tournament_size=3,seed=42
    real(dp), parameter :: lower(dim)=(/-3.0_dp,-3.0_dp/),upper(dim)=(/3.0_dp,3.0_dp/)
    real(dp), parameter :: pc=0.90_dp,pm=0.30_dp,sigma(dim)=0.04_dp*(upper-lower)
    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 genetic_search(objective,lower,upper,population_size,max_generations, &
        pc,pm,sigma,tournament_size,pbest,fbest,evaluations)
    write(*,'(/,a)') 'Genetic algorithm'
    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_genetic_algorithm