cmdla/10-12/GRS.jl

using Plots
using LinearAlgebra
using Printf

function GRS(f;
        rangeplt::Union{Nothing, Tuple{Real, Real}}=nothing,
        delta::Real=1e-6,
        MaxFeval::Integer=100,
        plt::Union{Plots.Plot, Nothing}=nothing,
        plotatend::Bool=true,
        Plotg::Integer=0
    )::Tuple{Real, String}

    # Plotg
    # 0 = nothing is plotted
    # 1 = the function value / gap are plotted
    # 2 = the function and the test points used are plotted

    resolution = 1000

    Interactive = false # if we pause at every iteration

    # reading and checking input- - - - - - - - - - - - - - - - - - - - - - - -
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    if rangeplt == nothing
        (fStar, _, _, rangeplt) = f(nothing)
    else
        fStar = -Inf
    end

    xm = rangeplt[1] # x_-
    xp = rangeplt[2] # x_+

    if xm > xp
        throw(ArgumentError("rangeplt is empty"))
        return (NaN, "empty")
    end

    if MaxFeval < 2
        throw(ArgumentError("At least two function computations are required"))
        return (NaN, "empty")
    end

    # initializations - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    r = (sqrt( 5 ) - 1)/2

    xmp = xm + (1 - r) * (xp - xm) # x'_-
    xpp = xm + r * (xp - xm)       # x'_+

    fxmp = f(xmp)[1]               # f(x'_-)
    fxpp = f(xpp)[1]               # f(x'_+)

    feval = 2

    if fxmp ≤ fxpp
        fx = fxmp
    else
        fx = fxpp
    end

    status = "optimal"

    if Plotg > 0
        gap = []
    end

    if Plotg == 1 && plt == nothing
        plt = plot(yscale = :log,
                   xlims=(0, 35),
                   ylims=(1e-15, Inf),
                   guidefontsize=16)
    elseif Plotg == 2 && plt == nothing
        plt = plot(legend = false)
    end

    print("Golden ratio search\n")
    if fStar > -Inf
       print("feval\trel gap\t\tx_-\t\tx_+\n")
    else
       print("feval\tfbest\t\tx_-\t\tx_+\n")
    end

    # main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    while xp - xm > delta

        # output statistics - - - - - - - - - - - - - - - - - - - - - - - - - -

        if fStar > -Inf
            gapk = (fx - fStar) / max(abs(fStar), 1)

            if Plotg == 1
                push!(gap, gapk)
            end
        else
            gapk = fx
        end
        @printf("%4d\t%1.4e\t%1.8e\t%1.8e\n", feval, gapk, xm, xp)

        if Plotg == 2
            xbot = xm - ( xp - xm ) / 20;
            xtop = xp + ( xp - xm ) / 20;

            if !isempty(plt.series_list)
                plt.series_list[end][:linealpha] = 0
            end
            xx = range(xbot, xtop, resolution)
            yy = map(v -> v[1], f.(xx))
            xlims!(plt, (xbot, xtop))
            old_ylims = ylims(plt)
            ylims!(plt, (old_ylims[1], max(fxmp, fxpp)))
            plot!(plt, xx, yy)

            # old_xticks = xticks(plt[1])
            new_xticks = ([xbot, xm, xmp, xpp, xp, xtop],
                [@sprintf("%1.1g", xbot), "x-", "x''-", "x''+", "x+", @sprintf("%1.1g", xtop)])
            # keep_indices = findall(x -> all(x .≠ new_xticks[1]), old_xticks[1])
            # merged_xticks = (old_xticks[1][keep_indices] ∪ new_xticks[1],
            #                  old_xticks[2][keep_indices] ∪ new_xticks[2])
            # xticks!(merged_xticks)
            xticks!(plt, new_xticks)
        end

        # stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -
        if feval > MaxFeval
            status = "stopped"
            break
        end

        # main logic- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        if fxmp <= fxpp
            xp = xpp
            xpp = xmp
            xmp = xm + (1 - r) * (xp - xm)

            fxpp = fxmp
            fx = fxmp
            fxmp = f(xmp)[1]
        else
            xm = xmp
            xmp = xpp
            xpp = xm + r * (xp - xm)
            
            fxmp = fxpp
            fx = fxpp
            fxpp = f(xpp)[1]
        end

        feval += 1

        if Plotg ≠ 0 && Interactive
            IJulia.clear_output(wait=true)
            display(plt)
            sleep(0.1)
            readline()
        end

        # iterate - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    end

    # end of main loop- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    # select final answer
    if fxmp ≤ fxpp
        x = xmp
    else
        x = xpp
    end

    if Plotg == 1
        plot!(plt,
            gap,
            linewidth=2)
    end

    if plotatend && Plotg ≠ 0
        display(plt)
    end

    (x, status)
end