text:
#Algorithms
abstract type Bandit end
abstract type BanditResult end
@with_kw mutable struct RandomBandit <: Bandit
    actiondistr::ActionDistr  = UniformActions()
    seed::Int               = 0
    n_iters::Int            = 100
    outs::Set{Symbol}       = Set{Symbol}()
end
metadata(x) = []
Base.string(b::RandomBandit) = "Random"
Base.string(b::Type{RandomBandit}) = "Random"
metadata(b::RandomBandit) = [:algorithm=>string(b), :seed=>b.seed, :n_iters=>b.n_iters]
struct RandomBanditResult <: BanditResult
    actions::Vector
    qs::Vector{Float64}
    ns::Vector{Int}
    cum_regret::Vector{Float64}
    plts::Vector{Plots.Plot}
    simple_regret::Vector{Float64}
end
RandomBanditResult() = RandomBanditResult([], Float64[], Int[], Float64[], Plots.Plot[], Float64[])
function POMDPs.solve(b::RandomBandit, G::ObjectiveFunc, rng::AbstractRNG=MersenneTwister(b.seed))
    actions = []
    qs = Float64[]
    cum_regret = 0
    g_max = maximum(G)
    result = RandomBanditResult()
    for n = 1:b.n_iters
        x = rand(rng, b.actiondistr)
        push!(actions, x)
        r = G(x, rng)
        push!(qs, r)

        if :simple_regret in b.outs
            i = indmax(qs)
            g_best = truth(G, actions[i])
            push!(result.simple_regret, g_max-g_best)
        end
        if :cum_regret in b.outs
            cum_regret += (g_max - r)
            push!(result.cum_regret, cum_regret)
        end
        if :plot in b.outs
            p = plot(G, b, actions, qs)
            push!(result.plts, p)
        end
    end
    append!(result.actions, actions)
    append!(result.qs, qs)
    append!(result.ns, ones(length(actions)))
    result
end
@recipe function plot(G::Objective1D, b::RandomBandit, actions, qs)
    title := string(b)
    ylim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        actions, qs
    end
end
@recipe function plot(G::Objective2D, b::RandomBandit, actions, qs)
    title := string(b)
    zlim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        zcolor := qs
        xs = map(x->x[1], actions)
        ys = map(x->x[2], actions)
        xs, ys
    end
end
function Plots.animate(result::BanditResult, filename="./result.gif"; fps=5, every=2)
    animate(result.plts, filename, fps=fps, every=every)
end

@with_kw mutable struct PWUCB <: Bandit
    actiondistr::ActionDistr  = UniformActions()
    seed::Int               = 0
    n_iters::Int            = 100
    k::Float64              = 0.5
    α::Float64              = 0.85
    ec::Float64             = 0.5
    q0::Float64             = 0
    n0::Float64             = 0
    outs::Set{Symbol}       = Set{Symbol}()
end
metadata(b::PWUCB) = [:algorithm=>"PW", :seed=>b.seed, :n_iters=>b.n_iters, :k=>b.k, :α=>b.α, :ec=>b.ec,
                       :q0=>b.q0, :n0=>b.n0]
Base.string(b::PWUCB) = "PW"
Base.string(b::Type{PWUCB}) = "PW"
struct PWUCBResult <: BanditResult
    actions::Vector
    qs::Vector{Float64}
    ns::Vector{Int}
    cum_regret::Vector{Float64}
    plts::Vector{Plots.Plot}
    simple_regret::Vector{Float64}
end
PWUCBResult() = PWUCBResult([], Float64[], Int[], Float64[], Plots.Plot[], Float64[])
function POMDPs.solve(b::PWUCB, G::ObjectiveFunc, rng::AbstractRNG=MersenneTwister(b.seed))
    actions = []
    qs = Float64[]
    ns = Float64[]
    result = PWUCBResult()
    cum_regret = 0
    g_max = maximum(G)
    for n = 1:b.n_iters
        if length(actions) < b.k*n^b.α
            a = rand(rng, b.actiondistr)
            push!(actions, a)
            push!(qs, b.q0)
            push!(ns, b.n0)
        end
        ebs = b.ec .* [sqrt(log(n)/na) for na in ns]
        ebs[isnan.(ebs)] = 0
        ucb_vals = [q+eb for (q,eb) in zip(qs,ebs)]

        i = indmax(ucb_vals)
        x = actions[i]
        r = G(x, rng)
        ns[i] += 1
        qs[i] += (r - qs[i]) / ns[i]

        if :simple_regret in b.outs
            i = indmax(qs)
            g_best = truth(G, actions[i])
            push!(result.simple_regret, g_max-g_best)
        end
        if :cum_regret in b.outs
            cum_regret += (g_max - r)
            push!(result.cum_regret, cum_regret)
        end
        if :plot in b.outs
            p = plot(G, b, actions, qs, ebs)
            push!(result.plts, p)
        end
    end
    append!(result.actions, actions)
    append!(result.qs, qs)
    append!(result.ns, ones(length(actions)))
    result
end
@recipe function plot(G::Objective1D, b::PWUCB, actions, qs, errbars)
    title := string(b)
    ylim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        title := string(b)
        err := errbars
        actions, qs
    end
end
@recipe function plot(G::Objective2D, b::PWUCB, actions, qs, errbars)
    title := string(b)
    zlim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        #err := errbars  #Figure out how to plot this...
        zcolor := qs
        xs = map(x->x[1], actions)
        ys = map(x->x[2], actions)
        xs, ys
    end
end

@with_kw mutable struct SBUCB <: Bandit
    actiondistr::ActionDistr  = UniformActions()
    seed::Int               = 0
    n_iters::Int            = 100
    k::Float64              = 0.5
    α::Float64              = 0.5
    ec::Float64             = 0.5
    q0::Float64             = 0
    n0::Float64             = 0
    outs::Set{Symbol}       = Set{Symbol}()
end
Base.string(b::SBUCB) = "SB"
Base.string(b::Type{SBUCB}) = "SB"
metadata(b::SBUCB) = [:algorithm=>string(b), :seed=>b.seed, :n_iters=>b.n_iters, :k=>b.k, :α=>b.α, :ec=>b.ec,
                       :q0=>b.q0, :n0=>b.n0]
struct SBUCBResult <: BanditResult
    actions::Vector
    qs::Vector{Float64}
    ns::Vector{Int}
    cum_regret::Vector{Float64}
    plts::Vector{Plots.Plot}
    simple_regret::Vector{Float64}
end
SBUCBResult() = SBUCBResult([], Float64[], Int[], Float64[], Plots.Plot[], Float64[])
function POMDPs.solve(b::SBUCB, G::ObjectiveFunc, rng::AbstractRNG=MersenneTwister(b.seed))
    actions = []
    qs = Float64[]
    ns = Float64[]
    radii = Float64[]
    result = SBUCBResult()
    cum_regret = 0
    g_max = maximum(G)
    for n = 1:b.n_iters
        a = rand(rng, b.actiondistr)
        if isempty(actions)
            d_nn = Inf
        else
            d_nn,i_nn = findmin(dist(a,aa) for aa in actions)
        end
        if d_nn > b.k/n^b.α
            push!(actions, a)
            push!(qs, b.q0)
            push!(ns, b.n0)
            push!(radii, 0)
        end #else, discard
        ebs = b.ec .* [sqrt(log(n)/na) for na in ns]
        ebs[isnan.(ebs)] = 0
        ucb_vals = [q+eb for (q,eb) in zip(qs,ebs)]

        i = indmax(ucb_vals)
        x = actions[i]
        r = G(x, rng)
        ns[i] += 1
        qs[i] += (r - qs[i]) / ns[i]

        if :cum_regret in b.outs
            cum_regret += (g_max - r)
            push!(result.cum_regret, cum_regret)
        end
        if :plot in b.outs
            fill!(radii, b.k/n^b.α)
            p = plot(G, b, actions, qs, ebs, radii)
            push!(result.plts, p)
        end
        if :simple_regret in b.outs
            i = indmax(qs)
            g_best = truth(G, actions[i])
            push!(result.simple_regret, g_max-g_best)
        end
    end
    append!(result.actions, actions)
    append!(result.qs, qs)
    append!(result.ns, ones(length(actions)))
    result
end
dist(x, y) = norm(x-y)
@recipe function plot(G::Objective1D, b::SBUCB, actions, qs, errbars, radii)
    title := string(b)
    ylim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        err := errbars
        actions, qs
    end
    @series begin
        fillalpha := 0.25
        [circle(actions[i],qs[i],radii[i]) for i=1:length(actions)]
    end
end
@recipe function plot(G::Objective2D, b::SBUCB, actions, qs, errbars, radii)
    title := string(b)
    zlim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        #err := errbars  #TODO: how to plot this?
        zcolor := qs
        xs = map(x->x[1], actions)
        ys = map(x->x[2], actions)
        xs, ys
    end
    @series begin
        fillalpha := 0.25
        [circle(actions[i]...,radii[i]) for i=1:length(actions)]
    end
end

#See Srinivas2010 and Dorrard2009
@with_kw mutable struct GPUCBGrid <: Bandit
    actiondistr::ActionDistr = UniformActions()
    grid                     = GridActions()
    seed::Int               = 0
    n_iters::Int            = 100
    n_sig::Float64          = 2.0  #number of stdev's to add to mean (optimism)
    mean_init::Float64      = 0.0
    log_length_scale        = -2.6
    log_signal_sigma        = 0.0
    log_obs_noise           = -1.0
    outs::Set{Symbol}       = Set{Symbol}()
end
Base.string(b::GPUCBGrid) = "GPUCB-Grid"
Base.string(b::Type{GPUCBGrid}) = "GPUCB-Grid"
metadata(b::GPUCBGrid) = [:algorithm=>string(b), :seed=>b.seed, :n_iters=>b.n_iters, :n_sig=>b.n_sig,
                          :mean_init=>b.mean_init, :log_length_scale=>b.log_length_scale,
                          :log_signal_sigma=>b.log_signal_sigma,
                          :log_obs_noise=>b.log_obs_noise]
struct GPUCBGridResult <: BanditResult
    actions::Vector
    qs::Vector{Float64}
    ns::Vector{Int}
    cum_regret::Vector{Float64}
    plts::Vector{Plots.Plot}
    simple_regret::Vector{Float64}
end
GPUCBGridResult() = GPUCBGridResult([], Float64[], Int[], Float64[], Plots.Plot[], Float64[])
function POMDPs.solve(b::GPUCBGrid, G::ObjectiveFunc, rng::AbstractRNG=MersenneTwister(b.seed))
    actions = initial_actions(b.actiondistr)
    qs = Float64[]
    result = GPUCBGridResult()
    cum_regret = 0
    g_max = maximum(G)
    gridpoints = grid_points(b.grid)
    mconst = MeanConst(b.mean_init)
    kern = SE(b.log_length_scale,b.log_signal_sigma)
    gp = GP(vec2mat(actions), qs, mconst, kern, b.log_obs_noise)
    if :plot in b.outs #plot needs these vars initialized
        m, Σ = predict_f(gp, gridpoints)
        ucb = b.n_sig*sqrt.(Σ)
        ucbmax,i = findmax(m + ucb)
    end
    for n = 1:b.n_iters
        if n == 1
            a = rand(rng, b.actiondistr)
            push!(actions, a)
            imax = 1
        else
            GaussianProcesses.fit!(gp, vec2mat(actions), qs)
            #optimize!(gp; method=Optim.BFGS())
            m, Σ = predict_f(gp, gridpoints)
            ucb = b.n_sig*sqrt.(Σ)
            ucbmax,imax = findmax(m + ucb)
            push!(actions, getindex_mat(gridpoints,imax))
        end
        x = actions[end]
        r = G(x, rng)
        push!(qs, r)

        if :cum_regret in b.outs
            cum_regret += (g_max - r)
            push!(result.cum_regret, cum_regret)
        end
        if :plot in b.outs
            p = plot(G, b, gridpoints, m, actions, qs, ucb, ucbmax)
            push!(result.plts, p)
        end
        if :simple_regret in b.outs
            if n == 1
                xbest = x
            else
                i = indmax(m)
                xbest = getindex_mat(gridpoints,i)
            end
            g_best = truth(G, xbest)
            push!(result.simple_regret, g_max-g_best)
        end
    end
    append!(result.actions, actions)
    append!(result.qs, qs)
    append!(result.ns, ones(length(actions)))
    result
end
initial_actions(A::DistributionActions) = Float64[]
initial_actions(A::Union{UniformActions2D,DistributionActions2D}) = Vector{Float64}[]
vec2mat(xs::Vector{Float64}) = xs[:,1:1]'
vec2mat(xs::Vector{Vector{Float64}}) = isempty(xs) ? Array{Float64,2}(2,0) : hcat(xs...)
getindex_mat(xs::Vector{Float64},i) = xs[i]
getindex_mat(xs::Vector{Vector{Float64}},i) = xs[i]
getindex_mat(xs::Array{Float64,2},i) = xs[:,i]
@recipe function plot(G::Objective1D, b::GPUCBGrid, gridpoints, m, actions, qs, ucb, ucbmax)
    title := string(b)
    ylim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        label := "observations"
        actions[1:end-1], qs[1:end-1]
    end
    @series begin
        seriestype := :path
        linestyle := :dash
        ribbon := ucb
        label := "predicted mean"
        gridpoints, m
    end
    @series begin
        seriestype := :scatter
        markershape := :star4
        label := "ucb max"
        actions[end:end], [ucbmax]
    end
end
@recipe function plot(G::Objective2D, b::GPUCBGrid, gridpoints, m, actions, qs, ucb, ucbmax)
    title := string(b)
    zlim := (G.g_min - 1.0, G.g_max + 1.0)
    layout := @layout grid(1,3)
    @series begin
        subplot := 1
        G
    end
    @series begin
        subplot := 1
        seriestype := :scatter
        label := "observations"
        xs = map(x->x[1], actions[1:end-1])
        ys = map(x->x[2], actions[1:end-1])
        xs, ys
    end
    @series begin
        subplot := 2
        seriestype := :heatmap
        title := "predicted mean"
        legend := true
        b.grid.xs, b.grid.ys, reshape(m, (b.grid.n,b.grid.n))'
    end
    @series begin
        subplot := 3
        seriestype := :heatmap
        title := "ucb"
        legend := true
        b.grid.xs, b.grid.ys, reshape(ucb, (b.grid.n,b.grid.n))'
    end
end

@with_kw mutable struct GPUCB <: Bandit
    actiondistr::ActionDistr  = UniformActions()
    grid                     = GridActions()  #for plotting only
    seed::Int               = 0
    n_iters::Int            = 100
    k::Int                  = 10
    n_sig::Float64          = 2.0  #number of stdev's to add to mean (optimism)
    mean_init::Float64      = 0.0
    log_length_scale        = -2.6
    log_signal_sigma        = 0.0
    log_obs_noise           = -1.0
    outs::Set{Symbol}       = Set{Symbol}()
end
Base.string(b::GPUCB) = "GPUCB"
Base.string(b::Type{GPUCB}) = "GPUCB"
metadata(b::GPUCB) = [:algorithm=>string(b), :seed=>b.seed, :n_iters=>b.n_iters, :k=>b.k, :n_sig=>b.n_sig,
                      :mean_init=>b.mean_init, :log_length_scale=>b.log_length_scale,
                      :log_signal_sigma=>b.log_signal_sigma,
                      :log_obs_noise=>b.log_obs_noise]
struct GPUCBResult <: BanditResult
    actions::Vector
    qs::Vector{Float64}
    ns::Vector{Int}
    cum_regret::Vector{Float64}
    plts::Vector{Plots.Plot}
    simple_regret::Vector{Float64}
end
GPUCBResult() = GPUCBResult([], Float64[], Int[], Float64[], Plots.Plot[], Float64[])
function POMDPs.solve(b::GPUCB, G::ObjectiveFunc, rng::AbstractRNG=MersenneTwister(b.seed))
    actions = initial_actions(b.actiondistr)
    qs = Float64[]
    result = GPUCBResult()
    cum_regret = 0
    g_max = maximum(G)
    gridpoints = grid_points(b.grid) #for plotitng only
    mconst = MeanConst(b.mean_init)
    kern = SE(b.log_length_scale,b.log_signal_sigma)
    gp = GP(vec2mat(actions), qs, mconst, kern, b.log_obs_noise)
    for n = 1:b.n_iters
        if n == 1
            a = rand(rng, b.actiondistr)
            push!(actions, a)
            imax = 1
        else
            GaussianProcesses.fit!(gp, vec2mat(actions), qs)  #fit over existing actions
            #optimize!(gp; method=Optim.BFGS())
            xs = vcat(actions, [rand(rng, b.actiondistr) for i=1:b.k])
            m, Σ = predict_f(gp, vec2mat(xs))  #over actions and k random samples from proposal distribution
            ucb = b.n_sig*sqrt.(Σ)
            ucbmax,imax = findmax(m + ucb)
            push!(actions, xs[imax])
        end
        x = actions[end]
        r = G(x, rng)
        push!(qs, r)

        if :cum_regret in b.outs
            cum_regret += (g_max - r)
            push!(result.cum_regret, cum_regret)
        end
        if :plot in b.outs
            if n != 1
                m´, Σ´ = predict_f(gp, gridpoints) #over plotting grid
                ucb´ = b.n_sig*sqrt.(Σ´)
                ucbmax´ = maximum(m´ + ucb´)
                p = plot(G, b, gridpoints, m´, actions, qs, ucb´, ucbmax´)
                push!(result.plts, p)
            end
        end
        if :simple_regret in b.outs
            if n == 1
                xbest = x
            else
                i = indmax(m)
                xbest = xs[i]
            end
            g_best = truth(G, xbest)
            push!(result.simple_regret, g_max-g_best)
        end
    end
    append!(result.actions, actions)
    append!(result.qs, qs)
    append!(result.ns, ones(length(actions)))
    result
end
@recipe function plot(G::Objective1D, b::GPUCB, xs, m, actions, qs, ucb, ucbmax)
    title := string(b)
    ylim := (G.g_min - 1.0, G.g_max + 1.0)
    @series begin
        G
    end
    @series begin
        seriestype := :scatter
        label := "observations"
        actions[1:end-1], qs[1:end-1]
    end
    @series begin
        seriestype := :path
        linestyle := :dash
        ribbon := ucb
        label := "predicted mean"
        xs, m
    end
    @series begin
        seriestype := :scatter
        markershape := :star4
        label := "ucb max"
        actions[end:end], [ucbmax]
    end
end
@recipe function plot(G::Objective2D, b::GPUCB, gridpoints, m, actions, qs, ucb, ucbmax)
    title := string(b)
    zlim := (G.g_min - 1.0, G.g_max + 1.0)
    layout := @layout grid(1,3)
    @series begin
        subplot := 1
        G
    end
    @series begin
        subplot := 1
        seriestype := :scatter
        label := "observations"
        xs = map(x->x[1], actions[1:end-1])
        ys = map(x->x[2], actions[1:end-1])
        xs, ys
    end
    @series begin
        subplot := 2
        seriestype := :heatmap
        title := "predicted mean"
        legend := true
        b.grid.xs, b.grid.ys, reshape(m, (b.grid.n,b.grid.n))'
    end
    @series begin
        subplot := 3
        title := "ucb"
        legend := true
        b.grid.xs, b.grid.ys, reshape(ucb, (b.grid.n,b.grid.n))'
    end
end


text:
using Test
using StaticArrays
using ClimateMachine.Atmos:
    AtmosProblem,
    NoReferenceState,
    AtmosPhysics,
    AtmosModel,
    DryModel,
    ConstantDynamicViscosity,
    AtmosLESConfigType,
    HBFVReconstruction
import ClimateMachine.DGMethods.FVReconstructions: FVConstant, FVLinear, width
using ClimateMachine.Orientations
import StaticArrays: SUnitRange
import ClimateMachine.BalanceLaws:
    Primitive, Prognostic, vars_state, number_states
using ClimateMachine.VariableTemplates: Vars
using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# lin_func, quad_func, third_func, fourth_func
#
# ```
# num_state_primitive::Int64
# pointwise values::Array{FT,   num_state_primitive by length(ξ)}
# integration values::Array{FT, num_state_primitive by length(ξ)}
# ```

function lin_func(ξ)
    return 1, [(2 * ξ .+ 1)';], [(ξ .^ 2 .+ ξ)';]
end

function quad_func(ξ)
    return 2,
    [(3 * ξ .^ 2 .+ 1)'; (2 * ξ .+ 1)'],
    [(ξ .^ 3 .+ ξ)'; (ξ .^ 2 .+ ξ)']
end

function third_func(ξ)
    return 1, [(4 * ξ .^ 3 .+ 1)';], [(ξ .^ 4 .+ ξ)';]
end

function fourth_func(ξ)
    return 2,
    [(5 * ξ .^ 4 .+ 1)'; (3 * ξ .^ 2 .+ 1)'],
    [(ξ .^ 5 .+ ξ)'; (ξ .^ 3 .+ ξ)']
end

@testset "Hydrostatic balanced linear reconstruction test" begin

    function initialcondition!(problem, bl, state, aux, coords, t, args...) end

    fv_recon! = FVLinear()
    stencil_width = width(fv_recon!)
    stencil_center = stencil_width + 1
    stencil_diameter = 2stencil_width + 1
    @test stencil_width == 1
    func = lin_func

    for FT in (Float64,)
        physics = AtmosPhysics{FT}(
            param_set;
            ref_state = NoReferenceState(),
            turbulence = ConstantDynamicViscosity(FT(0)),
            moisture = DryModel(),
        )
        model = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            problem = AtmosProblem(;
                physics = physics,
                init_state_prognostic = initialcondition!,
            ),
            orientation = FlatOrientation(),
        )

        vars_prim = Vars{vars_state(model, Primitive(), FT)}

        hb_recon! = HBFVReconstruction(model, fv_recon!)
        _grav = FT(grav(param_set))

        @test width(hb_recon!) == 1

        num_state_prognostic = number_states(model, Prognostic())
        num_state_primitive = number_states(model, Primitive())
        local_state_face_primitive = ntuple(Val(2)) do _
            MArray{Tuple{num_state_primitive}, FT}(undef)
        end

        # interior point reconstruction test
        grid = FT[0; 1; 3; 6]
        grid_c = (grid[2:end] + grid[1:(end - 1)]) / 2
        local_cell_weights =
            MArray{Tuple{stencil_diameter}, FT}(grid[2:end] - grid[1:(end - 1)])

        # linear profile for all variables expect pressure
        local_state_primitive = SVector(ntuple(Val(stencil_diameter)) do _
            MArray{Tuple{num_state_primitive}, FT}(undef)
        end...)

        # values at the cell centers       1     2      3
        _, uc, _ = func(grid_c)
        # values at the cell faces     0.5   1.5*   2.5*     3.5
        _, uf, _ = func(grid)
        for i_d in 1:stencil_diameter
            for i_p in 1:num_state_prognostic
                local_state_primitive[i_d][i_p] = uc[i_d]
            end
        end

        # pressure profile is updated to satisfy the discrete hydrostatic balance
        p_surf = FT(100) # at the bottom wall
        p_ref = p_surf
        for i_d in 1:stencil_diameter
            p_ref -=
                vars_prim(local_state_primitive[i_d]).ρ *
                _grav *
                local_cell_weights[i_d] / 2
            vars_prim(local_state_primitive[i_d]).p = p_ref
            p_ref -=
                vars_prim(local_state_primitive[i_d]).ρ *
                _grav *
                local_cell_weights[i_d] / 2
        end

        local_state_primitive_hb = copy(local_state_primitive)

        # interior point reconstruction
        hb_recon!(
            local_state_face_primitive[1],
            local_state_face_primitive[2],
            local_state_primitive,
            local_cell_weights,
        )
        # bottom face
        @test vars_prim(local_state_face_primitive[1]).ρ ≈ uf[2]
        @test vars_prim(local_state_face_primitive[1]).p ≈
              p_surf -
              vars_prim(local_state_primitive[1]).ρ *
              _grav *
              local_cell_weights[1]
        # top face
        @test vars_prim(local_state_face_primitive[2]).ρ ≈ uf[3]
        @test vars_prim(local_state_face_primitive[2]).p ≈
              p_surf -
              vars_prim(local_state_primitive[1]).ρ *
              _grav *
              local_cell_weights[1] -
              vars_prim(local_state_primitive[2]).ρ *
              _grav *
              local_cell_weights[2]

        # make sure the
        for i_d in 1:stencil_diameter
            @test all(
                local_state_primitive[i_d] ≈ local_state_primitive_hb[i_d],
            )
        end

        @info "Start boundary test"

        # boundary point reconstruction
        rng = SUnitRange(stencil_center, stencil_center)

        hb_recon!(
            local_state_face_primitive[1],
            local_state_face_primitive[2],
            local_state_primitive[rng],
            local_cell_weights[rng],
        )

        # bottom face
        @test vars_prim(local_state_face_primitive[1]).ρ ≈ uc[stencil_center]
        @test vars_prim(local_state_face_primitive[1]).p ≈
              vars_prim(local_state_primitive[stencil_center]).p +
              uc[stencil_center] * _grav * local_cell_weights[stencil_center] /
              2

        # top face
        @test vars_prim(local_state_face_primitive[2]).ρ ≈ uc[stencil_center]
        @test vars_prim(local_state_face_primitive[2]).p ≈
              vars_prim(local_state_primitive[stencil_center]).p -
              uc[stencil_center] * _grav * local_cell_weights[stencil_center] /
              2

    end

end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Maps a string key to a path within a volume.

    IoK8sApiCoreV1KeyToPath(;
        key=nothing,
        mode=nothing,
        path=nothing,
    )

    - key::String : The key to project.
    - mode::Int32 : Optional: mode bits to use on this file, must be a value between 0 and 0777. If not specified, the volume defaultMode will be used. This might be in conflict with other options that affect the file mode, like fsGroup, and the result can be other mode bits set.
    - path::String : The relative path of the file to map the key to. May not be an absolute path. May not contain the path element '..'. May not start with the string '..'.
"""
mutable struct IoK8sApiCoreV1KeyToPath <: SwaggerModel
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    mode::Any # spec type: Union{ Nothing, Int32 } # spec name: mode
    path::Any # spec type: Union{ Nothing, String } # spec name: path

    function IoK8sApiCoreV1KeyToPath(;key=nothing, mode=nothing, path=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        o
    end
end # type IoK8sApiCoreV1KeyToPath

const _property_map_IoK8sApiCoreV1KeyToPath = Dict{Symbol,Symbol}(Symbol("key")=>Symbol("key"), Symbol("mode")=>Symbol("mode"), Symbol("path")=>Symbol("path"))
const _property_types_IoK8sApiCoreV1KeyToPath = Dict{Symbol,String}(Symbol("key")=>"String", Symbol("mode")=>"Int32", Symbol("path")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1KeyToPath }) = collect(keys(_property_map_IoK8sApiCoreV1KeyToPath))
Swagger.property_type(::Type{ IoK8sApiCoreV1KeyToPath }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1KeyToPath[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1KeyToPath }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1KeyToPath[property_name]

function check_required(o::IoK8sApiCoreV1KeyToPath)
    (getproperty(o, Symbol("key")) === nothing) && (return false)
    (getproperty(o, Symbol("path")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1KeyToPath }, name::Symbol, val)
end


text:
using Images, ImageView
img = load("sample-images/cats-3061372_640.jpg")
img_channel_view = channelview(img)
imshow(img_channel_view)

# update all colors with a value over 0.7 to be 0.9:
using Images, ImageView
img = load("sample-images/cats-3061372_640.jpg")
img_channel_view = channelview(img)
img_channel_view[img_channel_view .> 0.7] = 0.9
imshow(img)


text:
function genTikhonov(A,mr,p,Aufloesung)
lrows,lcols = size(A);
L = zeros(lcols,lcols);
for r = 1:lcols-1
    L[r,r] = -2;
    L[r,r+1] = 1;
    if r <= (Aufloesung^2-Aufloesung)
        L[r,r+Aufloesung] = 1;
    end
end
L[lcols,lcols] = -2;

ATA = transpose(A)*A;
LTL = transpose(L)*L;
GenTikVek = (ATA+p*LTL)\transpose(A)*mr;

return GenTikVek;
end


text:
#############
#ad-hoc way of getting the analytical hessian:
function hessian2(at)
    hess = JuLIP.hessian_pos(at)
    return JuLIP._pos_to_dof(hess,at)
end

function hessian2(at,dfs)
    set_dofs!(at,dfs)
    return hessian2(at)
end

#function to plot and highlight different regions
function plot_at(cluster)
    N1 = cluster["N1"]
    for x in cluster.X[1:N1]
        plot(x[1],x[2],"g.")
    end

    for x in cluster.X[I_R2]
         plot(x[1],x[2],"b.")
    end

    for x in cluster.X[I_R3]
        plot(x[1],x[2],"r.")
    end

    for x in cluster.X[I_R4]
        plot(x[1],x[2],"k.")
    end
end

#function to plot a configuration specified by the triplet (Ubar,K,alpha)
#ncf is the NCFlex struct
function plot_at_ncf(ncf,Ubar,K,α)
    at = ncf.at

    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
    Uhat = v2d(U_CLE(at,α))

    set_dofs!(at,x0.+K*Uhat.+Ubar_d)

    for x in at.X[1:N1]
        plot(x[1],x[2],"g.")
    end

    for x in at.X[I_R2]
         plot(x[1],x[2],"b.")
    end

    for x in at.X[I_R3]
        plot(x[1],x[2],"r.")
    end

    for x in at.X[I_R4]
        plot(x[1],x[2],"k.")
    end

    set_dofs!(at,x0)
end

###############################################################################
#predictor and its derivatives, relies on defining the function
# u, ∇u = u_CLE(rac, clust, _, _);

function U_CLE(at,α)
    uu = u(1.0,α)
    return [ SVector{3,Float64}([uu[:,i][:];0.0]) for i in 1:size(uu)[2] ]
end

function V_CLE(at,α)
 #   duu = central_fdm(2,1)(β -> u(1.0,β),α)
    duu = ∇u(1.0,α)
    return [ SVector{3,Float64}([duu[:,i][:];0.0]) for i in 1:size(duu)[2] ]
end

function VV_CLE(at,α)
    dduu = central_fdm(7,1)(β -> ∇u(1.0,β),α)
    return [ SVector{3,Float64}([dduu[:,i][:];0.0]) for i in 1:size(dduu)[2] ]
end
;

### alternative ways of doing it, commented out
# U_CLE(at,α) = [[SVector{2,Float64}(crk.displacements(x[1], x[2], α, 0.0, 1.0));0.0] for x in at.X]

# function V_CLE(at,α)
#     vv = SVector{3,Float64}[]
#     for x in at.X
#         ff(β) = [SVector{2,Float64}(crk.displacements(x[1], x[2], β, 0.0, 1.0));0.0]
#         gg(β) = (central_fdm(7, 1))(ff,β)
#         push!(vv,gg(α))
#     end
#     return vv
# end

# function VV_CLE(at,α)
#     vv = SVector{3,Float64}[]
#     for x in at.X
#         ff(β) = [SVector{2,Float64}(crk.displacements(x[1], x[2], β, 0.0, 1.0));0.0]
#         hh(β) = (central_fdm(11, 2))(ff,β)
# #        hh(β) = (central_fdm(7, 1))(gg,β)
#         push!(vv,hh(α))
#     end
#     return vv
# end
# ;

# function U_CLE_2(at,α)
#     uu = u(1.0,α)
#     return [ SVector{3,Float64}([uu[:,i][:];0.0]) for i in 1:size(uu)[2] ]
# end

# function V_CLE_2(at,α)
#     duu = central_fdm(7,1)(β -> u(1.0,β),α)
#  #   duu = ∇u(1.0,α)
#     return [ SVector{3,Float64}([duu[:,i][:];0.0]) for i in 1:size(duu)[2] ]
# end

# function VV_CLE_2(at,α)
#     dduu(γ) = central_fdm(11,2)(β -> u(1.0,β),γ)
#     dduu = dduu(α) #central_fdm(7,1)(duu,α)
#     return [ SVector{3,Float64}([dduu[:,i][:];0.0]) for i in 1:size(dduu)[2] ]
# end

# U_CLE_3(at,α) = [[SVector{2,Float64}(crk.displacements(x[1], x[2], α, 0.0, 1.0));0.0] for x in at.X]

# function V_CLE_3(at,α)
#     vv = SVector{3,Float64}[]
#     for x in at.X
#         Du = crk.deformation_gradient(x[1], x[2], α, 0.0, 1.0)
#         ddu = SVector{3,Float64}(-(Du[1,1]-1.0),-Du[1,2],0.0)
#         push!(vv,ddu)
#     end
#     return vv
# end

# for aa in range(-1.5,1.5,length=50)
#     @assert maximum(norm.(U_CLE(clust,aa) .- U_CLE_2(clust,aa))) < 1e-7
#     @assert maximum(norm.(U_CLE_2(clust,aa) .- U_CLE_3(clust,aa))) < 1e-7
#     @assert maximum(norm.(U_CLE_3(clust,aa) .- U_CLE_4(clust,aa))) < 1e-7
# end

# for aa in range(-1.5,1.5,length=50)
#     @assert maximum(norm.(V_CLE(clust,aa) .- V_CLE_2(clust,aa))) < 1e-4
#     @assert maximum(norm.(V_CLE_2(clust,aa) .- V_CLE_3(clust,aa))) < 1e-4
#     @assert maximum(norm.(V_CLE_3(clust,aa) .- V_CLE_4(clust,aa))) < 1e-4
# end

#useful functions for shuffling degrees of freedom (d) and the vector format (v), also for indices themselves
d2v_m1(x) = [[x[2*i-1];x[2*i];0.0] for i in 1:(length(x)÷2)]
d2v(x) = [[x[3*i-2];x[3*i-1];x[3*i]] for i in 1:(length(x)÷3)]
v2d(x) = vcat([[y[1];y[2];y[3]] for y in x]...)

Iv2Id(II) = vcat([[3*x-2;3*x-1;3*x] for x in II]...)
;

#energy, gradient and hessian associated with our NCFlex scheme:
# Ubar is in dofs format and is of size 3*N1
function f_en_m1_1(at,Ubar,K,α)
    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
#    Ubar_v = d2v_m1(vcat(Ubar,zeros(2*(N4-N1))))
#    Ubar_d = v2d(Ubar_v)
    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
    Uhat = v2d(U_CLE(at,α))
    en = energy(at,x0.+K*Uhat.+Ubar_d)
    set_dofs!(at,x0)
    return en
end

function g_en_m1_1(at,Ubar,K,α)
    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
#     Ubar_v = d2v_m1(vcat(Ubar,zeros(2*(N4-N1))))
#     Ubar_d = v2d(Ubar_v)
    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
    Uhat = v2d(U_CLE(at,α))
    grad = gradient(at,x0.+K*Uhat.+Ubar_d)
    set_dofs!(at,x0)
    return grad
end

function g_en_m1_alpha_1(at,Ubar,K,α)
    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
#     Ubar_v = d2v_m1(vcat(Ubar,zeros(2*(N4-N1))))
#     Ubar_d = v2d(Ubar_v)
#    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
#    Uhat = v2d(U_CLE(at,α))
    Vhat = v2d(V_CLE(at,α))

    II_v = [at["I_R1"]; at["I_R2"]]
    II = Iv2Id(II_v)

    g = g_en_m1_1(at,Ubar,K,α)[II]
    return dot(g,K*Vhat[II])
end

function g_en_m1_tot_1(at,Ubar,K,α)
    II1 = Iv2Id(at["I_R1"])
    return [g_en_m1_1(at,Ubar,K,α)[II1]; g_en_m1_alpha_1(at,Ubar,K,α)]
end

function h_en_m1_1(at,Ubar,K,α)
    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
#     Ubar_v = d2v_m1(vcat(Ubar,zeros(2*(N4-N1))))
#     Ubar_d = v2d(Ubar_v)
    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
    Uhat = v2d(U_CLE(at,α))
    hess = hessian2(at,x0.+K*Uhat.+Ubar_d)
    set_dofs!(at,x0)
    return hess
end

function h_en_m1_tot_1(at,Ubar,K,α)
    N4 = length(at)
    N1 = at["N1"]
    x0 = copy(dofs(at))
#     Ubar_v = d2v_m1(vcat(Ubar,zeros(2*(N4-N1))))
#     Ubar_d = v2d(Ubar_v)
    Ubar_d = vcat(Ubar,zeros(3*(N4-N1)))
    Uhat = v2d(U_CLE(at,α))
    Vhat = v2d(V_CLE(at,α))
    VVhat = v2d(VV_CLE(at,α))

    II12_v = [at["I_R1"]; at["I_R2"]]
    II12 = Iv2Id(II12_v)

    II1 = Iv2Id(at["I_R1"])

    hess = spzeros(3*N1+1,3*N1+1)
    hess_base = h_en_m1_1(at,Ubar,K,α)
    g_base = g_en_m1_1(at,Ubar,K,α)
    hess[II1,II1] = hess_base[II1,II1]

    balpha = hess_base*(K*Vhat)
    Calpha = dot(balpha[II12],K*Vhat[II12]) + dot(g_base[II12],K*VVhat[II12])
    hess[II1[end]+1,II1] = balpha[II1]
    hess[II1,II1[end]+1] = balpha[II1]
    hess[II1[end]+1,II1[end]+1] = Calpha
    return hess
end

########

### the NCF struct, here denoted by NCF_1 to not lead to clashes with older code
### TO DO: get rid of the older code, as we want to integrate fully with JuLIP anyway
mutable struct NCF_1
    at
    f
    g
    h
    Ks
    αs
    Us
    data::Dict{Any,Any}
end

function NCF_1(at)
    f = (Ubar,K,α) -> f_en_m1_1(at,Ubar,K,α)
    g = (Ubar,K,α) -> g_en_m1_tot_1(at,Ubar,K,α)
    h = (Ubar,K,α) -> h_en_m1_tot_1(at,Ubar,K,α)

    Ks = nothing
    αs = nothing
    Us = nothing
    data = Dict{Any,Any}()
    return NCF_1(at,f,g,h,Ks,αs,Us,data)
end

###########

# various useful functions for dealing with the NCFlex scheme
function Newton_static_1(ncf,xx,K,α; show = false)
    at = ncf.at
    N1 = at["N1"]
    x = copy(xx)
    If = 1:3*N1
    show && @printf("------Newton Iteration------\n")
    nit = 0
    for nit = 0:15
        ∇E = ncf.g(x,K,α)[If]
        show && @printf("%d : %4.2e \n", nit, norm(∇E[If], Inf))
        norm(∇E[If], Inf) < 1e-7 && break
        x[If] = x[1:3*N1] - ncf.h(x,K,α)[If, If] \ ∇E[If]
    end
    nit == 12 && warn("the Newton iteration did not converge")
    show && @printf("----------------------------\n")
    return x
end

function Newton_flex_1(ncf,xx,K;show=false)
    at = ncf.at
    N1 = at["N1"]
    x = copy(xx)
    If = 1:3*N1
    show && @printf("------Newton Iteration------\n")
    nit = 0
    for nit = 0:15
        ∇E = ncf.g(x[If],K,x[end])
        show && @printf("%d : %4.2e \n", nit, norm(∇E, Inf))
        norm(∇E, Inf) < 1e-6 && break
        x = x - ncf.h(x[If],K,x[end]) \ ∇E
    end
    nit == 12 && warn("the Newton iteration did not converge")
    show && @printf("----------------------------\n")
    return x
end

function alpha_to_K_1(ncf,α;K = 1.0, show=false)
    at = ncf.at
    N1 = at["N1"]
    Ubar0 = zeros(3*N1)
    fk = K -> ncf.g(Ubar0,K,α)[3*N1+1]
    gk(K) = central_fdm(2,1)(fk,K)
    show && @printf("------Newton Iteration------\n")
    nit = 0
    for nit = 0:15
        #∇E = g_en_m1_tot(x[If],K,x[end])
        show && @printf("%d : %4.2e \n", nit, norm(fk(K)))
        norm(fk(K)) < 1e-5 && break
        K = K - fk(K)/gk(K)
    end
    nit == 12 && warn("the Newton iteration did not converge")
    show && @printf("----------------------------\n")
    return K
end

function preparation_part_1(ncf;ll = 5,a1=-0.6,a2=-0.5)
    N1 = ncf.at["N1"]

    #step1: estimate (K_-,K_+) using CLE displacements only
    ks = []
    αs = range(a1,a2,length=ll)

    for α in αs
        #println(α)
        if length(ks) == 0
            push!(ks,alpha_to_K_1(ncf_t,α,K=1.0,show=true))
        else
            push!(ks,alpha_to_K_1(ncf_t,α,K = ks[end],show=true))
        end
    end
    println("(1/3) CLE-only estimation done")
    ncf.data["K_est"] = extrema(ks)
    return ks, αs
end

# function preparation_part_1a(ncf;a1=-0.4,b1=0.03,k1=0.3,k2=0.9)
#     N1 = ncf.at["N1"]

#     αs = range(a1,a1+b1,length=2)
#     ks = []
#     for α in αs
#         fffk(x) = (ncf_t.g(Ubar0,x[1],α)[3*N1+1])
#         push!(ks,find_zero(fffk, (k1,k2), FalsePosition()))
#     end
# #    ncf.Ks = [ks[1]]
# #    ncf.αs = [αs[1]]
#     return ks, αs
# end

function preparation_part_2(ncf,k,α)
    N1 = ncf.at["N1"]
    ff0(x) = ncf.f(x,k,α)
    gg0! = (gg, x) -> copyto!(gg, ncf.g(x,k,α)[1:3*N1])
    xbar0 = Optim.minimizer(optimize(ff0, gg0!, zeros(3*N1),ConjugateGradient(linesearch = BackTracking(order=2,maxstep=Inf)),Optim.Options(show_trace=true,g_tol=1e-3, f_tol=1e-32)))
    println("(2/3) First static equilibrium found")
    return xbar0
end

# function finding_k_1!(ncf,xbar,K,α_0)
#     N1 = ncf.at["N1"]
# #    xbar2 = Newton_static(ncf,xbar,K,α_0,show=false)
#     ff0(x) = ncf.f(x,K,α_0)
#     gg0! = (gg, x) -> copyto!(gg, ncf.g(x,K,α_0)[1:3*N1])
#     xbar2 = Optim.minimizer(optimize(ff0, gg0!, xbar,ConjugateGradient(linesearch = BackTracking(order=2,maxstep=Inf)),Optim.Options(show_trace=false,g_tol=1e-5, f_tol=1e-32)))
# #    xbar[1:3*N1] = xbar2[1:3*N1]
#     println("done")
#     return xbar2, ncf.g(xbar2,K,α_0)[end]
# end

#this function computes a new static equilibrium with K and alpha_0 fixed and xbar the initial guess
# it returns the new static equilibrium as well as f_alpha (to check if the new equilibrium is also a flex equilibrium)
function find_k(ncf,xbar,K,α_0)
    N1 = ncf.at["N1"]
#    xbar2 = Newton_static_1(ncf,xbar,K,α_0,show=true)
    ff0(x) = ncf.f(x,K,α_0)
    gg0! = (gg, x) -> copyto!(gg, ncf.g(x,K,α_0)[1:3*N1])
    xbar2 = Optim.minimizer(optimize(ff0, gg0!, xbar,ConjugateGradient(linesearch = BackTracking(order=2,maxstep=Inf)),Optim.Options(show_trace=true,g_tol=1e-3, f_tol=1e-32)))
    xbar3 = Newton_static_1(ncf,xbar2,K,α_0,show=true)
#    xbar[1:3*N1] = xbar2[1:3*N1]
#    print("done ")
    return xbar3, ncf.g(xbar2,K,α_0)[end]
end

#this version adjusts xbar instead of returning the solution as a separate vector
function find_k!(ncf,xbar,K,α_0)
    N1 = ncf.at["N1"]
#    xbar2 = Newton_static_1(ncf,xbar,K,α_0,show=true)
    ff0(x) = ncf.f(x,K,α_0)
    gg0! = (gg, x) -> copyto!(gg, ncf.g(x,K,α_0)[1:3*N1])
    xbar2 = Optim.minimizer(optimize(ff0, gg0!, xbar,ConjugateGradient(linesearch = BackTracking(order=2,maxstep=Inf)),Optim.Options(show_trace=true,g_tol=1e-3, f_tol=1e-32)))
    xbar3 = Newton_static_1(ncf,xbar2,K,α_0,show=true)
    xbar[1:3*N1] = xbar2[1:3*N1]
#    print("done ")
    return ncf.g(xbar2,K,α_0)[end]
end

# for this to work the choice of k1 and k2 is crucial with the corresponding f_alpha of different
# signs. Not easy to guess it, so not very efficient at the moment...
# if we choose a bad alpha for which we get the jump in f_alpha such that it misses zero, then this will never converge
function preparation_part_3(ncf,k1,k2,xbar0,α)
    xbar1 = copy(xbar0)
    fff = K-> find_k!(ncf,xbar1,K,α)
    K_best = find_zero(fff,(k1,k2),FalsePosition(),atol = 1e-5,verbose=true,maxevals=10)
    println("(3/3) First flex equilibrium found")
    return K_best, xbar1
end

############################################################################################################

# continuation routine for the CLE-only prediction
function simple_continuation_CLE(ncf ; dsmin = 0.000001, dsmax = 0.0003, ds= 0.00004,
            pMax = 4.1, maxSteps=100,theta=0.5,
            tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x, y)/length(x),
            linsolver = DefaultLS(),linearAlgo = BorderingBLS())
    N1 = ncf.at["N1"]
    function FF_m1(x, p)
        @unpack K0 = p
        return [ncf.g(Ubar0,K0,x[1])[end]]
    end
    function JJ_m1(x,p)
        @unpack K0 = p
        return hcat(ncf.h(Ubar0,K0,x[1])[end,end])
    end

    optcont1 = ContinuationPar(dsmin = dsmin, dsmax = dsmax, ds= ds, pMax = pMax, maxSteps=maxSteps,theta=theta,
        newtonOptions = NewtonPar(tol = 1e-8,linsolver = linsolver),saveSolEveryStep = 0,doArcLengthScaling = false);

    x0 = [αs[1]]
    par = (K0 = ks[1],)

    iter1 = BK.ContIterable(FF_m1, JJ_m1, x0, par, (@lens _.K0),
        optcont1; plot = false, verbosity = 0,
        tangentAlgo = tangentAlgo,dotPALC = dotPALC, linearAlgo = linearAlgo)
        #tangentAlgo = SecantPred(),dotPALC = (x, y) -> dot(x, y))

    Kdot_count = 0
    Sts = []
    dαs = []
    dKs = []

    cstate1 = nothing
    cstate2 = nothing

    next = nothing

    if ~isnothing(ncf.Ks) && length(ncf.Ks) > 1
        x0 = [ncf.αs[end-1]]
        x1 = [ncf.αs[end]]
        k0 = ncf.Ks[end-1]
        k1 = ncf.Ks[end]
        state = BK.iterateFromTwoPoints(iter1, x0, k0, x1, k1);
        next = iterate(iter1,state[1])
    else
        next = iterate(iter1)
    end

    while next !== nothing
        (i,state) = next

        ########iteration
      #  println(state.stepsizecontrol)
        print(state.step)
        print(" ")
        if state.step == 0 && state.tau.u[end] < 0.0
      #      println(state.z_old)
            state.tau.u = -state.tau.u
            state.tau.p = -state.tau.p
    #        state.stepsizecontrol = false
        end
        if ~isnothing(cstate1) && cstate1.tau.p*state.tau.p < 0.0
            if abs(state.tau.p) > 3e-5
                state.step = cstate1.step
                state.z_pred  = cstate1.z_pred
                state.tau = cstate1.tau
                state.z_old = cstate1.z_old
                state.stepsizecontrol = false
                state.ds = cstate1.ds/2
                #println("HA")
            else
                state.stepsizecontrol=true
                state.ds = iter1.contParams.ds
                Kdot_count+= 1
                println(" ")
                println("Fold point found")
                push!(ncf.αs, getx(state)[end])
                push!(ncf.Ks, getp(state))
                push!(Sts,state.step)
                push!(dαs,state.tau.u[end])
                push!(dKs,state.tau.p)
            #    push!(ncf.Us,getx(state)[1:(2*N1)])
            end
        else
            push!(ncf.αs, getx(state)[end])
            push!(ncf.Ks, getp(state))
            push!(Sts,state.step)
            push!(dαs,state.tau.u[end])
            push!(dKs,state.tau.p)
          #  push!(ncf.Us,getx(state)[1:(2*N1)])
        end
         cstate1 = copy(state)
        if Kdot_count > 6
            state.stopcontinuation = true
        end
        #################
        next = iterate(iter1,state)
    end

    ncf.data["Sts"] = Sts
    ncf.data["dαs"] = dαs
    ncf.data["dKs"] = dKs
    println(" ")
    println("End")
end

###########################################################

# full continuation routine compatible with JuLIP

function simple_continuation_1(ncf ; dsmin = 0.000001, dsmax = 0.0003, ds= 0.00004,
            pMax = 4.1, maxSteps=100,theta=0.5,
            tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x, y)/length(x),
            linsolver = DefaultLS(),linearAlgo = BorderingBLS())
    N1 = ncf.at["N1"]
    function FF_m1(x, p)
        @unpack K0 = p
        return ncf.g(x[1:3*N1],K0,x[end])
#        return [ncf.g(Ubar0,K0,x[1])[end]]
    end
    function JJ_m1(x,p)
        @unpack K0 = p
        return ncf.h(x[1:3*N1],K0,x[end])
#        return hcat(ncf.h(Ubar0,K0,x[1])[end,end])
    end

    optcont1 = ContinuationPar(dsmin = dsmin, dsmax = dsmax, ds= ds, pMax = pMax, maxSteps=maxSteps,theta=theta,
        newtonOptions = NewtonPar(tol = 1e-5,linsolver = linsolver),saveSolEveryStep = 0,doArcLengthScaling = false);

#     x0 = [αs[1]]
#     par = (K0 = ks[1],)
    x0 = [ncf.Us[1];ncf.αs[1]]
    par = (K0 = ncf.Ks[1],)

    iter1 = BK.ContIterable(FF_m1, JJ_m1, x0, par, (@lens _.K0),
        optcont1; plot = false, verbosity = 0,
        tangentAlgo = tangentAlgo,dotPALC = dotPALC, linearAlgo = linearAlgo)
        #tangentAlgo = SecantPred(),dotPALC = (x, y) -> dot(x, y))

    Kdot_count = 0
    Sts = []
    dαs = []
    dKs = []

    cstate1 = nothing
    cstate2 = nothing

    next = nothing

    if ~isnothing(ncf.Ks) && length(ncf.Ks) > 1
        x0 = [ncf.Us[end-1];ncf.αs[end-1]]
        x1 = [ncf.Us[end];ncf.αs[end]]
        k0 = ncf.Ks[end-1]
        k1 = ncf.Ks[end]

#         x0 = [ncf.αs[end-1]]
#         x1 = [ncf.αs[end]]
#         k0 = ncf.Ks[end-1]
#         k1 = ncf.Ks[end]
        state = BK.iterateFromTwoPoints(iter1, x0, k0, x1, k1);
        next = iterate(iter1,state[1])
    else
        next = iterate(iter1)
    end

    while next !== nothing
        (i,state) = next

        ########iteration
      #  println(state.stepsizecontrol)
        print(state.step)
        print(" ")
        if state.step == 0 && state.tau.u[end] < 0.0
      #      println(state.z_old)
            state.tau.u = -state.tau.u
            state.tau.p = -state.tau.p
    #        state.stepsizecontrol = false
        end
        if ~isnothing(cstate1) && cstate1.tau.p*state.tau.p < 0.0
            if abs(state.tau.p) > 3e-5
                state.step = cstate1.step
                state.z_pred  = cstate1.z_pred
                state.tau = cstate1.tau
                state.z_old = cstate1.z_old
                state.stepsizecontrol = false
                state.ds = cstate1.ds/2
                #println("HA")
            else
                state.stepsizecontrol=true
                state.ds = iter1.contParams.ds
                Kdot_count+= 1
                println(" ")
                println("Fold point found")
                push!(ncf.αs, getx(state)[end])
                push!(ncf.Ks, getp(state))
                push!(Sts,state.step)
                push!(dαs,state.tau.u[end])
                push!(dKs,state.tau.p)
            #    push!(ncf.Us,getx(state)[1:(2*N1)])
            end
        else
            push!(ncf.αs, getx(state)[end])
            push!(ncf.Ks, getp(state))
            push!(Sts,state.step)
            push!(dαs,state.tau.u[end])
            push!(dKs,state.tau.p)
            push!(ncf.Us,getx(state)[1:(3*N1)])
        end
         cstate1 = copy(state)
        if Kdot_count > 6
            state.stopcontinuation = true
        end
        #################
        next = iterate(iter1,state)
    end

    ncf.data["Sts"] = Sts
    ncf.data["dαs"] = dαs
    ncf.data["dKs"] = dKs
    println(" ")
    println("End")
end


text:
ui = plot_dict(recParams)
title = """
Reconstruction parameters
"""
content!(w, "div#content", title*ui)


text:
using QPIM
using Test

@testset "QPIM.jl" begin
    # Write your tests here.
end


text:
using IL_simulations
using Test

@testset "IL_simulations.jl" begin
    # Write your tests here.
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true) # still stakes 7+seconds, as pyimport is not sped up.
module PySide

##################################################
## initialize
using PyCall
using Mustache
using DataFrames

import Base.getindex, Base.setindex!, Base.push!, Base.pop!, Base.delete!

include("utils.jl")
include("qtutils.jl")
include("qtextras.jl")
include("data-frame-model.jl")

## replace this
##include("pyqtgraph.jl")
## call reload(Pkg.dir("PySide", "src", "pyqtgraph.jl"))
##      using PyQtGraph

export Qt, QtCore, QtSvg, QtWebkit
export qconnect, qemit, qcall, qt_enum
export qexec
export qnew_class, qnew_class_instance, qset_method
export qinvoke

export get_value, set_value, get_items, set_items, change_slot
export get_width, get_height, get_size, set_size

## widgets/controls
export Icon, Pixmap, Label, PushButton, Button,
       Slider, SpinBox,
       LineEdit, TextEdit,
       CheckBox, RadioButton, RadioGroup,
       ComboBox,
       Action, ActionGroup, MenuBar, Menu,
       GraphicsScene, GraphicsView

## models
export StandardItemModel, DataFrameModel

## views
export TableView

## container like
export MainWindow, DockWidget, Widget,
       VBoxLayout, HBoxLayout, FormLayout, GridLayout,
       TabWidget, StackedLayout, Splitter,
       ButtonGroup, GroupBox

## dialogs
export MessageBox, InputDialog, FileDialog

## methods
export windowTitle, setWindowTitle,
       text,setText,
       toPlainText, setPlainText,
       html, setHtml,
       setPlaceholderText,
       value, setValue,
       setModel,
       setIcon, setPixmap, load,
       setLayout,
       addWidget, addLayout, addTab, addRow, insertRow, setWidget,
       addAction, addMenu, addSeparator,
       setCentralWidget,
       setOrientation

export setFocus, raise

function __init__()

    global const Qt = pyimport("PySide.QtGui")
    global const QtCore = pyimport("PySide.QtCore")
    global const QtSvg = pyimport("PySide.QtSvg")
    global const QtWebKit = pyimport("PySide.QtWebKit")
    global const QtNamespace = QtCore["Qt"]
    sys = pyimport("sys")

    ## combining enums require us to work in python:
    if PyCall.pyversion < v"3"
        PyCall.pyeval("execfile(x, globals())", x = Pkg.dir("PySide", "tpl", "imports.tpl"))
    else
        PyCall.pyeval("exec(open(x, 'rb').read(), globals())", x = Pkg.dir("PySide", "tpl", "imports.tpl"))
    end

    ## use PyCall's start gui
    pygui(:qt)
    PyCall.eventloops[:qt] = PyCall.qt_eventloop("PySide", 50e-3)
    app = Qt[:QApplication](sys[:argv])

end

end


text:
module ApproximateGPs

using Reexport
@reexport using AbstractGPs
@reexport using GPLikelihoods
using Distributions
using LinearAlgebra
using Statistics
using StatsBase
using FastGaussQuadrature
using SpecialFunctions
using ChainRulesCore
using FillArrays
using KLDivergences

using AbstractGPs:
    AbstractGP,
    FiniteGP,
    LatentFiniteGP,
    ApproxPosteriorGP,
    At_A,
    diag_At_A,
    Xt_A_X,
    Xt_A_Y,
    diag_Xt_A_X

export SparseVariationalApproximation, Centered, NonCentered
export DefaultQuadrature, Analytic, GaussHermite, MonteCarlo

include("utils.jl")
include("expected_loglik.jl")
include("sparse_variational.jl")

using ForwardDiff

export LaplaceApproximation
export build_laplace_objective, build_laplace_objective!
export approx_lml  # TODO move to AbstractGPs, see https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/issues/221
include("laplace.jl")

include("deprecations.jl")

end


text:
# Kirsten Van Huffel & Tristan Vanneste

module BeesAlgorithm

# you have to import everything you need for your module to work
# if you use a new package, don't forget to add it in the package manager

# export all functions that are relevant for the user
export ArtificialBeeColonization, sphere, ackley, rosenbrock, branin, rastrigine

"""
    initialize_population(D::Number, bounds_lower::Vector, bounds_upper::Vector, Np::Number)

This function generates 'Np' random solutions (food sources) within the domain
of the variables to form an initial population for the ABC algorithm.

Input
- 'D': number of decision variables in the function to be minimized
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables
- 'Np': number of food sources/employed bees/onlooker bees

Output
- 'population': a random population of the size 'Np' with 'D' variables

## Examples

```julia-repl
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> D=4
julia> n=9
julia> initialize_population(D, bounds_lower, bounds_upper, n)
9-element Array{Any,1}:
 [-2, 2, 5, -5]
 [2, 1, -5, 1]
 [4, 1, 2, 2]
 [-3, -4, -2, -2]
 [1, 2, 4, 3]
 [4, -2, 3, -4]
 [1, 1, 0, -1]
 [2, -2, -2, -5]
 [-3, 2, 3, -3]
```
"""
function initialize_population(D::Number, bounds_lower::Vector, bounds_upper::Vector, Np::Number)
    population = []
    for i in 1:Np
        food_source = collect(rand(bounds_lower[i]:bounds_upper[i]) for i in 1:D)
        push!(population, food_source)
    end
    return population
end

"""
    compute_objective(solutions, f::Function)

Compute the objective values for an objective function 'f'. This objective value is minimized in the ABC algorithm.

Input
- 'solutions': population of solutions or a single solution
- 'f': the function that you want to use for computing objective values

Output
- 'output': objective values

## Examples

```julia-repl
julia> D=4
julia> n=9
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> population = initialize_population(D, bounds_lower, bounds_upper, n)
julia> compute_objective(population, sphere)
9-element Array{Any,1}:
 34
 75
 18
 51
 42
 51
 31
 42
 59

 julia> compute_objective(population, ackley)
 9-element Array{Any,1}:
  8.836638915350669
 11.587599478917705
  6.914978162949289
 10.20776787857683
  9.538956614465143
 10.20776787857683
  8.538932726524155
  9.538956614465143
 10.722262633846869
```
"""
function compute_objective(solutions, f::Function)
    return f.(solutions)
end

"""
    compute_fitness(objective_values)

This functions computes the fitness of each solution in the population.
The fitness is computed as ``1/(1+objective values)`` if f(x) >= 0 and
as ``1 + absolute(objective values)`` if f(x) < 0.
The bigger the objective values the smaller the fitness values.

Input
- 'objective_values': objective values

Output
- 'fitness_values': fitness values corresponding to the input values (objective values)

## Examples

```julia-repl
julia> objective_values = [4 5 -1]
julia> compute_fitness(objective_values)
3-element Array{Any,1}:
 0.2
 0.16666666666666666
 2
```
"""
function compute_fitness(objective_values)

    fitness_values = zeros(length(objective_values),1)
    fitness_values[objective_values .>= 0] = 1 ./ (1 .+ objective_values[objective_values .>= 0])
    fitness_values[objective_values .< 0] = 1 .+ abs.(objective_values[objective_values .< 0])
    return fitness_values

end

"""
    foodsource_info_prob(fitness_values)

This function measures the food source information as probabilities.
The food source information is computed as following: ``0.9 * (fitness / maximum(fitness )) + 0.1``.

Input
- 'fitness_values': fitness values

Output
- 'probabilities': probabilities of the food source information

## Examples

```julia-repl
julia> fitness_values = [0.2 0.16 1]
julia> foodsource_info_prob(fitness_values)
3-element Array{Any,1}:
 0.28
 0.24400000000000002
 1.0
```
"""
function foodsource_info_prob(fitness_values)

    probabilities = 0.9 .* (fitness_values ./ maximum(fitness_values)) .+ 0.1
    return probabilities
end

"""
    create_newsolution(solution::Vector, population, bounds_lower::Vector, bounds_upper::Vector)

Creates new solution by changing one variable using a random partner solution.
Firstly, a random variable in the solution to be modified will be selected.
Secondly,  a random partner in the current population is chosen, and a random
variable of this partner solution is picked. Based on the value of these two variables,
a new value is interpolated and assigned to the first random variable of the original solution.

Input
- 'solutions': solution
- 'population': population of solutions
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables

Output
- 'solution_new': new solution with one variable changed using a random variable in a random partner solution

## Examples

```julia-repl
julia> solution = [4 0 1 4]
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> D=4
julia> n=9
julia> population = initialize_population(D, bounds_lower, bounds_upper, n)
julia> create_newsolution(solution, population, bounds_lower, bounds_upper)
1×4 Array{Float64,2}:
 -1.49748  0.0  1.0  4.0
```
"""
function create_newsolution(solution::Vector, population, bounds_lower::Vector, bounds_upper::Vector)
    # select random variable to change
    randomvar1_index = rand(1:size(solution)[1], 1)

    # select partner solution to generate new solution
    randompartner_index = rand(1:size(population)[1], 1)

    # select random variable in partner solution to interpolate with

    randompartner = population[randompartner_index, :][1]
    randomvar2_index = rand(1:length(randompartner), 1)

    # create new food location
    phi = rand()*2-1 #random number between -1 and 1
    global solution_new = float(deepcopy(solution))
    a = solution[randomvar1_index]
    b = randompartner[randomvar2_index]
    solution_new[randomvar1_index] = a + phi*(a - b)

    # check if lower bound is violated
    if solution_new[randomvar1_index] < bounds_lower[randomvar1_index]
        solution_new[randomvar1_index] = bounds_lower[randomvar1_index]
    end

    # check if upper bound is violated
    if solution_new[randomvar1_index] > bounds_upper[randomvar1_index]
        solution_new[randomvar1_index] = bounds_upper[randomvar1_index]
    end
return solution_new
end

"""
    employed_bee_phase(population, bounds_lower::Vector, bounds_upper::Vector,
                            trial::Array, Np::Number, f::Function)

In this function, the employed bee phase is implemented.
Employed bees try to identify better food source than the one they were associated previously.
A new solution is generated using a partner solution.
Thereafter, greedy selection is performed, meaning that a new solution only  will be accepted if it is better than the current solution.
Every bee in the swarm will explores one food source. All solutions get an opportunity to generate a new solution in the employed bee phase.

Input
- 'population': population of solutions
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables
- 'trial': Vector of current trials of solutions in the population
- 'Np': number of food sources/employed bees/onlooker bees
- 'f': the function that you want to use for computing objective values

Output
- 'population new evolved': newly generated population of solutions
- 'fitness new evolved': fitness values of the new population
- 'objective new evolved': objective values of the new population
- 'trial': updated trials of new solutions in the new population
    When original solution has failed to generate better solution, trial counter is increased by 1 unit
    When better solution has been found, the trial counter for this new solution is set to zero

## Examples

```julia-repl
julia> trial = zeros(size(population)[1])
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> D=4
julia> N=9
julia> population = initialize_population(D, bounds_lower, bounds_upper, n)
julia> population_new_evolved, fitness_new_evolved, objective_new_evolved, trial = employed_bee_phase(population, bounds_lower, bounds_upper, trial,n ,sphere)
julia> population_new_evolved
9-element Array{Any,1}:
 [5.0, -5.0, 5.0, 2.748420414594148]
 [2, 4, 0, -4]
 [-5.0, 2.0, 0.0, -1.1761590179019343]
 [1, 5, 5, -2]
 [2.0, -1.0, -3.0, 0.75942632274392]
 [1, 5, -1, 1]
 [0.0, 3.1592990065770192, -2.0, -3.0]
 [-4, 3, 3, 2]
 [3, 2, -3, -3]
```
"""
function employed_bee_phase(population, bounds_lower::Vector, bounds_upper::Vector,
                                trial::Array, Np::Number, f::Function)
    population_new = []

    # create new food sources
    for i in 1:Np
        solution = population[i, :][1]
        solution_new = solution
        while solution_new == solution
            solution_new = create_newsolution(solution, population, bounds_lower, bounds_upper)
        end
        append!(population_new, [solution_new])
    end

    # evaluate fitness old and new population
    objective_values_old = compute_objective(population, f)
    fitness_old = compute_fitness(objective_values_old)
    objective_values_new = compute_objective(population_new, f)
    fitness_new = compute_fitness(objective_values_new)

    # perform greedy selection
    population_new_evolved = []
    fitness_new_evolved = []
    objective_new_evolved = []

    for j in 1:Np
        if fitness_new[j] > fitness_old[j]
            append!(population_new_evolved, [population_new[j]])
            append!(fitness_new_evolved, fitness_new[j])
            append!(objective_new_evolved, objective_values_new[j])
            trial[j] = 0
        else
            append!(population_new_evolved, [population[j]])
            append!(fitness_new_evolved, fitness_old[j])
            append!(objective_new_evolved, objective_values_old[j])
            trial[j] += 1
        end
    end

    return population_new_evolved, fitness_new_evolved, objective_new_evolved, trial
end

"""
    onlooker_bee_phase(population, bounds_lower::Vector,
                                    bounds_upper::Vector,
                                    trial::Array, Np::Number,
                                    f::Function)

This function implements the onlooker bee phase.
In the onlooker bee phase, a food source is selected for further exploitation with a probability related to the nectar amount,
i.e. a solution with higher fitness will have a higher probability to be chosen.
Fitter solutions may undergo multiple onlooker bee explorations. As in the employed bee phase,
new solutions are generated using a partner solution and greedy selection is performed.
In contrast to the employed bee phase, not every food source will be explored,
since every onlooker bee will explore a certain food source with a certain probability (depending on nectar amount).

Input
- 'population': population of solutions
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables
- 'trial': Vector of current trials of solutions in the population
- 'Np': number of food sources/employed bees/onlooker bees
- 'f': the function that you want to use for computing objective values

Output
- 'population new evolved': newly generated population of solutions
- 'fitness new evolved': fitness values of the new population
- 'objective new evolved': objective values of the new population
- 'trial': updated trials of new solutions in the new population
    When original solution has failed to generate better solution, trial counter is increased by 1 unit
    When better solution has been found, the trial counter for this new solution is set to zero

## Examples

```julia-repl
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> D=4
julia> N=9
julia> population = initialize_population(D, bounds_lower, bounds_upper, N)
julia> trial = zeros(size(population)[1])
julia> population_new_evolved, fitness_new_evolved, objective_new_evolved, trial = onlooker_bee_phase(population, bounds_lower, bounds_upper, trial,N ,sphere)
julia> population_new_evolved
9-element Array{Any,1}:
 [-1.0, -5.0, 0.5313751547457599, 5.0]
 [-1, -3, -4, -2]
 [0, 1, -3, 5]
 [1.0, -4.0, -5.0, 1.8770812747894876]
 [-0.944880850158663, -1.2865870734899967, 0.04370050469571929, 1.0]
 [4, -3, 5, 5]
 [0.0, 1.0208794819836884, 0.0, -3.0]
 [-2, 2, 5, -1]
 [-5, -3, -4, 5]
```
"""
function onlooker_bee_phase(population, bounds_lower::Vector, bounds_upper::Vector,
                    trial::Array, Np::Number, f::Function)
    m = 0 # onlooker bee
    n = 1 # food source

    objective_values = compute_objective(population,f)
    fitness = compute_fitness(objective_values)
    # first calculate the probability values
    proba = foodsource_info_prob(fitness)

    while m <= Np # we want for every onlooker bee a new solution
        r = rand()
        if r <= proba[n]
            solution = population[n, :][1] # solution n
            objective_values_old = compute_objective([solution], f)
            fitness_old = compute_fitness(objective_values_old)[1]

            solution_new = solution
            while solution_new == solution
                solution_new = create_newsolution(solution, population, bounds_lower, bounds_upper)
            end

            objective_values_new = compute_objective([solution_new], f)
            fitness_new = compute_fitness(objective_values_new)[1]
            if fitness_new > fitness_old # if this get accepted
                population[n, :] = [solution_new]
                trial[n]=0
            else
                trial[n] += 1
            end
            m = m + 1
        end
        # if the rand < proba is not sattisfied
        n = n +1
        if n > Np
            n = 1
        end
    end
    objective_new_evolved = compute_objective(population,f)
    fitness_new_evolved = compute_fitness(objective_new_evolved)

    return population, fitness_new_evolved, objective_new_evolved, trial
end

"""
    scouting_phase(population, bounds_lower::Vector, bounds_upper::Vector,
                    D::Number ,trials::Array, fitness, objective, limit::Number, f::Function)

This function implements the scouting phase.
If the value of the trial counter for a certain solution is greater than fixed limit, then a solution can enter the scout phase.
The latter food source is then considered as exhausted and will therefore be abandoned by the bees.
The current optimal solution will be stored separately to avoid losing it in the scouting phase.
After discarding the exhausted solution, a new random solution is generated and the trial counter of this solution is reset to zero.

Input
- 'population': population of solutions
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables
- 'D': number of decision variables in the function to be minimized
- 'trials': current trial of solutions
- 'fitness': fitness values
- 'objective': objective values
- 'limit': limit value
- 'f': the function that you want to use for computing objective values

Output
- 'population new evolved': newly generated population of solutions
- 'fitness new evolved': fitness values of the new population
- 'objective new evolved': objective values of the new population
- 'trial': updated trials of new solutions in the new population
    When original solution has failed to generate better solution, trial counter is increased by 1 unit
    When better solution has been found, the trial counter for this new solution is set to zero

## Examples

```julia-repl
julia> trial = ones(size(population)[1])
julia> bounds_lower = [-5,-5,-5,-5]
julia> bounds_upper = [5,5,5,5]
julia> D=4
julia> N=9
julia> population = initialize_population(D, bounds_lower, bounds_upper, n)
julia> population
9-element Array{Any,1}:
 [3, 1, 2, -2]
 [-5, 4, 4, 2]
 [-4, 5, -2, -4]
 [5, 5, -5, 0]
 [-2, -4, 5, 5]
 [1, 2, 5, -2]
 [1, -1, 5, -5]
 [-2, -5, -4, -5]
 [5, -1, -1, -1]
julia> objective= compute_objective(population,sphere)
julia> fitness = compute_fitness(objective)
julia> population_new_evolved, fitness_new_evolved, objective_new_evolved, trial  = scouting_phase(population, bounds_lower, bounds_upper, D, trial, fitness, objective, 0, sphere)
julia> population_new_evolved
9-element Array{Any,1}:
 [3, 1, 2, -2]
 [-5, 4, 4, 2]
 [-2.959166963468496, 2.2098469985145712, -0.2708630367391418, 4.89872222904274]
 [5, 5, -5, 0]
 [-2, -4, 5, 5]
 [1, 2, 5, -2]
 [1, -1, 5, -5]
 [-2, -5, -4, -5]
 [5, -1, -1, -1]
```
"""
function scouting_phase(population, bounds_lower::Vector, bounds_upper::Vector,
                        D::Number ,trials::Array, fitness, objective, limit::Number, f::Function)

    # check whether the trial vector exceed the limit value and importantly where
    index_exceed = trials .> limit
    if sum(index_exceed) >= 1 # there is minimal one case where we exceed the limit
        if sum(maximum(trials) .== trials) > 1 # multiple cases have the same maximum so chose randomly

            possible_scoutings = findall(trials .== maximum(trials))
            idx = rand(1:size(possible_scoutings)[1])
            scouting_array = possible_scoutings[idx]
        else # only one array has a maximum => chose this one

            scouting_array = argmax(trials)

        end
        pop = population[scouting_array]
        fit = fitness[scouting_array]
        obj = objective[scouting_array]
        trail = trials[scouting_array]

        #creating random population
        sol_new = bounds_lower + (bounds_upper-bounds_lower) .* rand(D) # -5 *(10*rand)
        new_obj = compute_objective([sol_new],f)
        new_fit = compute_fitness(new_obj)

        # replacing the new population
        population[scouting_array] = sol_new

        fitness[scouting_array] = new_fit[1]
        objective = float(objective)
        objective[scouting_array] = new_obj[1]
        trials[scouting_array] = 0

    end
    return population, fitness, objective, trials
end

"""
    ArtificialBeeColonization(D::Number, bounds_lower::Vector, bounds_upper::Vector,
                                S::Number, T::Number, limit::Number, f::Function)

This functions runs the Artificial Bee Colony Algorithm.
First a initialize population is made. For each iteration (T) the ABC algorithm undergoes the employed bee phase,
onlooker bee phase and possibly the scouting phase.
During the process, the solution with the highest fitness value is stored separately and updated after every iteration and scouting phase.

Input
- 'D': number of decision variables in the function to be minimized
- 'bounds_lower': lower bounds of variables
- 'bounds_upper': upper bounds of variables
- 'S': swarm size (S/2 = Np, the number of employed bees, onlooker bees and scouting bees)
- 'T': number of cycles
- 'limit': decides when scouts phase needs to be executed (often taken Np*D)
- 'f': the function that you want to use for computing objective values (this function will be minimized)

Output
- 'optimal_solution': gives a vector of the size D with the optimal solution
- 'populations': all populations that were computed during the ABC algorithm
- 'best fitness tracker': a vector with for each iteration the best fitness value so far encountered.

## Examples

```julia-repl
julia> S = 24
julia> bounds_lower = [-100,-100,-100,-100]
julia> bounds_upper = [100,100,100,100]
julia> D=4
julia> limit = D * (S/2)
julia> T = 500
julia> optimal_solution,populations,best_fitness_tracker = ArtificialBeeColonization(D, bounds_lower, bounds_upper, S, T, limit, sphere)
julia> optimal_solution
4-element Array{Float64,1}:
  1.4537737816170574e-9
 -9.379333377223156e-9
 -2.478809612571768e-9
  3.9716934385662925e-11
```
"""
function ArtificialBeeColonization(D::Number, bounds_lower::Vector, bounds_upper::Vector,
                                    S::Number, T::Number, limit::Number, f::Function)
    @assert D > 0 "D must be positive"
    @assert bounds_lower <= bounds_upper "Lower bounds must be smaller than upper bounds"
    @assert length(bounds_lower) == length(bounds_upper) == D  "The length of the lower bounds must be equal to the length of the upper bounds and the number of decision variables"
    @assert iseven(S) "The particle swarm size must be an even number"
    @assert S > 0 "Particle swarm size must be a positive number"
    @assert T > 0 "Number of cycles must be positive"

    Np = Int8(S/2) # number of food sources/employed bees/onlooker bees

    # initialize population
    population = initialize_population(D, bounds_lower, bounds_upper, Np)

    # calculate objective values and fitness values for population
    objective_values = compute_objective(population, f)
    fitness_values = compute_fitness(objective_values)

    # initialize trial vector for population
    trial = zeros(Np, 1)
    best_fitness = 0
    best_fitness_tracker = []
    optimal_solution = []
    populations = []
    for iterations in 1:T

        ## EMPLOYED BEE PHASE
        population, fitness_values, objective_values, trial = employed_bee_phase(population,
                                                                    bounds_lower, bounds_upper, trial, Np, f)

        ## ONLOOKER BEE PHASE
        population, fitness_values, objective_values, trial = onlooker_bee_phase(population,
                                                                    bounds_lower, bounds_upper, trial, Np, f)

        ## SCOUTING PHASE
        if maximum(fitness_values) > best_fitness
            best_fitness = maximum(fitness_values)
            ind = argmax(fitness_values)
            optimal_solution = population[ind]

        end

        population, fitness_values, objective_values, trial = scouting_phase(population, bounds_lower, bounds_upper,
                                                                    D,trial, fitness_values, objective_values, limit, f)

        if maximum(fitness_values) > best_fitness
            best_fitness = maximum(fitness_values)
            ind = argmax(fitness_values)
            optimal_solution = population[ind]

        end
        populations = append!(populations, [population])
        best_fitness_tracker = append!(best_fitness_tracker, best_fitness)
    end

    return optimal_solution, populations, best_fitness_tracker
end

"""
    sphere(x)

This is computing the sphere function value for the vector x.

Input
- 'x': input vector for the sphere function

Output
- 'output': output value for the sphere function

## Examples

```julia-repl
julia> sphere(4)
16
julia> sphere([4,5])
41
```
"""
function sphere(x)
    return sum(x.^2)
end

"""
    ackley(x; a=20, b=0.2, c=2π)

This is computing the ackley function value for the vector x.

Input
- 'x': input vector for the ackley function

Output
- 'output': output value for the ackley function

## Examples

```julia-repl
julia> ackley(4)
11.013420717655569
julia> ackley([4,5])
11.913518152857637
```
"""
function ackley(x; a=20, b=0.2, c=2π)
    d = length(x)
    return -a * exp(-b*sqrt(sum(x.^2)/d)) -
        exp(sum(cos.(c .* x))/d) + a + exp(1)
end

"""
    rosenbrock(x; a=1, b=5)

This is computing the rosenbrock function value for the vector x.
Watch out! This function always needs a 2-element Array as input

Input
- 'x': input vector for the rosenbrock function

Output
- 'output': output value for the rosenbrock function

## Examples

```julia-repl
julia> rosenbrock([4,5])
12109
```
"""
function rosenbrock(x; a=1, b=100)
    # 2 dimensions!
    @assert length(x) == 2 "The rosenbrock function is restricted to two decision variables"
    return (a-x[1])^2 + b*(x[2]-x[1]^2)^2
end

"""
    branin((x1, x2); a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)

This is computing the branin function value for the vector x.
Watch out! This function always needs a 2-element Array as input

Input
- 'x': input vector for the branin function

Output
- 'output': output value for the branin function

## Examples

```julia-repl
julia> branin([4,5])
14.608661704375713
```
"""
function branin(x; a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)
    # 2 dimensions!
    @assert length(x) == 2 "The branin function is restricted to two decision variables"
    return a * (x[2] - b * x[1]^2 + c * x[1] - r)^2 + s * (1 - t) * cos(x[1]) + s
end
"""
    rastrigine(x; A=10)

This is computing the rastrigine function value for the input values of x.

Input

- 'x': input vector for the rastrigine function

Output

- 'output': output value for the rastrigine function

## Examples

```julia-repl
julia> rastrigine(4)
16.0
julia> rastrigine([4,5])
41.0
```
"""
function rastrigine(x; A=10)
    return length(x) * A + sum(x.^2 .- A .* cos.(2pi .* x))
end

end

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
function _f(a::Array, x::Function)
    s = copy(a)
    n = length(s)
    p = 1
    i = 0
    while p <= n
        i += 1
        value = s[p]
        s[p] = x(value)
        p += value
    end
    return i
end

function puzzle05(path::String=joinpath(@__DIR__, "..", "data/05.txt"))
    input = parse.(Int, readlines(path))
    solution = map(x->_f(input, x), [x -> x+1, x -> x + (x≥3 ? -1 : 1)])
    return solution
end


text:
@testset "Particle Swarm" begin
    # TODO: Run on MultivariateProblems.UnconstrainedProblems?
    Random.seed!(100)

    function f_s(x::Vector)
        (x[1] - 5.0)^4
    end

    function rosenbrock_s(x::Vector)
        (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    end

    initial_x = [0.0]
    upper = [100.0]
    lower = [-100.0]
    n_particles = 4
    options = Optim.Options(iterations=100)
    res = Optim.optimize(f_s, initial_x, ParticleSwarm(lower, upper, n_particles),
                         options)
    @test norm(Optim.minimizer(res) - [5.0]) < 0.1

    initial_x = [0.0, 0.0]
    lower = [-20., -20.]
    upper = [20., 20.]
    n_particles = 5
    options = Optim.Options(iterations=300)
    res = Optim.optimize(rosenbrock_s, initial_x, ParticleSwarm(lower, upper, n_particles),
                             options)
    @test norm(Optim.minimizer(res) - [1.0, 1.0]) < 0.1
    options = Optim.Options(iterations=300, show_trace=true, extended_trace=true, store_trace=true)
    res = Optim.optimize(rosenbrock_s, initial_x, ParticleSwarm(lower, upper, n_particles), options)
    @test summary(res) == "Particle Swarm"
    res = Optim.optimize(rosenbrock_s, initial_x, ParticleSwarm(n_particles = n_particles), options)
    @test summary(res) == "Particle Swarm"
end


text:
module NASAPowerDownloader

using HTTP
using JSON

function download_regional(year, regional, mode, timeout=30)

    try
        r = HTTP.request(
            "GET",
            "https://power.larc.nasa.gov/api/temporal/$(mode)/regional?start=$(year)0101&end=$(year)1231&latitude-min=$(regional.latitude_min)&latitude-max=$(regional.latitude_max)&longitude-min=$(regional.longitude_min)&longitude-max=$(regional.longitude_max)&community=re&parameters=ALLSKY_SFC_SW_DWN,T2M,T2M_MIN,T2M_MAX,RH2M,WS10M,WS10M_MIN,WS10M_MAX,WD10M,PS&time-standard=utc&format=json&header=true",
            readtimeout = timeout
        )
        j = JSON.parse(String(r.body))
        return j
    catch e
        # Too Many Requests
        if (e.status == 429) || (typeof(e) == HTTP.TimeoutRequest.ReadTimeoutError)
            sleep(5)
            # println("Retrying download ♻")
            return download_regional(year, regional, mode, timeout)
        else
            print("Error $(e)")
            exit(-1)
        end
    end
end

function download_point(year, point, mode, timeout=30)

    try
        r = HTTP.request(
            "GET",
            "https://power.larc.nasa.gov/api/temporal/$(mode)/point?start=$(year)0101&end=$(year)1231&latitude=$(point.latitude)&longitude=$(point.longitude)&community=re&parameters=ALLSKY_SFC_SW_DWN&time-standard=utc&format=json&header=true",
            readtimeout = timeout
        )
        j = JSON.parse(String(r.body))
        return j
    catch e
        # Too Many Requests
        if (e.status == 429) || (typeof(e) == HTTP.TimeoutRequest.ReadTimeoutError)
            sleep(5)
            # println("Retrying download ♻")
            return download_point(year, point, mode, timeout)
        else
            print("Error $(e)")
            exit(-1)
        end
    end
end
end


text:
for (elty, relty, ext) in ((:ElInt, :ElInt, :i),
                           (:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:ComplexF32, :Float32, :c),
                           (:ComplexF64, :Float64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"),
                       (:DistMultiVec, "DistMultiVec_"))
        @eval begin
            # Bernoulli
            function bernoulli!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1, p::Real = 0.5)
                ElError(ccall(($(string("ElBernoulli", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt, Float64),
                    A.obj, m, n, Float64(p)))
                return A
            end
            bernoulli(::Type{$mat{$elty}}, m::Integer, n::Integer = 1, p::Real = 0.5) =
                bernoulli!($mat($elty), m, n, p)

            # Gaussian
            function gaussian!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1,
                               mean::Number = 0, stddev::Number = 1)
                ElError(ccall(($(string("ElGaussian", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt, $elty, $relty),
                    A.obj, m, n, mean, stddev))
                return A
            end
            gaussian(::Type{$mat{$elty}}, m::Integer, n::Integer = 1,
                mean::Number = 0, stddev::Number = 1) =
                    gaussian!($mat($elty), m, n, mean, stddev)

            # Ones
            function ones!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1)
                ElError(ccall(($(string("ElOnes", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt),
                    A.obj, m, n))
                return A
            end
            ones(::Type{$mat{$elty}}, m::Integer, n::Integer = 1) = ones!($mat($elty), m, n)

            # Uniform
            function uniform!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1,
                              center::Number = 0, radius::Real = 1)
                ElError(ccall(($(string("ElUniform", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt, $elty, $relty),
                    A.obj, m, n, center, radius))
                return A
            end
            uniform(::Type{$mat{$elty}}, m::Integer, n::Integer = 1, center::Number = 0, radius::Real = 1) =
                uniform!($mat($elty), m, n, center, radius)

            # Zeros
            function zeros!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1)
                ElError(ccall(($(string("ElZeros", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt),
                    A.obj, m, n))
                return A
            end
            zeros(::Type{$mat{$elty}}, m::Integer, n::Integer = 1) = zeros!($mat($elty), m, n)
        end

        if elty == :ComplexF32 || elty == :ComplexF64
            # Uniform
            @eval begin
                function foxLi!(A::$mat{$elty}, n::Integer = size(A, 1), ω::Real = 1.0)
                    ElError(ccall(($(string("ElFoxLi", sym, ext)), libEl), Cuint,
                        (Ptr{Cvoid}, ElInt, $relty),
                        A.obj, n, ω))
                    return A
                end
                foxLi(::Type{$mat{$elty}}, m::Integer, n::Integer = 1, ω::Real = 1.0) = foxLi!($mat($elty), m, n, ω)
            end
        end
    end
end

for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:ComplexF32, :Float32, :c),
                           (:ComplexF64, :Float64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"),
                       (:SparseMatrix, "Sparse_"),
                       (:DistSparseMatrix, "DistSparse_"))
        @eval begin
            # Helmholtz
            function helmholtz!(A::$mat{$elty}, nx::Integer; shift::Number = 0)
                ElError(ccall(($(string("ElHelmholtz1D", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, $elty),
                    A.obj, nx, shift))
                return A
            end

            function helmholtz!(A::$mat{$elty}, nx::Integer, ny::Integer; shift::Number = 0)
                ElError(ccall(($(string("ElHelmholtz2D", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt, $elty),
                    A.obj, nx, ny, shift))
                return A
            end

            function helmholtz!(A::$mat{$elty}, nx::Integer, ny::Integer, nz::Integer; shift::Number = 0)
                ElError(ccall(($(string("ElHelmholtz3D", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt, ElInt, ElInt, $elty),
                    A.obj, nx, ny, nz, shift))
                return A
            end
            helmholtz(::Type{$mat{$elty}}, n::Integer...; shift::Number = 0) = helmholtz!($mat($elty), n..., shift = shift)
        end
    end
end

for (elty, ext) in ((:ComplexF32, :c),
                    (:ComplexF64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"))
        @eval begin
            # Fourier
            function fourier!(A::$mat{$elty}, n::Integer)
                ElError(ccall(($(string("ElFourier", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, ElInt),
                    A.obj, n))
                return A
            end
            fourier(::Type{$mat{$elty}}, n::Integer) = fourier!($mat($elty), n)
        end
    end
end


text:
using CSDP, JuMP
optim = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

@testset "Spectral radius" begin
    @test spectral_radius(I(2)) == 1.0
    A = [1 2; 3 0]
    @test spectral_radius(A) == 3.0
    A = [1 0 1; 0 1 0; 1 0 1]
    @test spectral_radius(A) == 2.0
end

E1, E2, E3, E4 = [1 0.; 0 1.], [1 0.; 0 -1.], [-1 0.; 0 1.], [-1 0.; 0 -1.]
D1 = [0. 0.]; D2, D3, D4 = copy(D1), copy(D1), copy(D1)
X = [[E1, D1], [E2, D2], [E3, D3], [E4, D4]]

Σ1 = [[1. -1.; 2. 1.], [1. -1.; 1. 1.], [3. -2.; 1. 0.], [-2. 0.; 1. -1.]]
G1 = automaton_constructor(Σ1)
s1 = discreteswitchedsystem(Σ1, G1, X)
Σ2 = [[1.34 -1.34; 1. 0.], [1.34 -0.67; 1. 0.], [0.67 -1.34; 1. 0.], [0.67 -0.67; 1. 0.]]
G2 = automaton_constructor(Σ2)
s2 = discreteswitchedsystem(Σ2, G2, X)
Σ3 = [[1. -5.; 2. 0.], [2. -3.; 1. 0.], [0 -1.; 0. 1.], [2. 3.; -1. -4.]]
G3 = automaton_constructor(Σ3)
s3 = discreteswitchedsystem(Σ3, G3, X)

γsdp1 = sdpbound_γ(s1)
gammasdp1 = sdpbound_gamma(s1)
γsos1 = sosbound_γ(s1, 2)
gammasos1 = sosbound_gamma(s1, 2)
γsdp2 = sdpbound_γ(s2, optimizer=optim, verbose=2)
γsos2 = sosbound_γ(s2, 2, optimizer=optim, verbose=2)
γsdp3 = sdpbound_γ(s3, optimizer=optim, verbose=4)
γsos3 = sosbound_γ(s3, 2, optimizer=optim, verbose=4)

@testset "SDP program" begin
    @test γsdp1 == gammasdp1
    @test γsos1 == gammasos1
    @test (γsdp1 ≥ 2 && γsdp1 ≤ 2.1)
    @test (γsdp2 ≥ 1.1 && γsdp2 ≤ 1.2)
    @test (γsdp3 ≥ 3.1 && γsdp3 ≤ 3.2)
end

@testset "SOS program" begin
    @test (γsos1 ≥ 2 && γsos1 ≤ 2.1)
    @test (γsos2 ≥ 1.1 && γsos2 ≤ 1.2)
    @test (γsos3 ≥ 3.1 && γsos2 ≤ 3.2)
end

@testset "(C)JSR calculation" begin
    @test jsr(s1) == jsr(Σ1)
    @test cjsr(s1) == cjsr(Σ1, G1)
end


text:
# created with the help of SnoopCompile.jl
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(internal_rate_of_return),Vector{Int64}})
    Base.precompile(Tuple{typeof(internal_rate_of_return),Vector{Float64}})
    Base.precompile(Tuple{typeof(internal_rate_of_return),Vector{Int64},Vector{Float64}})
    Base.precompile(Tuple{typeof(internal_rate_of_return),Vector{Float64},Vector{Float64}})

    Base.precompile(Tuple{typeof(present_value),Float64,Vector{Int64},Vector{Float64}})
    Base.precompile(Tuple{typeof(present_value),Float64,Vector{Float64},Vector{Float64}})
    Base.precompile(Tuple{typeof(present_value),Float64,Vector{Float64},Vector{Int64}})

    Base.precompile(Tuple{typeof(duration),Float64,Vector{Float64}})
    Base.precompile(Tuple{typeof(duration),Float64,Vector{Int64}})

    Base.precompile(Tuple{typeof(convexity),Float64,Vector{Float64}})
    Base.precompile(Tuple{typeof(convexity),Float64,Vector{Int64}})

end

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
# list

export ZeCommandList, execute!

mutable struct ZeCommandList
    handle::ze_command_list_handle_t

    context::ZeContext
    device::ZeDevice

    function ZeCommandList(ctx::ZeContext, dev::ZeDevice, ordinal=1; flags=0)
        desc_ref = Ref(ze_command_list_desc_t(;
            commandQueueGroupOrdinal=ordinal-1, flags,
        ))
        handle_ref = Ref{ze_command_list_handle_t}()
        zeCommandListCreate(ctx, dev, desc_ref, handle_ref)
        obj = new(handle_ref[], ctx, dev)
        finalizer(obj) do obj
            zeCommandListDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_command_list_handle_t}, list::ZeCommandList) = list.handle

Base.close(list::ZeCommandList) = zeCommandListClose(list)

Base.reset(list::ZeCommandList) = zeCommandListReset(list)

"""
    ZeCommandList(dev::ZeDevice, ...) do list
        append_...!(list)
    end

Create a command list for device `dev`, passing in a do block that appends operations.
The list is then closed and can be used immediately, e.g. for execution.

"""
function ZeCommandList(f::Base.Callable, args...; kwargs...)
    list = ZeCommandList(args...; kwargs...)
    f(list)
    close(list)
    return list
end

execute!(queue::ZeCommandQueue, lists::Vector{ZeCommandList}, fence=nothing) =
    zeCommandQueueExecuteCommandLists(queue, length(lists), lists, something(fence, C_NULL))

"""
    execute!(queue::ZeCommandQueue, ...) do list
        append_...!(list)
    end

Create a command list for the device that owns `queue`, passing in a do block that appends
operations. The list is then closed and executed on the queue.
"""
function execute!(f::Base.Callable, queue::ZeCommandQueue, fence=nothing; kwargs...)
    list = ZeCommandList(f, queue.context, queue.device, queue.ordinal; kwargs...)
    execute!(queue, [list])
end


text:
using InsarLOS
using Test

@testset "InsarLOS.jl" begin

    dbfile = joinpath(@__DIR__, "S1A_IW_SLC__1SDV_20150316T005120_20150316T005148_005050_006570_8AB8.SAFE.db.1")
    lat = 31.0
    lon = -103.0
    xyz = [0.72732057375709702; 0.60360651252993913; -0.32659449018789627]
    @test isapprox(InsarLOS.calculate_los_xyz(lat, lon, dbfile=dbfile), xyz, atol=1e-3)
    enu = [0.572897473168904; 0.107232543160571; -0.812582098574662]
    @test isapprox(InsarLOS.rotate_xyz_to_enu(xyz, lat, lon), enu, atol=1e-3)

    # Combo function
    @test isapprox(InsarLOS.get_los_enu([lat, lon], dbfile=dbfile), enu, atol=1e-3)

    lat = 31.0
    lon = -102.0
    xyz = [0.78237765815518601; 0.55250630461773964; -0.28744039969575685]
    @test isapprox(InsarLOS.calculate_los_xyz(lat, lon, dbfile=dbfile), xyz, atol=1e-3)
    enu = [0.650408309211683; 0.115737820558811; -0.750715517490279]
    @test isapprox(InsarLOS.rotate_xyz_to_enu(xyz, lat, lon), enu, atol=1e-3)

    # InsarLOS._compute_xyz_ground(31., -102.)
    # xyzground = [-1.1378e6; -5.3529e6; 3.26624e6];

#     enumat = hcat([0.75734; 0.20000; -0.62163], [0.70140; 0.19368; -0.68559])
#     llmat = hcat([31.; -103.], [31.; -102.])
#     @show enumat, llmat
#     @test isapprox(InsarLOS.get_los_enu(llmat, dbfile=dbfile), enumat, atol=1e-3)
end


text:
#blend
#copy
#filter
#get
#loadPixels
#set
#updatePixels


text:
"""
    mgradient(img; mode=:beucher, dims=coords_spatial(img), r=1)
    mgradient(img, se; mode=:beucher)

Calculate morphological gradient of the image using given mode.

There are three widely used modes[1]:

- `:beucher`: the default mode. It calculates the arithmetic difference between the dilation
  and the erosion -- `dilate(img, se) - erode(img, se)`.
- `:internal`: also known as _half-gradient by erosion_. It calculates the arithmetic
  difference between the original image and its erosion -- `img - erode(img, se)`.
- `:external`: also known as _half-gradient by dilation_. It calculates the arithmetic
  difference between dilation and the original image -- `dilate(img, se) - se`.

$(_docstring_se)

# Examples

```jldoctest; setup = :(using ImageMorphology)
julia> img = falses(7, 7); img[3:5, 3:5] .= true; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img)) # default mode :beucher always creates a two-pixel wide boundary
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  1  1  1  1  1  0
 0  1  1  1  1  1  0
 0  1  1  0  1  1  0
 0  1  1  1  1  1  0
 0  1  1  1  1  1  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img; mode=:internal)) # half-gradient -- the boundary is internal to original image
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img; mode=:external)) # half-gradient -- the boundary is external to original image
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  1  1  1  1  1  0
 0  1  0  0  0  1  0
 0  1  0  0  0  1  0
 0  1  0  0  0  1  0
 0  1  1  1  1  1  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  1  1  1  1  1  0
 0  1  1  0  1  1  0
 0  1  1  1  1  1  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
```

The beucher operator is a self-complementary operator in the sense that `mgradient(img, se;
mode=:beucher) == mgradient(complement.(img), se; mode=:beucher)`. When `r>1`, it is usually
called _thick gradient_. If a line segment is used as `se`, then the gradient becomes the
_directional gradient_.

## See also

- [`mgradient!`](@ref) is the in-place version of this function.
- [`mlaplacian`](@ref) for the laplacian operator.
- `ImageBase.FiniteDiff` also provides a few finite difference operators, including `fdiff`,
  `fgradient`, etc.

## References

- [1] Rivest, Jean-Francois, Pierre Soille, and Serge Beucher. "Morphological gradients."
  Journal of Electronic Imaging 2.4 (1993): 326-336.
"""
function mgradient(img; dims=coords_spatial(img), r=nothing, mode=:beucher)
    return mgradient(img, strel_box(img, dims; r); mode)
end
function mgradient(img::AbstractArray{T}, se; mode=:beucher) where {T}
    out = similar(img, maybe_floattype(T))
    buffer = _make_gradient_buffer(img, mode)
    return mgradient!(out, img, se, buffer; mode)
end

"""
    mgradient!(out, img, buffer; [dims], [r], [mode])
    mgradient!(out, img, se, buffer; [mode])

The in-place version of [`mgradient`](@ref) with input image `img` and output image `out`.

The `buffer` array is required for `:beucher` mode. For `:internal` and `:external` modes,
`buffer` is not needed and can be `nothing`.
"""
function mgradient!(out, img, buffer; dims=coords_spatial(img), r=nothing, mode=:beucher)
    return mgradient!(out, img, strel_box(img, dims; r), buffer; mode)
end

function mgradient!(out, img, se, @nospecialize(buffer); mode=:beucher)
    require_symmetric_strel(se)
    @debug "calculate mgradient using $mode mode"
    if mode == :beucher
        isnothing(buffer) && throw(ArgumentError("buffer array is required for mode :beucher"))
        _beucher_gradient!(out, img, se, buffer)
    elseif mode == :internal
        _internal_gradient!(out, img, se)
    elseif mode == :external
        _external_gradient!(out, img, se)
    end
    return out
end
mgradient!(out::AbstractArray{<:Color3}, img, se, buffer; kwargs...) = throw(ArgumentError("color image is not supported"))

_make_gradient_buffer(img, mode) = mode == :beucher ? similar(img, maybe_floattype(eltype(img))) : nothing

function _internal_gradient!(out, img, se)
    erode!(out, img, se)
    @. out = img - out
    return out
end

function _external_gradient!(out, img, se)
    dilate!(out, img, se)
    @. out = out - img
    return out
end

function _beucher_gradient!(out, img, se, buffer)
    dilate!(out, img, se)
    erode!(buffer, img, se)
    @. out = out - buffer
    return out
end