"""
    $TYPEDSIGNATURES

Gradient sampling algorthm.
"""
Base.@kwdef struct GradientSampling{Tf} <: NonSmoothOptimizer{Tf}
    m::Int64
    Œ≤::Tf = 1e-4
    Œ≥::Tf = 0.5
    œµ_opt::Tf = 1e-6
    ŒΩ_opt::Tf = 1e-6
    Œ∏_œµ::Tf = 0.1
    Œ∏_ŒΩ::Tf = 0.1
    ls_maxit::Int64 = 70
    nppopt::NearestPointPolytope{Tf}
    # Œ≤::Tf
    # Œ≥::Tf
    # œµ_opt::Tf
    # ŒΩ_opt::Tf
    # Œ∏_œµ::Tf
    # Œ∏_ŒΩ::Tf
    # ls_maxit::Int64
    # nppopt::NearestPointPolytope{Tf}
end
# function GradientSampling(; m, Tf)
#     return GradientSampling{Tf}(
#         m,
#         1e-4,
#         0.5,
#         1e-6,
#         1e-6,
#         0.1,
#         0.1,
#         70,
#         NearestPointPolytope{Tf}(),
#     )
# end

GradientSampling(initial_x::AbstractVector{Tf}; œµ_opt=1e-6, ŒΩ_opt=1e-6) where Tf = GradientSampling{Tf}(
    m = length(initial_x) * 2,
    nppopt = NearestPointPolytope{Tf}(),
    œµ_opt = œµ_opt,
    ŒΩ_opt = ŒΩ_opt,
)

Base.@kwdef mutable struct GradientSamplingState{Tf} <: OptimizerState{Tf}
    x::Vector{Tf}
    ‚àÇg·µ¢s::Matrix{Tf}
    œµ‚Çñ::Tf
    ŒΩ‚Çñ::Tf
    nppstate::NearestPointPolytopeState{Tf}
    k::Int64 = 1
end

function initial_state(gs::GradientSampling{Tf}, initial_x::Vector{Tf}, pb) where {Tf}
    ‚àÇg·µ¢s = zeros(Tf, length(initial_x), gs.m + 1)
    return GradientSamplingState(;
        x = initial_x,
        ‚àÇg·µ¢s,
        œµ‚Çñ = Tf(0.1),
        ŒΩ‚Çñ = Tf(0.1),
        nppstate = QPS.initial_state(‚àÇg·µ¢s),
    )
end


#
### Printing
#
print_header(gs::GradientSampling) = println("**** GradientSampling algorithm\nm = $(gs.m)")

display_logs_header_post(gs::GradientSampling) =
    print("||g·µè||     œµ‚Çñ       ŒΩ‚Çñ        it_ls")
function display_logs_post(os, gs::GradientSampling)
    @printf "%.3e  %.1e  %.1e  %2i" os.additionalinfo.g·µè_norm os.additionalinfo.œµ‚Çñ os.additionalinfo.ŒΩ‚Çñ os.additionalinfo.it_ls
end


#
### GradientSampling method
#
"""
    update_iterate!(state, gs::GradientSampling, pb)

NOTE: each iteration is costly. This can be explored with NonSmoothProblems.to.
On the maxquadBGLS problem
 ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ
                                                            Time                    Allocations
                                                   ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ   ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ
                 Tot / % measured:                      103ms /  93.9%           9.50MiB /  96.1%

 Section                                   ncalls     time    %tot     avg     alloc    %tot      avg
 ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ
 update_iterate!                              100   95.6ms   98.7%   956Œºs   8.92MiB   97.8%  91.4KiB
   GS 2. minimum norm (sub)gradient           100   83.3ms   86.0%   833Œºs   3.56MiB   39.0%  36.4KiB
   GS 1. sampling points, eval gradients      100   7.39ms    7.6%  73.9Œºs   3.85MiB   42.2%  39.4KiB
   GS 4. Update parameters                    100   4.24ms    4.4%  42.4Œºs   1.40MiB   15.3%  14.3KiB
   GS 5. diff check                           100    288Œºs    0.3%  2.88Œºs    117KiB    1.3%  1.17KiB
   GS 3. Termination                          100   41.1Œºs    0.0%   411ns     0.00B    0.0%    0.00B
 build_optimstate                             100   1.19ms    1.2%  11.9Œºs    207KiB    2.2%  2.07KiB
 CV check                                     100   47.0Œºs    0.0%   470ns     0.00B    0.0%    0.00B
 ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ‚îÄ
"""
function update_iterate!(
    state::GradientSamplingState{Tf},
    gs::GradientSampling,
    pb,
) where {Tf}
    iteration_status = iteration_completed
    ‚àÇg·µ¢s = state.‚àÇg·µ¢s

    ## 1. Sample m points in ùîπ(x, œµ‚Çñ)
    samplegradients!(‚àÇg·µ¢s, pb, state.x, state.œµ‚Çñ)

    ## 2. Find minimal norm element of convex hull at gradients of previous points.
    QPS.initialize_state!(state.nppstate, ‚àÇg·µ¢s)
    QPS.nearest_point_polytope!(state.nppstate, gs.nppopt, ‚àÇg·µ¢s; show_trace = false)
    g·µè = state.nppstate.x

    g·µè_norm = norm(g·µè)


    ## 3. termination
    if g·µè_norm ‚â§ gs.ŒΩ_opt && state.œµ‚Çñ ‚â§ gs.œµ_opt
        iteration_status = problem_solved
        @info "termination condition satisfied"
    end

    ## 4. Update parameters
    ŒΩ_next = 0.0
    œµ_next = 0.0
    t‚Çñ = 1.0
    it_ls = 0
    if g·µè_norm ‚â§ state.ŒΩ‚Çñ
        ŒΩ_next = gs.Œ∏_ŒΩ * state.ŒΩ‚Çñ
        œµ_next = gs.Œ∏_œµ * state.œµ‚Çñ
        t‚Çñ = 0.0
        @info "reducing sampling size"
    else
        # This test is not costly, and may help detect difficult cases
        gtd = dot(g·µè, ‚àÇF_elt(pb, state.x))
        if gtd <= 0
            @warn "not descent direction, gtd = $gtd"
        end

        ŒΩ_next = state.ŒΩ‚Çñ
        œµ_next = state.œµ‚Çñ

        f‚Çñ = F(pb, state.x)
        while !(F(pb, state.x - t‚Çñ * g·µè) < f‚Çñ - gs.Œ≤ * t‚Çñ * g·µè_norm^2)
            t‚Çñ *= gs.Œ≥
            it_ls += 1

            if it_ls > gs.ls_maxit
                @warn "GradientSampling(): linesearch exceeded $(gs.ls_maxit) iterations, no suitable steplength found."
                break
            end
            if (norm(g·µè) * t‚Çñ < 10 * eps(Tf))
                @warn("Linesearch reached numerical precision")
                break
            end
        end
    end

    x_next = state.x - t‚Çñ * g·µè
    if !is_differentiable(pb, x_next)
        @error(
            "Gradient sampling: F not differentiable at next point, portion to be implemented."
        )
    end

    state.œµ‚Çñ = œµ_next
    state.ŒΩ‚Çñ = ŒΩ_next
    state.x = x_next
    state.k += 1

    return (;
        œµ‚Çñ = state.œµ‚Çñ,
        ŒΩ‚Çñ = state.ŒΩ‚Çñ,
        it_ls,
        g·µè_norm,
        F = 2 + it_ls,              # orcale calls
        ‚àÇF_elt = gs.m + 1,
        is_differentiable = 1,
    ),
    iteration_status
end

get_minimizer_candidate(state::GradientSamplingState) = state.x


# function find_minimumnormelt_CHP(‚àÇg·µ¢s::Matrix{Tf}) where {Tf}
#     set = CHP.SimplexShadow(‚àÇg·µ¢s)
#     x0 = zeros(Tf, size(‚àÇg·µ¢s, 2))

#     showtermination = true
#     showtrace = true
#     showls = false
#     Œ±, str = CHP.optimize(set, x0; showtermination, showtrace, showls, maxiter = 13)
#     return Œ±
# end

function samplegradients!(‚àÇg·µ¢s, pb, x, œµ‚Çñ)
    n = size(‚àÇg·µ¢s, 1)
    nsamples = size(‚àÇg·µ¢s, 2) - 1

    for i = 1:nsamples
        ‚àÇg·µ¢ = @view ‚àÇg·µ¢s[:, i]
        ‚àÇg·µ¢ .= rand(MvNormal(zeros(n), ScalMat(n, 1.0)))
        ‚àÇg·µ¢ .*= œµ‚Çñ * rand()^(1 / n) / norm(‚àÇg·µ¢)
        ‚àÇg·µ¢ .+= x
        ‚àÇg·µ¢ .= ‚àÇF_elt(pb, ‚àÇg·µ¢)
    end
    ‚àÇg·µ¢s[:, end] .= ‚àÇF_elt(pb, x)
    return
end
