using LinearAlgebra
using Printf
using Plots

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

    # Plotg 
    # 1 = the function value / gap are plotted
    # 2 = the function and the model (if used) are plotted
    # all the rest: nothing is 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

    (fxm, f1xm, _) = f(xm)
    if( f1xm > 0 )
        throw(DomainError(rangeplt, "f'(x_-) must be ≤ 0"))
        return (NaN, "empty")
    end

    (fxp, f1xp, _) = f(xp)
    if( f1xp < 0 )
        throw(DomainError(rangeplt, "f'(x_+) must be ≥ 0"))
        return (NaN, "empty")
    end

    if (sfgrd < 0) || (sfgrd ≥ 0.5)
        throw(DomainError(sfgrd, "sfgrd must be in [ 0 , 1/2 )"))
        return (NaN, "empty")
   end

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

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

    feval = 2
    status = "optimal"

    if f1xm ≥ -eps
        return (xm, status)
    end

    if f1xp ≤ eps
        return (xp, status)
    end

    fbest = min(fxp, fxm)
    f1x = min(-f1xp, f1xm)

    if Plotg == 1
        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

    if iszero(sfgrd)
        println("Dichotomic search")
    else
        @printf("Dichotomic search with safeguarded interpolation (%1.4f)\n", sfgrd)
    end

    if fStar > -Inf
        println("feval\trel gap\t\tx_-\t\tx_+\t\tx\t\tf'(x)")
    else
        println("feval\tfbest\t\tx_-\t\tx_+\t\tx\t\tf'(x)")
    end

    # main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    x = 0

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

        # main logic- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        if iszero(sfgrd)
            # compute the new point as the middle of the inteval
            x = (xm + xp) / 2
        else
            # compute the new point by safeguarded quadratic interpolation
            sfxp = xp - (xp - xm) * sfgrd
            sfxm = xm + (xp - xm) * sfgrd
            x = (xm * f1xp - xp * f1xm) / (f1xp - f1xm)
            x = max(sfxm, min(sfxp, x))
        end

        (fx, f1x, _) = f(x) # compute f(x) and f'(x)

        if fx < fbest
            fbest = fx
        end

        feval += 1

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

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

            if Plotg == 1
                push!(gap, gapk)
            end
        else
            gapk = fbest
        end

        @printf("%4d\t%1.4e\t%1.8e\t%1.8e\t%1.8e\t%1.4e\n", feval, gapk, xm, xp, x, f1x)

        if Plotg == 2
            xmp = xm - (xp - xm) / 20
            xpp = xp + (xp - xm) / 20

            for e in plt.series_list
                e[:linealpha] = 0
            end

            xx = range(xmp, xpp, resolution)
            yy = map(v -> v[1], f.(xx))

            xlims!(plt, (xmp, xpp))
            old_ylims = ylims(plt)
            ylims!(plt, (old_ylims[1], max(fxm, fxp)))
            plot!(plt, xx, yy)

            if !iszero(sfgrd)
                a = (f1xp - f1xm) / (2 * (xp - xm))
                b = (xp * f1xm - xm * f1xp) / (xp - xm)
                # a xm^2 + b xm + c == fxm   ==>
                # c == fxm - a xm^2 - b xm
                c = fxm - a * xm^2 - b * xm

                xlims!(plt, (xmp, xpp))
                new_xticks = ([xmp, xm, sfxm, sfxp, xp, xpp],
                    [@sprintf("%1.1g", xmp), "x-", "sx-", "sx+", "x+", @sprintf("%1.1g", xpp)])

                xticks!(plt, new_xticks)
                plot!(plt, xx, @. a*xx^2 + b*xx + c)
            else
                new_xticks = ([xmp, xm, x, xp, xpp],
                    [@sprintf("%1.1g", xmp), "x-", "x", "x+", @sprintf("%1.1g", xpp)])

                xticks!(plt, new_xticks)
            end
        end
        # stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -
        if abs(f1x) ≤ eps
            break
        end

        # restrict the interval based on sign of the derivative in xn - - - - -
        if f1x < 0
            xm = x
            fxm = fx
            f1xm = f1x
        else
            xp = x
            fxp = fx
            f1xp = f1x
        end

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

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

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

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

    (x, status)
end