! ch13_golden_section.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! One-dimensional minimization by golden-section search.
! Change the objective function and driver settings below.
! The objective must return a finite real value and be unimodal on [a,b].
! Derivatives are not required. Several local minima invalidate the usual
! single-minimum interpretation of the final interval.
!
! Compile (Fortran 95):
!   gfortran -std=f95 -Wall -Wextra -O2 ch13_golden_section.f95 -o ch13_golden_section
! Run: ./ch13_golden_section (Windows: ch13_golden_section.exe)

module golden_section_methods
    implicit none
    private
    integer, parameter :: real64 = selected_real_kind(15,307)
    public :: real64, search_result, golden_section

    type :: search_result
        real(real64) :: xmin, fmin, a, b
        integer :: iterations, evaluations
        logical :: converged
        character(len=90) :: status
    end type search_result


contains

    subroutine golden_section(f, a0, b0, tol_x, max_iter, 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(real64), intent(in) :: a0, b0, tol_x
        integer, intent(in) :: max_iter
        type(search_result), intent(out) :: result
        real(real64) :: a, b, rho, x1, x2, f1, f2

        if (.not. finite_value(a0) .or. .not. finite_value(b0)) &
            stop 'Interval endpoints must be finite.'
        if (b0 <= a0) stop 'Require b > a.'
        if (.not. finite_value(b0-a0)) stop 'Interval width must be finite.'
        if (.not. finite_value(tol_x) .or. tol_x <= 0.0_real64) &
            stop 'Tolerance must be finite and positive.'
        if (max_iter < 0) stop 'Iteration limit must be nonnegative.'

        a = a0
        b = b0
        rho = (sqrt(5.0_real64)-1.0_real64)/2.0_real64
        x1 = b-rho*(b-a)
        x2 = a+rho*(b-a)
        f1 = checked_value(f,x1)
        f2 = checked_value(f,x2)
        result%iterations = 0
        result%evaluations = 2
        result%status = 'Maximum iteration count reached.'

        do while (b-a > tol_x .and. result%iterations < max_iter)
            if (.not. (a < x1 .and. x1 < x2 .and. x2 < b)) then
                result%status = 'Stopped: interval points cannot be resolved in floating point.'
                exit
            end if
            ! Reuse one interior point and its stored function value.
            if (f1 < f2) then
                b = x2
                x2 = x1
                f2 = f1
                x1 = b-rho*(b-a)
                f1 = checked_value(f,x1)
            else
                a = x1
                x1 = x2
                f1 = f2
                x2 = a+rho*(b-a)
                f2 = checked_value(f,x2)
            end if
            result%iterations = result%iterations+1
            result%evaluations = result%evaluations+1
        end do

        result%converged = b-a <= tol_x
        if (result%converged) result%status = 'Interval-width tolerance reached.'
        result%a = a
        result%b = b
        result%xmin = a+(b-a)/2.0_real64
        result%fmin = checked_value(f,result%xmin)
        result%evaluations = result%evaluations+1
        ! Interval width is a numerical localization, not a rigorous error
        ! bound: rounding of function values can impede comparisons.
    end subroutine golden_section

    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(real64), intent(in) :: x
        real(real64) :: 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(real64), intent(in) :: x
        ! False for infinities and NaNs under ordinary IEEE arithmetic.
        finite_value = abs(x) <= huge(x)
    end function finite_value
end module golden_section_methods

! Replace this function to solve a different problem.
module example_objective
    use golden_section_methods, only: real64
    implicit none
contains
    function objective(x) result(value)
        real(real64), intent(in) :: x
        real(real64) :: value
        value = (x-2.0_real64)**2+1.0_real64
    end function objective
end module example_objective

program ch13_golden_section
    use golden_section_methods, only: real64, search_result, golden_section
    use example_objective, only: objective
    implicit none

    ! Editable settings. Example minimum: x = 2, f(x) = 1.
    real(real64), parameter :: a = -1.0_real64, b = 4.0_real64
    real(real64), parameter :: tol_x = 1.0e-6_real64
    integer, parameter :: max_iter = 200
    type(search_result) :: answer

    call golden_section(objective,a,b,tol_x,max_iter,answer)
    write(*,'(/,a)') 'Golden-section search'
    write(*,'(a,es23.15)') 'Estimated minimum: x = ',answer%xmin
    write(*,'(a,es23.15)') 'Objective value:   f = ',answer%fmin
    write(*,'(a,es23.15,a,es23.15,a)') 'Final interval: [',answer%a,', ',answer%b,']'
    write(*,'(a,es14.6)') 'Interval width: ',answer%b-answer%a
    write(*,'(a,i0)') 'Iterations:     ',answer%iterations
    write(*,'(a,i0)') 'Function calls: ',answer%evaluations
    write(*,'(a,a)') 'Status: ',trim(answer%status)
end program ch13_golden_section