! ch13_simulated_annealing.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! Bounded one-dimensional minimization by simulated annealing.
! Edit the objective and driver settings below. No derivatives are required.
! Repeat with different seeds and cooling schedules to assess reliability.
! The best point found is not proof of a global minimum.
!
! Compile (Fortran 95):
! gfortran -std=f95 -Wall -Wextra -O2 ch13_simulated_annealing.f95 -o ch13_simulated_annealing
! Run: ./ch13_simulated_annealing (Windows: ch13_simulated_annealing.exe)

module annealing_methods
    implicit none
    private
    integer, parameter :: dp=selected_real_kind(15,307)
    public :: dp, annealing_result, simulated_annealing, initialize_random
    type :: annealing_result
        real(dp) :: xbest, fbest, temperature, acceptance_fraction
        integer :: levels, trials, evaluations, outside
        character(len=80) :: status
    end type annealing_result
contains
    subroutine initialize_random(seed)
        integer, intent(in) :: seed
        integer, allocatable :: seeds(:)
        integer :: n, i
        call random_seed(size=n)
        allocate(seeds(n))
        ! Bounded arithmetic avoids integer overflow when forming seed values.
        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)
        ! Reproducible on the same compiler/runtime; sequences are compiler-dependent.
    end subroutine initialize_random

    subroutine simulated_annealing(f,a,b,x0,T0,T_min,cooling,moves,step_size,max_levels,result)
        interface
            function f(x) result(value)
                integer, parameter :: dp=selected_real_kind(15,307)
                real(dp), intent(in) :: x
                real(dp) :: value
            end function f
        end interface
        real(dp), intent(in) :: a,b,x0,T0,T_min,cooling,step_size
        integer, intent(in) :: moves,max_levels
        type(annealing_result), intent(out) :: result
        real(dp) :: x,fx,x_trial,f_trial,df,T,u
        integer :: j,accepted
        logical :: accept

        if (.not. finite_value(a) .or. .not. finite_value(b)) stop 'Bounds must be finite.'
        if (b<=a) stop 'Require upper bound > lower bound.'
        if (.not. finite_value(x0) .or. x0<a .or. x0>b) stop 'Starting point outside bounds.'
        if (.not. finite_value(T0) .or. T0<=0.0_dp) stop 'Invalid initial temperature.'
        if (.not. finite_value(T_min) .or. T_min<=0.0_dp .or. T_min>=T0) stop 'Invalid final temperature.'
        if (.not. finite_value(cooling) .or. cooling<=0.0_dp .or. cooling>=1.0_dp) stop 'Invalid cooling factor.'
        if (.not. finite_value(step_size) .or. step_size<=0.0_dp) stop 'Invalid step size.'
        if (moves<=0 .or. max_levels<=0) stop 'Search counts must be positive.'

        x=x0
        fx=checked_value(f,x)
        result%xbest=x
        result%fbest=fx
        result%levels=0
        result%trials=0
        result%outside=0
        result%evaluations=1
        accepted=0
        T=T0
        do while (T>T_min .and. result%levels<max_levels)
            do j=1,moves
                result%trials=result%trials+1
                call random_number(u)
                x_trial=x+step_size*(2.0_dp*u-1.0_dp)
                ! Reject outside proposals; do not clip them to the bounds.
                if (.not. finite_value(x_trial) .or. x_trial<a .or. x_trial>b) then
                    result%outside=result%outside+1
                    cycle
                end if
                f_trial=checked_value(f,x_trial)
                result%evaluations=result%evaluations+1
                df=f_trial-fx
                accept=df<=0.0_dp
                if (.not. accept) then
                    call random_number(u)
                    ! Avoid log(0); uphill acceptance is exp(-df/T).
                    if (u==0.0_dp) then
                        accept=.true.
                    else
                        accept=log(u)<-df/T
                    end if
                end if
                if (accept) then
                    x=x_trial
                    fx=f_trial
                    accepted=accepted+1
                    if (fx<result%fbest) then
                        result%xbest=x
                        result%fbest=fx
                    end if
                end if
            end do
            result%levels=result%levels+1
            T=cooling*T
        end do
        result%temperature=T
        result%acceptance_fraction=real(accepted,dp)/real(result%trials,dp)
        result%status='Maximum temperature-level count reached.'
        if (T<=T_min) result%status='Cooling schedule completed.'
    end subroutine simulated_annealing

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

module annealing_example
    use annealing_methods, only: dp
    implicit none
contains
    function objective(x) result(value)
        real(dp), intent(in) :: x
        real(dp) :: value
        value=x**4-20.0_dp*x*x+24.0_dp*x+5.0_dp
    end function objective
end module annealing_example

program ch13_simulated_annealing
    use annealing_methods, only: dp, annealing_result, simulated_annealing, initialize_random
    use annealing_example, only: objective
    implicit none
    ! Editable problem and search settings.
    real(dp), parameter :: a=-5.0_dp,b=4.0_dp,x0=2.0_dp
    real(dp), parameter :: T0=100.0_dp,T_min=1.0e-3_dp,cooling=0.95_dp,step_size=0.5_dp
    integer, parameter :: moves_per_temp=100,max_levels=1000,seed=1
    ! integer :: clock_seed  ! Uncomment when using clock-based seeding.
    type(annealing_result) :: answer

    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 simulated_annealing(objective,a,b,x0,T0,T_min,cooling, &
        moves_per_temp,step_size,max_levels,answer)
    write(*,'(/,a)') 'Simulated annealing'
    write(*,'(a,es23.15)') 'Best point found:     ',answer%xbest
    write(*,'(a,es23.15)') 'Best objective value: ',answer%fbest
    write(*,'(a,i0)') 'Temperature levels:   ',answer%levels
    write(*,'(a,i0)') 'Trial moves:          ',answer%trials
    write(*,'(a,i0)') 'Function calls:       ',answer%evaluations
    write(*,'(a,f8.2,a)') 'Accepted moves:       ',100.0_dp*answer%acceptance_fraction,' %'
    write(*,'(a,i0)') 'Out-of-bounds trials: ',answer%outside
    write(*,'(a,es14.6)') 'Temperature at stop:  ',answer%temperature
    write(*,'(a,a)') 'Status: ',trim(answer%status)
    write(*,'(a)') 'Repeat with other seeds; global optimality is not guaranteed.'
end program ch13_simulated_annealing