text:
#####
##### mixtures of distributions
#####

"""
$(TYPEDEF)

Mixture of two distributions.

# Fields

$(FIELDS)
"""
struct Mix{T,L1,L2} <: SamplingLogDensity
    "Mixture weight on `ℓ1` (an object that supports [`weight_and_gradient`](@ref)."
    α::T
    "The log density with weight `α`."
    ℓ1::L1
    "The log density with weight `1-α`."
    ℓ2::L2
end

"""
$(SIGNATURES)

Return `α(x)` and `∇α(x)`, where `α` is the mixture weight on the first distribution.
"""
weight_and_gradient(α::Real, x) = α, zero(x)

"""
$(SIGNATURES)

Return `α(x)`, where `α` is the mixture weight on the first distribution. When not defined,
falls back to [`weight_and_gradient`](@ref).
"""
weight(α, x) = first(weight_and_gradient(α, x))

"""
$(SIGNATURES)

A *mixture* of two densities: `ℓ1(x)` with probability `α(x)`, and `ℓ2(x)` with probability
`1-α(x)`. `α` needs to implement [`weight_and_gradient`](@ref).
"""
function mix(α, ℓ1, ℓ2)
    @argcheck dimension(ℓ1) == dimension(ℓ2)
    Mix(α, ℓ1, ℓ2)
end

hypercube_dimension(ℓ::Mix) = max(hypercube_dimension(ℓ.ℓ1), hypercube_dimension(ℓ.ℓ2)) + 1

function hypercube_transform(ℓ::Mix, x)
    @unpack α, ℓ1, ℓ2 = ℓ
    if x[end] < weight(α, x[1:dimension(ℓ)])
        hypercube_transform(ℓ1, @view x[1:hypercube_dimension(ℓ1)])
    else
        hypercube_transform(ℓ2, @view x[1:hypercube_dimension(ℓ2)])
    end
end

dimension(ℓ::Mix) = dimension(ℓ.ℓ1)

function logdensity_and_gradient(ℓ::Mix, x)
    @unpack α, ℓ1, ℓ2 = ℓ
    ℓ1x, ∇ℓ1x = logdensity_and_gradient(ℓ1, x)
    ℓ2x, ∇ℓ2x = logdensity_and_gradient(ℓ2, x)
    αx, ∇αx = weight_and_gradient(α, x)
    ℓx = logaddexp(log(αx) + ℓ1x, log1p(-αx) + ℓ2x)
    f(x, y) = iszero(x) ? zero(x + y) : x * y # avoiding NaN in corner cases
    ∇ℓx = @. f((∇αx + αx * ∇ℓ1x), exp(ℓ1x - ℓx)) + f((-∇αx + (1 - αx) * ∇ℓ2x), exp(ℓ2x - ℓx))
    ℓx, ∇ℓx
end

struct DirectionalWeight{V}
    "The direction vector."
    y::V
end

"""
$(SIGNATURES)

Represents the weight `α(x) = logistic(dot(x, y))`. The Euclidean norm of `y` is related to
the slope of the gradient of `α` in the given direction.
"""
function directional_weight(y)
    @argcheck norm(y, 2) > 0
    DirectionalWeight(y)
end

function weight_and_gradient(α::DirectionalWeight, x)
    @unpack y = α
    αx = logistic(dot(x, y))
    ∇αx = (αx * (1 - αx)) .* y
    αx, ∇αx
end


text:
create(c::Client, ::Type{Message}, ch::DiscordChannel; kwargs...) = create_message(c, ch.id; kwargs...)

retrieve(c::Client, ::Type{Message}, ch::DiscordChannel, message::Integer) = get_channel_message(c, ch.id, message)
retrieve(
    c::Client,
    ::Type{Message},
    ch::DiscordChannel;
    pinned::Bool=false,
    kwargs...,
) = if pinned
        get_pinned_messages(c, ch.id)
    else
        get_channel_messages(c, ch.id; kwargs...)
    end

update(c::Client, m::Message; kwargs...) = edit_message(c, m.channel_id, m.id; kwargs...)

delete(c::Client, m::Message; pinned::Bool=false) = if pinned
        delete_pinned_channel_message(c, m.channel_id, m.id)
    else
        delete_message(c, m.channel_id, m.id)
    end
delete(c::Client, ms::Vector{Message}) = bulk_delete_messages(c, ms[1].channel_id; messages = map(m -> m.id, ms))

ext: jl


text:
using Zygote: pullback
using ProximalOperators

function ProximalOperators.gradient(f, x)
    fx, pb = pullback(f, x)
    grad = pb(one(fx))[1]
    return grad, fx
end

function ProximalOperators.gradient!(grad, f, x)
    y, fx = gradient(f, x)
    grad .= y
    return fx
end


text:
#!/usr/bin/env julia

"""
Tests for the [`Digits`](@ref) common module.
"""
module TestDigits

include(joinpath("..", "src", "Digits.jl"))

using Test

@testset "Digits" begin
    @testset "is_palindrome" begin
        @test Digits.is_palindrome(1)
        @test !Digits.is_palindrome(10)
        @test Digits.is_palindrome(11)
        @test Digits.is_palindrome(22)
        @test Digits.is_palindrome(UInt8(22))
        @test Digits.is_palindrome(Int16(22))
        @test Digits.is_palindrome(BigInt(22))
        @test Digits.is_palindrome(404)
        @test Digits.is_palindrome(18, base = 8)
        @test Digits.is_palindrome(24, base = 23)
        @test !Digits.is_palindrome(2, base = 2)
        @test !Digits.is_palindrome(146)
        @test Digits.is_palindrome(146, base = 8)
        @test !Digits.is_palindrome(1212)
        @test Digits.is_palindrome(2332)
        @test Digits.is_palindrome(72427)
        @test Digits.is_palindrome(0x72427, base = 16)
        @test Digits.is_palindrome(724427)
        @test Digits.is_palindrome(0x724427, base = 16)
        @test Digits.is_palindrome(722444227)
        @test Digits.is_palindrome(313, base = 2)
        @test !Digits.is_palindrome(513513)
        @test !Digits.is_palindrome(169803, base = 8)
        @test !Digits.is_palindrome(551133)
        @test !Digits.is_palindrome(0x551133, base = 16)
        @test !Digits.is_palindrome(5133150)
        @test !Digits.is_palindrome(5133150, base = 23)
    end
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
using FiniteDifferences

"""
line search wrapper
NOTE: The stepsize is defined as x - gamma * d

Returns the step size gamma and the Lipschitz estimate L
"""
function line_search_wrapper(
    line_search,
    t,
    f,
    grad!,
    x,
    d,
    gradient,
    dual_gap,
    L,
    gamma0,
    linesearch_tol,
    step_lim,
    gamma_max,
)
    if line_search isa Agnostic
        gamma = 2 // (2 + t)
    elseif line_search isa Goldenratio # FIX for general d
        gamma, _ = segment_search(
            f,
            grad!,
            x,
            d,
            gamma_max,
            linesearch_tol=linesearch_tol,
            inplace_gradient=true,
        )
    elseif line_search isa Backtracking # FIX for general d
        gamma, _ = backtrackingLS(
            f,
            gradient,
            x,
            d,
            gamma_max,
            linesearch_tol=linesearch_tol,
            step_lim=step_lim,
        )
    elseif line_search isa Nonconvex
        gamma = 1 / sqrt(t + 1)
    elseif line_search isa Shortstep
        gamma = min(max(fast_dot(gradient, d) * inv(L * norm(d)^2), 0), gamma_max)
    elseif line_search isa RationalShortstep
        gamma = min(max(fast_dot(gradient, d) * inv(L * fast_dot(d, d)), 0), gamma_max)
    elseif line_search isa FixedStep
        gamma = min(gamma0, gamma_max)
    elseif line_search isa Adaptive
        gamma, L = adaptive_step_size(f, grad!, gradient, x, d, L, gamma_max=gamma_max)
    end
    return gamma, L
end

"""
Slight modification of
Adaptive Step Size strategy from https://arxiv.org/pdf/1806.05123.pdf

Note: direction is opposite to the improving direction
norm(gradient, direction) > 0
TODO:
- make emphasis aware and optimize
"""
function adaptive_step_size(
    f,
    grad!,
    gradient,
    x,
    direction,
    L_est;
    eta=0.9,
    tau=2,
    gamma_max=1,
    upgrade_accuracy=false,
)
    #If there is no initial smoothness estimate
    #try to build one from the definition.
    if norm(direction) > 0
        if isnothing(L_est) || !isfinite(L_est)
            epsilon_step = min(1.0e-3, gamma_max)
            gradient_stepsize_estimation = similar(gradient)
            grad!(gradient_stepsize_estimation, x - epsilon_step * direction)
            L_est = norm(gradient - gradient_stepsize_estimation) / (epsilon_step * norm(direction))
        end
        M = eta * L_est
        if !upgrade_accuracy
            dot_dir = fast_dot(gradient, direction)
            ndir2 = norm(direction)^2
        else
            direction = big.(direction)
            x = big.(x)
            dot_dir = fast_dot(big.(gradient), direction)
            ndir2 = norm(direction)^2
        end

        gamma = min(max(dot_dir / (M * ndir2), 0.0), gamma_max)
        while f(x - gamma * direction) - f(x) > -gamma * dot_dir + gamma^2 * ndir2 * M / 2
            M *= tau
            gamma = min(max(dot_dir / (M * ndir2), 0.0), gamma_max)
        end
        return gamma, M
    else
        return 0, L_est
    end
end

# simple backtracking line search (not optimized)
# TODO:
# - code needs optimization

function backtrackingLS(
    f,
    grad_direction,
    x,
    d,
    gamma_max;
    line_search=true,
    linesearch_tol=1e-10,
    step_lim=20,
    lsTau=0.5,
)
    gamma = gamma_max * one(lsTau)
    i = 0

    dot_gdir = fast_dot(grad_direction, d)
    if dot_gdir ≤ 0
        @warn "Non-improving"
        return zero(gamma), i
    end

    oldVal = f(x)
    newVal = f(x - gamma * d)
    while newVal - oldVal > -linesearch_tol * gamma * dot_gdir
        if i > step_lim
            if oldVal - newVal >= 0
                return gamma, i
            else
                return zero(gamma), i
            end
        end
        gamma *= lsTau
        newVal = f(x - gamma * d)
        i += 1
    end
    return gamma, i
end

# simple golden-ratio based line search (not optimized)
# based on boostedFW paper code and adapted for julia
# TODO:
# - code needs optimization.
# In particular, passing a gradient container instead of allocating

function segment_search(
    f,
    grad,
    x,
    d,
    gamma_max;
    line_search=true,
    linesearch_tol=1e-10,
    inplace_gradient=true,
)
    # restrict segment of search to [x, y]
    y = x - gamma_max * d
    left, right = copy(x), copy(y)

    if inplace_gradient
        gradient = similar(d)
        grad(gradient, x)
        dgx = fast_dot(d, gradient)
        grad(gradient, y)
        dgy = fast_dot(d, gradient)
    else
        gradient = grad(x)
        dgx = fast_dot(d, gradient)
        gradient = grad(y)
        dgy = fast_dot(d, gradient)
    end

    # if the minimum is at an endpoint
    if dgx * dgy >= 0
        if f(y) <= f(x)
            return one(eltype(d)), y
        else
            return zero(eltype(d)), x
        end
    end

    # apply golden-section method to segment
    gold = (1 + sqrt(5)) / 2
    improv = Inf
    while improv > linesearch_tol
        old_left, old_right = left, right
        new = left + (right - left) / (1 + gold)
        probe = new + (right - new) / 2
        if f(probe) <= f(new)
            left, right = new, right
        else
            left, right = left, probe
        end
        improv = norm(f(right) - f(old_right)) + norm(f(left) - f(old_left))
    end

    x_min = (left + right) / 2

    # compute step size gamma
    gamma = zero(eltype(d))
    if line_search
        for i in eachindex(d)
            if d[i] != 0
                gamma = (x[i] - x_min[i]) / d[i]
                break
            end
        end
    end

    return gamma, x_min
end

##############################
### emphasis macro
##############################

macro emphasis(Emphasis, ex)
    return esc(quote
        if $Emphasis === memory
            @. $ex
        else
            $ex
        end
    end)
end

##############################
### Visualization etc
##############################

"""
plot_results

Given a series of list, generate subplots.
list_data_y -> contains a list of a list of lists (where each list refers to a subplot, and a list of lists refers to the y-values of the series inside a subplot).
list_data_x -> contains a list of a list of lists (where each list refers to a subplot, and a list of lists refers to the x-values of the series inside a subplot).
So if we have one plot with two series, these might look like:
    list_data_y = [[[1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6]]]
    list_data_x = [[[1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6]]]

And if we have two plots, each with two series, these might look like:
    list_data_y = [[[1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6]], [[7, 8, 9, 10, 11, 12], [7, 8, 9, 10, 11, 12]]]
    list_data_x = [[[1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6]], [[7, 8, 9, 10, 11, 12], [7, 8, 9, 10, 11, 12]]]

list_label -> contains the labels for the series that will be plotted,
which has to have a length equal to the number of series that are being plotted:
    list_label = ["Series 1", "Series 2"]

list_axis_x -> contains the labels for the x-axis that will be plotted,
which has to have a length equal to the number of subplots:
    list_axis_x = ["x-axis plot 1", "x-axis plot 1"]

list_axis_y -> Same as list_axis_x but for the y-axis

xscalelog -> A list of values indicating the type of axes to use in each subplot,
must be equal to the number of subplots:
    xscalelog = [:log, :identity]

yscalelog -> Same as xscalelog but for the y-axis

"""
function plot_results(
    list_data_y,
    list_data_x,
    list_label,
    list_axis_x,
    list_axis_y;
    filename=nothing,
    xscalelog=nothing,
    yscalelog=nothing,
    legend_position=nothing,
    list_style=fill(:solid, length(list_label)),
    list_color=get_color_palette(:auto, plot_color(:white)),
    list_markers=[
        :circle,
        :rect,
        :utriangle,
        :diamond,
        :hexagon,
        :+,
        :x,
        :star5,
        :cross,
        :xcross,
        :dtriangle,
        :rtriangle,
        :ltriangle,
        :pentagon,
        :heptagon,
        :octagon,
        :star4,
        :star6,
        :star7,
        :star8,
        :vline,
        :hline,
    ],
    number_markers_per_line=10,
)
    line_width = 3.0
    marker_size = 5.0
    transparency_markers = 0.45
    font_size_axis = 12
    font_size_legend = 9
    gr()
    plt = nothing
    list_plots = Plots.Plot{Plots.GRBackend}[]
    #Plot an appropiate number of plots
    for i in 1:length(list_data_x)
        for j in 1:length(list_data_x[i])
            if isnothing(xscalelog)
                xscale = :identity
            else
                xscale = xscalelog[i]
            end
            if isnothing(yscalelog)
                yscale = :log
            else
                yscale = yscalelog[i]
            end
            if isnothing(legend_position)
                position_legend = :best
                legend_display = true
            else
                position_legend = legend_position[i]
                if isnothing(position_legend)
                    legend_display = false
                else
                    legend_display = true
                end
            end
            if j == 1
                if legend_display
                    plt = plot(
                        list_data_x[i][j],
                        list_data_y[i][j],
                        label="",
                        xaxis=xscale,
                        yaxis=yscale,
                        ylabel=list_axis_y[i],
                        xlabel=list_axis_x[i],
                        legend=position_legend,
                        yguidefontsize=font_size_axis,
                        xguidefontsize=font_size_axis,
                        legendfontsize=font_size_legend,
                        width=line_width,
                        linestyle=list_style[j],
                        color=list_color[j],
                        grid=true,
                    )
                else
                    plt = plot(
                        list_data_x[i][j],
                        list_data_y[i][j],
                        label="",
                        xaxis=xscale,
                        yaxis=yscale,
                        ylabel=list_axis_y[i],
                        xlabel=list_axis_x[i],
                        yguidefontsize=font_size_axis,
                        xguidefontsize=font_size_axis,
                        width=line_width,
                        linestyle=list_style[j],
                        color=list_color[j],
                        grid=true,
                    )
                end
            else
                if legend_display
                    plot!(
                        list_data_x[i][j],
                        list_data_y[i][j],
                        label="",
                        width=line_width,
                        linestyle=list_style[j],
                        color=list_color[j],
                        legend=position_legend,
                    )
                else
                    plot!(
                        list_data_x[i][j],
                        list_data_y[i][j],
                        label="",
                        width=line_width,
                        linestyle=list_style[j],
                        color=list_color[j],
                    )
                end
            end
            if xscale == :log
                indices =
                    round.(
                        Int,
                        10 .^ (range(
                            log10(1),
                            log10(length(list_data_x[i][j])),
                            length=number_markers_per_line,
                        )),
                    )
                scatter!(
                    list_data_x[i][j][indices],
                    list_data_y[i][j][indices],
                    markershape=list_markers[j],
                    markercolor=list_color[j],
                    markersize=marker_size,
                    markeralpha=transparency_markers,
                    label=list_label[j],
                    legend=position_legend,
                )
            else
                scatter!(
                    view(
                        list_data_x[i][j],
                        1:length(list_data_x[i][j])÷number_markers_per_line:length(
                            list_data_x[i][j],
                        ),
                    ),
                    view(
                        list_data_y[i][j],
                        1:length(list_data_y[i][j])÷number_markers_per_line:length(
                            list_data_y[i][j],
                        ),
                    ),
                    markershape=list_markers[j],
                    markercolor=list_color[j],
                    markersize=marker_size,
                    markeralpha=transparency_markers,
                    label=list_label[j],
                    legend=position_legend,
                )
            end
        end
        push!(list_plots, plt)
    end
    fp = plot(list_plots..., layout=length(list_plots))
    plot!(size=(600, 400))
    if filename !== nothing
        savefig(fp, filename)
    end
    return fp
end

function plot_trajectories(
    data,
    label;
    filename=nothing,
    xscalelog=false,
    legend_position=:topright,
    lstyle=fill(:solid, length(data)),
)
    # theme(:dark)
    # theme(:vibrant)
    gr()

    x = []
    y = []
    pit = nothing
    pti = nothing
    dit = nothing
    dti = nothing
    offset = 2
    xscale = xscalelog ? :log : :identity
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][1] for j in offset:length(trajectory)]
        y = [trajectory[j][2] for j in offset:length(trajectory)]
        if i == 1
            pit = plot(
                x,
                y,
                label=label[i],
                xaxis=xscale,
                yaxis=:log,
                ylabel="Primal",
                legend=legend_position,
                yguidefontsize=8,
                xguidefontsize=8,
                legendfontsize=8,
                width=1.3,
                linestyle=lstyle[i],
            )
        else
            plot!(x, y, label=label[i], width=1.3, linestyle=lstyle[i])
        end
    end
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][5] for j in offset:length(trajectory)]
        y = [trajectory[j][2] for j in offset:length(trajectory)]
        if i == 1
            pti = plot(
                x,
                y,
                label=label[i],
                legend=false,
                xaxis=xscale,
                yaxis=:log,
                yguidefontsize=8,
                xguidefontsize=8,
                width=1.3,
                linestyle=lstyle[i],
            )
        else
            plot!(x, y, label=label[i], width=1.3, linestyle=lstyle[i])
        end
    end
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][1] for j in offset:length(trajectory)]
        y = [trajectory[j][4] for j in offset:length(trajectory)]
        if i == 1
            dit = plot(
                x,
                y,
                label=label[i],
                legend=false,
                xaxis=xscale,
                yaxis=:log,
                ylabel="Dual Gap",
                xlabel="Iterations",
                yguidefontsize=8,
                xguidefontsize=8,
                width=1.3,
                linestyle=lstyle[i],
            )
        else
            plot!(x, y, label=label[i], width=1.3, linestyle=lstyle[i])
        end
    end
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][5] for j in offset:length(trajectory)]
        y = [trajectory[j][4] for j in offset:length(trajectory)]
        if i == 1
            dti = plot(
                x,
                y,
                label=label[i],
                legend=false,
                xaxis=xscale,
                yaxis=:log,
                xlabel="Time",
                yguidefontsize=8,
                xguidefontsize=8,
                width=1.3,
                linestyle=lstyle[i],
            )
        else
            plot!(x, y, label=label[i], width=1.3, linestyle=lstyle[i])
        end
    end
    fp = plot(pit, pti, dit, dti, layout=(2, 2)) # layout = @layout([A{0.01h}; [B C; D E]]))
    plot!(size=(600, 400))
    if filename !== nothing
        savefig(fp, filename)
    end
    return fp
end

function plot_sparsity(data, label; filename=nothing, xscalelog=false, legend_position=:topright)
    # theme(:dark)
    # theme(:vibrant)
    gr()

    x = []
    y = []
    ps = nothing
    ds = nothing
    offset = 2
    xscale = xscalelog ? :log : :identity
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][6] for j in offset:length(trajectory)]
        y = [trajectory[j][2] for j in offset:length(trajectory)]
        if i == 1
            ps = plot(
                x,
                y,
                label=label[i],
                xaxis=xscale,
                yaxis=:log,
                ylabel="Primal",
                legend=legend_position,
                yguidefontsize=8,
                xguidefontsize=8,
                legendfontsize=8,
            )
        else
            plot!(x, y, label=label[i])
        end
    end
    for i in 1:length(data)
        trajectory = data[i]
        x = [trajectory[j][6] for j in offset:length(trajectory)]
        y = [trajectory[j][4] for j in offset:length(trajectory)]
        if i == 1
            ds = plot(
                x,
                y,
                label=label[i],
                legend=false,
                xaxis=xscale,
                yaxis=:log,
                ylabel="Dual",
                yguidefontsize=8,
                xguidefontsize=8,
            )
        else
            plot!(x, y, label=label[i])
        end
    end

    fp = plot(ps, ds, layout=(1, 2)) # layout = @layout([A{0.01h}; [B C; D E]]))
    plot!(size=(600, 200))
    if filename !== nothing
        savefig(fp, filename)
    end
    return fp
end

##############################################################
# simple benchmark of elementary costs of oracles and
# critical components
##############################################################

# TODO: add actual use of T for the rand(n)

function benchmark_oracles(f, grad!, x_gen, lmo; k=100, nocache=true)
    x = x_gen()
    sv = sizeof(x) / 1024^2
    println("\nSize of single atom ($(eltype(x))): $sv MB\n")
    to = TimerOutput()
    @showprogress 1 "Testing f... " for i in 1:k
        x = x_gen()
        @timeit to "f" temp = f(x)
    end
    @showprogress 1 "Testing grad... " for i in 1:k
        x = x_gen()
        temp = similar(x)
        @timeit to "grad" grad!(temp, x)
    end
    @showprogress 1 "Testing lmo... " for i in 1:k
        x = x_gen()
        @timeit to "lmo" temp = compute_extreme_point(lmo, x)
    end
    @showprogress 1 "Testing dual gap... " for i in 1:k
        x = x_gen()
        gradient = similar(x)
        grad!(gradient, x)
        v = compute_extreme_point(lmo, gradient)
        @timeit to "dual gap" begin
            dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
        end
    end
    @showprogress 1 "Testing update... (Emphasis: blas) " for i in 1:k
        x = x_gen()
        gradient = similar(x)
        grad!(gradient, x)
        v = compute_extreme_point(lmo, gradient)
        gamma = 1 / 2
        @timeit to "update (blas)" @emphasis(blas, x = (1 - gamma) * x + gamma * v)
    end
    @showprogress 1 "Testing update... (Emphasis: memory) " for i in 1:k
        x = x_gen()
        gradient = similar(x)
        grad!(gradient, x)
        v = compute_extreme_point(lmo, gradient)
        gamma = 1 / 2
        # TODO: to be updated to broadcast version once data structure ScaledHotVector allows for it
        @timeit to "update (memory)" @emphasis(memory, x = (1 - gamma) * x + gamma * v)
    end
    if !nocache
        @showprogress 1 "Testing caching 100 points... " for i in 1:k
            @timeit to "caching 100 points" begin
                cache = [gen_x() for _ in 1:100]
                x = gen_x()
                gradient = similar(x)
                grad!(gradient, x)
                v = compute_extreme_point(lmo, gradient)
                gamma = 1 / 2
                test = (x -> fast_dot(x, gradient)).(cache)
                v = cache[argmin(test)]
                val = v in cache
            end
        end
    end
    print_timer(to)
    return nothing
end

"""
`isequal` without the checks. Assumes a and b have the same axes.
"""
function _unsafe_equal(a::AbstractArray, b::AbstractArray)
    if a === b
        return true
    end
    @inbounds for idx in eachindex(a)
        if a[idx] != b[idx]
            return false
        end
    end
    return true
end

_unsafe_equal(a, b) = isequal(a, b)

fast_dot(A, B) = dot(A, B)

fast_dot(B::SparseArrays.SparseMatrixCSC, A::Matrix) = conj(fast_dot(A, B))

function fast_dot(A::Matrix{T1}, B::SparseArrays.SparseMatrixCSC{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    (m, n) = size(A)
    if (m, n) != size(B)
        throw(DimensionMismatch("Size mismatch"))
    end
    s = zero(T)
    if m * n == 0
        return s
    end
    rows = SparseArrays.rowvals(B)
    vals = SparseArrays.nonzeros(B)
    @inbounds for j in 1:n
        for ridx in SparseArrays.nzrange(B, j)
            i = rows[ridx]
            v = vals[ridx]
            s += v * conj(A[i, j])
        end
    end
    return s
end

"""
Check if the gradient using finite differences matches the grad! provided.
"""
function check_gradients(grad!, f, gradient, num_tests=10, tolerance=1.0e-5)
    for i in 1:num_tests
        random_point = rand(length(gradient))
        grad!(gradient, random_point)
        if norm(grad(central_fdm(5, 1), f, random_point)[1] - gradient) > tolerance
            @warn "There is a noticeable difference between the gradient provided and
            the gradient computed using finite differences."
        end
    end
end

"""
    trajectory_callback(storage)

Callback pushing the state at each iteration to the passed storage.
The state data is only the 5 first fields, usually:
`(t,primal,dual,dual_gap,time)`
"""
function trajectory_callback(storage)
    return function push_trajectory!(data)
        return push!(storage, Tuple(data)[1:5])
    end
end

# default printing for LineSearchMethod is just showing the type
Base.print(io::IO, ls::LineSearchMethod) = print(io, split(string(typeof(ls)), ".")[end])

"""
    MonotonousStepSize{F}

Represents a monotonous open-loop step size.
Contains a halving factor `N` increased at each iteration until there is primal progress
`gamma = 2 / (t + 2) * 2^(-N)`
"""
mutable struct MonotonousStepSize{F} <: LineSearchMethod
    domain_oracle::F
    factor::Int
end

MonotonousStepSize(f::F) where {F<:Function} = MonotonousStepSize{F}(f, 0)
MonotonousStepSize() = MonotonousStepSize(x -> true)

Base.print(io::IO, ::MonotonousStepSize) = print(io, "MonotonousStepSize")

function line_search_wrapper(
    line_search::MonotonousStepSize,
    t,
    f,
    grad!,
    x,
    d,
    gradient,
    dual_gap,
    L,
    gamma0,
    linesearch_tol,
    step_lim,
    gamma_max,
)
    gamma = 2.0^(1 - line_search.factor) / (2 + t)
    xnew = x - gamma * d
    f0 = f(x)
    while !line_search.domain_oracle(xnew) || f(xnew) > f0
        line_search.factor += 1
        gamma = 2.0^(1 - line_search.factor) / (2 + t)
        @. xnew = x - gamma * d
    end
    return gamma, L
end

"""
    MonotonousNonConvexStepSize{F}

Represents a monotonous open-loop non-convex step size.
Contains a halving factor `N` increased at each iteration until there is primal progress
`gamma = 1 / sqrt(t + 1) * 2^(-N)`
"""
mutable struct MonotonousNonConvexStepSize{F} <: LineSearchMethod
    domain_oracle::F
    factor::Int
end

MonotonousNonConvexStepSize(f::F) where {F<:Function} = MonotonousNonConvexStepSize{F}(f, 0)
MonotonousNonConvexStepSize() = MonotonousNonConvexStepSize(x -> true)

Base.print(io::IO, ::MonotonousNonConvexStepSize) = print(io, "MonotonousNonConvexStepSize")

function line_search_wrapper(
    line_search::MonotonousNonConvexStepSize,
    t,
    f,
    grad!,
    x,
    d,
    gradient,
    dual_gap,
    L,
    gamma0,
    linesearch_tol,
    step_lim,
    gamma_max,
)
    gamma = 2.0^(-line_search.factor) / sqrt(1 + t)
    xnew = x - gamma * d
    f0 = f(x)
    while !line_search.domain_oracle(xnew) || f(xnew) > f0
        line_search.factor += 1
        gamma = 2.0^(-line_search.factor) / sqrt(1 + t)
        @. xnew = x - gamma * d
    end
    return gamma, L
end

function print_callback(data, format_string; print_header=false, print_footer=false)
    print_formatted(fmt, args...) = @eval @printf($fmt, $(args...))
    if print_header || print_footer
        temp = strip(format_string, ['\n'])
        temp = replace(temp, "%" => "")
        temp = replace(temp, "e" => "")
        temp = replace(temp, "i" => "")
        temp = replace(temp, "s" => "")
        temp = split(temp, " ")
        len = 0
        for i in temp
            len = len + parse(Int, i)
        end
        lenHeaderFooter = len + 2 + length(temp) - 1
        if print_footer
            line = "-"^lenHeaderFooter
            @printf("%s\n\n", line)
        end
        if print_header
            line = "-"^lenHeaderFooter
            @printf("\n%s\n", line)
            s_format_string = replace(format_string, "e" => "s")
            s_format_string = replace(s_format_string, "i" => "s")
            print_formatted(s_format_string, data...)
            @printf("%s\n", line)
        end
    else
        print_formatted(format_string, data...)
    end
end

"""
    momentum_iterate(iter::MomentumIterator) -> ρ

Method to implement for a type `MomentumIterator`.
Returns the next momentum value `ρ` and updates the iterator internal state.
"""
function momentum_iterate end

"""
    ExpMomentumIterator{T}

Iterator for the momentum used in the variant of Stochastic Frank-Wolfe.
Momentum coefficients are the values of the iterator:
`ρ_t = 1 - num / (offset + t)^exp`

The state corresponds to the iteration count.

Source:
Stochastic Conditional Gradient Methods: From Convex Minimization to Submodular Maximization
Aryan Mokhtari, Hamed Hassani, Amin Karbasi, JMLR 2020.
"""
mutable struct ExpMomentumIterator{T}
    exp::T
    num::T
    offset::T
    iter::Int
end

ExpMomentumIterator() = ExpMomentumIterator(2/3, 4.0, 8.0, 0)

function momentum_iterate(em::ExpMomentumIterator)
    em.iter += 1
    return 1 - em.num / (em.offset + em.iter)^(em.exp)
end

"""
    ConstantMomentumIterator{T}

Iterator for momentum with a fixed damping value, always return the value and a dummy state.
"""
struct ConstantMomentumIterator{T}
    v::T
end

momentum_iterate(em::ConstantMomentumIterator) = em.v

# batch sizes

"""
    batchsize_iterate(iter::BatchSizeIterator) -> b

Method to implement for a batch size iterator of type `BatchSizeIterator`.
Calling `batchsize_iterate` returns the next batch size and typically update the internal state of `iter`.
"""
function batchsize_iterate end

"""
    ConstantBatchIterator(batch_size)

Batch iterator always returning a constant batch size.
"""
struct ConstantBatchIterator
    batch_size::Int
end

batchsize_iterate(cbi::ConstantBatchIterator) = cbi.batch_size

"""
    IncrementBatchIterator(starting_batch_size, max_batch_size, [increment = 1])

Batch size starting at starting_batch_size and incrementing by `increment` at every iteration.
"""
mutable struct IncrementBatchIterator
    starting_batch_size::Int
    max_batch_size::Int
    increment::Int
    iter::Int
    maxreached::Bool
end

function IncrementBatchIterator(starting_batch_size::Int, max_batch_size::Int, increment::Int)
    return IncrementBatchIterator(starting_batch_size, max_batch_size, increment, 0, false)
end

function IncrementBatchIterator(starting_batch_size::Int, max_batch_size::Int)
    return IncrementBatchIterator(starting_batch_size, max_batch_size, 1, 0, false)
end

function batchsize_iterate(ibi::IncrementBatchIterator)
    if ibi.maxreached
        return ibi.max_batch_size
    end
    new_size = ibi.starting_batch_size + ibi.iter * ibi.increment
    ibi.iter += 1
    if new_size > ibi.max_batch_size
        ibi.maxreached = true
        return ibi.max_batch_size
    end
    return new_size
end


text:
function pipefix(numbers::Vector{Int64})::Vector{Int64}
  min,max = extrema(numbers)
  range(min, stop=max)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Overload for writing.
function moi_to_object(
    foo::Nonlinear, name_map::Dict{MOI.VariableIndex, String}
)
    node_list = OrderedObject[]
    foo_object = convert_expr_to_mof(foo.expr, node_list, name_map)
    return OrderedObject(
        "head" => "ScalarNonlinearFunction",
        "root" => foo_object,
        "node_list" => node_list,
    )
end

# Overload for reading.
function function_to_moi(
    ::Val{:ScalarNonlinearFunction},
    object::T,
    name_map::Dict{String, MOI.VariableIndex}
) where {T <: Object}
    node_list = T.(object["node_list"])
    expr = convert_mof_to_expr(object["root"], node_list, name_map)
    return Nonlinear(expr)
end

function lift_variable_indices(expr::Expr)
    if expr.head == :ref && length(expr.args) == 2 && expr.args[1] == :x
        return expr.args[2]
    else
        for (index, arg) in enumerate(expr.args)
            expr.args[index] = lift_variable_indices(arg)
        end
    end
    return expr
end

lift_variable_indices(arg) = arg  # Recursion fallback.

function extract_function_and_set(expr::Expr)
    if expr.head == :call  # One-sided constraint or foo-in-set.
        @assert length(expr.args) == 3
        if expr.args[1] == :in
            # return expr.args[2], expr.args[3]
            error("Constraints of the form foo-in-set aren't supported yet.")
        elseif expr.args[1] == :(<=)
            return expr.args[2], MOI.LessThan(expr.args[3])
        elseif expr.args[1] == :(>=)
            return expr.args[2], MOI.GreaterThan(expr.args[3])
        elseif expr.args[1] == :(==)
            return expr.args[2], MOI.EqualTo(expr.args[3])
        end
    elseif expr.head == :comparison  # Two-sided constraint.
        @assert length(expr.args) == 5
        if expr.args[2] == expr.args[4] == :(<=)
            return expr.args[3], MOI.Interval(expr.args[1], expr.args[5])
        elseif expr.args[2] == expr.args[4] == :(>=)
            return expr.args[3], MOI.Interval(expr.args[5], expr.args[1])
        end
    end
    error("Oops. The constraint $(expr) wasn't recognised.")
end

function write_nlpblock(
    object::T, model::Model, name_map::Dict{MOI.VariableIndex, String}
) where {T <: Object}
    # TODO(odow): is there a better way of checking if the NLPBlock is set?
    nlp_block = try
        MOI.get(model, MOI.NLPBlock())
    catch ex
        @assert isa(ex, KeyError)
        return  # No NLPBlock set.
    end
    MOI.initialize(nlp_block.evaluator, [:ExprGraph])
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    if nlp_block.has_objective
        objective = MOI.objective_expr(nlp_block.evaluator)
        objective = lift_variable_indices(objective)
        sense = MOI.get(model, MOI.ObjectiveSense())
        object["objective"] = T(
            "sense" => moi_to_object(sense),
            "function" => moi_to_object(Nonlinear(objective), name_map)
        )
    end
    for (row, bounds) in enumerate(nlp_block.constraint_bounds)
        constraint = MOI.constraint_expr(nlp_block.evaluator, row)
        (func, set) = extract_function_and_set(constraint)
        func = lift_variable_indices(func)
        push!(
            object["constraints"],
            T(
                "function" => moi_to_object(Nonlinear(func), name_map),
                "set"      => moi_to_object(set, name_map)
            )
        )
    end
end

#=
Expr:
    2 * x + sin(x)^2 + y

Tree form:
                    +-- (2)
          +-- (*) --+
          |         +-- (x)
    (+) --+
          |         +-- (sin) --+-- (x)
          +-- (^) --+
          |         +-- (2)
          +-- (y)

MOF format:

    {
        "head": "nonlinear",
        "root": {"head": "node", "index": 4},
        "node_list": [
            {
                "head": "*", "args": [
                    {"head": "real", "value": 2},
                    {"head": "variable", "name": "x"}
                ]
            },
            {
                "head": "sin",
                "args": [
                    {"head": "variable", "name", "x"}
                ]
            }
            {
                "head": "^",
                "args": [
                    {"head": "node", "index": 2},
                    {"head": "real", "value": 2}
                ]
            },
            {
                "head": "+",
                "args": [
                    {"head": "node", "index": 1},
                    {"head": "node", "index": 3},
                    {"head": "variable", "name": "y"}
                ]
            }
        ]
    }
=#

"""
    ARITY

The arity of a nonlinear function. One of:
 - `Nary` if the function accepts one or more arguments
 - `Unary` if the function accepts exactly one argument
 - `Binary` if the function accepts exactly two arguments
 - `Ternary` if the function accepts exactly three arguments.
"""
@enum ARITY Nary Unary Binary Ternary

# A nice error message telling the user they supplied the wrong number of
# arguments to a nonlinear function.
function validate_arguments(function_name, arity::ARITY, num_arguments::Int)
    if ((arity == Nary && num_arguments < 1) ||
        (arity == Unary && num_arguments != 1) ||
        (arity == Binary && num_arguments != 2) ||
        (arity == Ternary && num_arguments != 3))
        error("The function $(function_name) is a $(arity) function, but you " *
              "have passed $(num_arguments) arguments.")
    end
end

"""
    SUPPORTED_FUNCTIONS

A vector of string-symbol pairs that map the MathOptFormat string representation
(i.e, the value of the `"head"` field) to the name of a Julia function (in
Symbol form).
"""
const SUPPORTED_FUNCTIONS = Pair{String, Tuple{Symbol, ARITY}}[
    # ==========================================================================
    # The standard arithmetic functions.
    # The addition operator: +(a, b, c, ...) = a + b + c + ...
    # In the unary case, +(a) = a.
    "+"     => (:+, Nary),
    # The subtraction operator: -(a, b, c, ...) = a - b - c - ...
    # In the unary case, -(a) = -a.
    "-"     => (:-, Nary),
    # The multiplication operator: *(a, b, c, ...) = a * b * c * ...
    # In the unary case, *(a) = a.
    "*"     => (:*, Nary),
    # The division operator. This must have exactly two arguments. The first
    # argument is the numerator, the second argument is the denominator:
    # /(a, b) = a / b.
    "/"     => (:/, Binary),
    # ==========================================================================
    # N-ary minimum and maximum functions.
    "min"   => (:min, Nary),
    "max"   => (:max, Nary),
    # ==========================================================================
    # floor(x) = ⌊x⌋
    "floor"   => (:floor, Unary),
    # ceil(x) = ⌈x⌉
    "ceil"   => (:ceil, Unary),
    # ==========================================================================
    # The absolute value function: abs(x) = (x >= 0 ? x : -x).
    "abs"   => (:abs, Unary),
    # ==========================================================================
    # Log- and power-related functions.
    # A binary function for exponentiation: ^(a, b) = a ^ b.
    "^"     => (:^, Binary),
    # The natural exponential function: exp(x) = e^x.
    "exp"   => (:exp, Unary),
    # The base-e log function: y = log(x) => e^y = x.
    "log"   => (:log, Unary),
    # The base-10 log function: y = log10(x) => 10^y = x.
    "log10" => (:log10, Unary),
    # The square root function: sqrt(x) = √x = x^(0.5).
    "sqrt"  => (:sqrt, Unary),
    # ==========================================================================
    # Boolean operators
    # A && B = A and B
    # "&&" => (:&&, Binary),
    # A || B = A or B
    # "||" => (:||, Binary),
    # ==========================================================================
    # Comparison operators
    "<" => (:<, Binary),
    "<=" => (:<=, Binary),
    ">" => (:>, Binary),
    ">=" => (:>=, Binary),
    "==" => (:(==), Binary),
    "!=" => (:!=, Binary),
    # ==========================================================================
    "ifelse" => (:ifelse, Ternary),
    # "not" => (:!, Unary),
    # ==========================================================================
    # The unary trigonometric functions. These must have exactly one argument.
    "cos"   => (:cos, Unary),
    "cosh"  => (:cosh, Unary),
    "acos"  => (:acos, Unary),
    "acosh" => (:acosh, Unary),
    "sin"   => (:sin, Unary),
    "sinh"  => (:sinh, Unary),
    "asin"  => (:asin, Unary),
    "asinh" => (:asinh, Unary),
    "tan"   => (:tan, Unary),
    "tanh"  => (:tanh, Unary),
    "atan"  => (:atan, Unary),
    # "atan2" => (:atan, Binary),
    "atanh" => (:atanh, Unary)
]

# An internal helper dictionary that maps function names in Symbol form to their
# MathOptFormat string representation.
const FUNCTION_TO_STRING = Dict{Symbol, Tuple{String, ARITY}}()

# An internal helper dictionary that maps function names in their MathOptFormat
# string representation to the symbol representing the Julia function.
const STRING_TO_FUNCTION = Dict{String, Tuple{Symbol, ARITY}}()

# Programatically add the list of supported functions to the helper dictionaries
# for easy of look-up later.
for (mathoptformat_string, (julia_symbol, num_arguments)) in SUPPORTED_FUNCTIONS
    FUNCTION_TO_STRING[julia_symbol] = (mathoptformat_string, num_arguments)
    STRING_TO_FUNCTION[mathoptformat_string] = (julia_symbol, num_arguments)
end

"""
    convert_mof_to_expr(
        node::T, node_list::Vector{T}, name_map::Dict{String, MOI.VariableIndex}
    )

Convert a MathOptFormat node `node` into a Julia expression given a list of
MathOptFormat nodes in `node_list`. Variable names are mapped through `name_map`
to their variable index.
"""
function convert_mof_to_expr(
    node::T, node_list::Vector{T}, name_map::Dict{String, MOI.VariableIndex}
) where {T <: Object}
    head = node["head"]
    if head == "real"
        return node["value"]
    elseif head == "complex"
        return Complex(node["real"], node["imag"])
    elseif head == "variable"
        return name_map[node["name"]]
    elseif head == "node"
        return convert_mof_to_expr(node_list[node["index"]], node_list, name_map)
    else
        if !haskey(STRING_TO_FUNCTION, head)
            error("Cannot convert MOF to Expr. Unknown function: $(head).")
        end
        (julia_symbol, arity) = STRING_TO_FUNCTION[head]
        validate_arguments(head, arity, length(node["args"]))
        expr = Expr(:call, julia_symbol)
        for arg in node["args"]
            push!(expr.args, convert_mof_to_expr(arg, node_list, name_map))
        end
        return expr
    end
end

"""
    convert_expr_to_mof(
        node::T,
        node_list::Vector{T},
        name_map::Dict{MOI.VariableIndex, String}
    )

Convert a Julia expression into a MathOptFormat representation. Any intermediate
nodes that are required are appended to `node_list`. Variable indices are mapped
through `name_map` to their string name.
"""
function convert_expr_to_mof(
    expr::Expr,
    node_list::Vector{T},
    name_map::Dict{MOI.VariableIndex, String},
) where {T <: Object}
    if expr.head != :call
        error("Expected an expression that was a function. Got $(expr).")
    end
    function_name = expr.args[1]
    if !haskey(FUNCTION_TO_STRING, function_name)
        error("Cannot convert Expr to MOF. Unknown function: $(function_name).")
    end
    (mathoptformat_string, arity) = FUNCTION_TO_STRING[function_name]
    validate_arguments(function_name, arity, length(expr.args) - 1)
    node = T("head" => mathoptformat_string, "args" => T[])
    for arg in @view(expr.args[2:end])
        push!(node["args"], convert_expr_to_mof(arg, node_list, name_map))
    end
    push!(node_list, node)
    return T("head" => "node", "index" => length(node_list))
end

# Recursion end for variables.
function convert_expr_to_mof(
    variable::MOI.VariableIndex,
    ::Vector{T},
    name_map::Dict{MOI.VariableIndex, String},
) where {T <: Object}
    return T("head" => "variable", "name" => name_map[variable])
end

# Recursion end for real constants.
function convert_expr_to_mof(
    value::Real, ::Vector{T}, name_map::Dict{MOI.VariableIndex, String}
) where {T <: Object}
    return T("head" => "real", "value" => value)
end

# Recursion end for complex numbers.
function convert_expr_to_mof(
    value::Complex, ::Vector{T}, ::Dict{MOI.VariableIndex, String}
) where {T <: Object}
    return T("head" => "complex", "real" => real(value), "imag" => imag(value))
end

# Recursion fallback.
function convert_expr_to_mof(
    fallback, ::Vector{<:Object}, ::Dict{MOI.VariableIndex, String}
)
    error("Unexpected $(typeof(fallback)) encountered: $(fallback).")
end


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

"""
    S = compSet([C₁,C₂,…,Cₖ])
    S = compSet([ψ₁,ψ₂,…,ψₖ])

Create a *component set* 'compSet' by providing a
vector of *canonical triplets*.

# Examples
```@example
 using ISA
 𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
 𝐶₁ = AMFMtriplet(t->1.0,t->10*t,0.1)
 𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
 𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
```
Another way to create a *component set* 'compSet'
by providing a vector of numerical *AM--FM components*.

```@example
 using ISA
 𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
 ψ₀ = AMFMcomp(𝐶₀)
 𝐶₁ = AMFMtriplet(t->1.0,t->10*t,0.1)
 ψ₁ = AMFMcomp(𝐶₁)
 𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
 ψ₂ = AMFMcomp(𝐶₂)
 𝑆 = compSet([ψ₀,ψ₁,ψ₂])
```
"""
struct compSet
  S::Vector{AMFMtriplet}
end

# CONSTRUCTIONS
function compSet(V::Vector{AMFMcomp})
  return compSet([ψ.C  for ψ ∈ V])
end

# DISPLAY
Base.show(io::IO, x::compSet) = print(io, "component set")

"""
    z = AMFMmodel(S)
    z = AMFMmodel([ψ₁,ψ₂,…,ψₖ])
    z = AMFMmodel([C₁,C₂,…,Cₖ])

Create a 'AMFMmodel' paramtertized by a 'compSet'.

# Examples
```@example
using ISA
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->10*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
```
Another convenient way to create a 'AMFMmodel'
by providing a vector of numerical *AM--FM components*.

```@example
 using ISA
 𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
 ψ₀ = AMFMcomp(𝐶₀)
 𝐶₁ = AMFMtriplet(t->1.0,t->10*t,0.1)
 ψ₁ = AMFMcomp(𝐶₁)
 𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
 ψ₂ = AMFMcomp(𝐶₂)
 z = AMFMmodel([ψ₀,ψ₁,ψ₂])
```
Another convenient way to create a 'AMFMmodel'
by providing a vector of *canonical triplets*.

```@example
using ISA
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->10*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
z = AMFMmodel([𝐶₀,𝐶₁,𝐶₂])
```
"""
mutable struct AMFMmodel
  S::compSet
end

# CONSTRUCTIONS
function AMFMmodel(W::Vector{AMFMcomp})
        return AMFMmodel(compSet(W))
end
function AMFMmodel(X::Vector{AMFMtriplet})
  return  AMFMmodel(compSet(X))
end

# METHODS
function (z::AMFMmodel)(t::Vector{<:Real})
  out = zeros(size(t))
  for k∈1:length(z.S.S)
    ψ = AMFMcomp(z.S.S[k])
    out += ψ(t)
  end
  return out
end
function (z::AMFMmodel)(t::Real)
  return (z::AMFMmodel)([t])[1]
end
function (z::AMFMmodel)(t::StepRangeLen)
  return (z::AMFMmodel)(collect(t))
end
function (z::AMFMmodel)(t::UnitRange)
  return (z::AMFMmodel)(collect(t))
end

# DISPLAY
Base.show(io::IO, x::AMFMmodel) = print(io, "AM--FM model")


text:
# src path
<<<<<<< HEAD
src_path = "C:\\Users\\koike\\Desktop\\git\\general_2d_NS_LES\\src\\"

=======
src_path = "C:\\Users\\hidee\\Desktop\\git\\general_2d_NS_LES\\src_muscl\\"
>>>>>>> c15c8afc37995e56aa93b21bda87f4a9cb53eeb5

# main (変更しないこと)
src_read="read_grid.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="output.jl"
include(src_path*src_read)

src_read="advection_term.jl"
include(src_path*src_read)
src_read="boundary.jl"
include(src_path*src_read)
src_read="converge.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_read="muscl.jl"
include(src_path*src_read)
src_read="rhs.jl"
include(src_path*src_read)
src_read="value_setup.jl"
include(src_path*src_read)
src_read="turb_val.jl"
include(src_path*src_read)
src_read="viscos_pturb.jl"
include(src_path*src_read)

src_read="allocation.jl"
include(src_path*src_read)
src_read="cal_time_step.jl"
include(src_path*src_read)
src_main="main.jl"
include(src_path*src_main)


text:
# There are many 3 primary stages of the experiment where representation
# is important: at transfection, at the bottleneck(s), and at
# sequencing.

function scan_rep_space(filepath; debug=false, quiet=false)

    if !debug
        parameters = Dict{Symbol, Vector}(
            :representation => map(x->round(Int, x), logspace(0, 4, 13)),
            :bottleneck_representation => map(x->round(Int, x),  logspace(0,4,13)),
            :seq_depth => map(x->round(Int, x), logspace(0, 4, 13))
        )
        num_runs = 10
    else
        parameters = Dict{Symbol, Vector}(
            :representation => map(x->round(Int, x), logspace(0, 2, 2)),
            :bottleneck_representation => map(x->round(Int, x),  logspace(0,2,2)),
            :seq_depth => [10^2]
        )
        num_runs = 1
    end

    const overlap = intersect(fieldnames(FacsScreen), fieldnames(GrowthScreen))
    # custom function for handling both growth and a facs screen in the
    # same relational datastructure
    flatten_overlap = (setup, lib) -> begin
        local results = Any[typeof(setup)]
        for name in overlap
            push!(results, getfield(setup, name))
        end
        results
    end

    methods = [venn, auprc]
    measures = [:inc, :dec, :incdec]
    genetypes = [:sigmoidal, :linear, :all]
    test_method_wrapper = (bc_counts, genes) -> test_methods(genes, methods, measures, genetypes)

    results = []
    for screentype in [FacsScreen(), GrowthScreen()]
        for crisprtype in [CRISPRi(), CRISPRn()]
            if typeof(screentype) == GrowthScreen
                max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
                    :inactive => (0.83, Delta(0.0)),
                    :negcontrol => (0.05, Delta(0.0)),
                    :increasing => (0.02, TruncatedNormal(0.1, 0.1, 0.025, 1)),
                    :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
                )
                lib = Library(max_phenotype_dists, crisprtype)
            else
                lib = Library(crisprtype)
            end

            runs = build_parameter_space(screentype, parameters, num_runs)

            before = time()
            result = pmap(args -> run_exp(args[1], lib, test_method_wrapper; run_idx=args[2], flatten_func=flatten_overlap), runs)
            (!quiet) && println("$(time() - before) seconds")
            result = DataFrame(permutedims(hcat(result...), [2, 1]))
            result[:crisprtype] = typeof(crisprtype)
            push!(results, result)
        end
    end
    results = vcat(results...)
    hierarchy = vcat([hcat(item...) for item in IterTools.product(map(Symbol, methods), measures, genetypes)]...)
    new_names = [[:method, :measure, :genetype, :score, :screen]...; overlap...; :run; :crisprtype]
    results = construct_hierarchical_label(hierarchy, results, new_names)
    CSV.write(filepath, results)
end


text:
const water_fraction_to_M = 1.0m^3*m^-3 * 1kg*L^-1 / 18.0g*mol^-1

using Base: tail

abstract type AbstractLayerInterpolator end

struct LinearLayerInterpolator <: AbstractLayerInterpolator
    layer::Int
    frac::Float64
end

struct LinearLayerInterpolators{H} <: AbstractLayerInterpolator
    range::LinearLayerInterpolator
    increment::LinearLayerInterpolator
    height::H
end
LinearLayerInterpolators(env, height) =
    LinearLayerInterpolators(range_interpolator(env, height), increment_interpolator(env, height), height)

# Use recursion so the compiler does this with tuples. Huge performance improvement.
@inline layer_sizes(env::MicroclimPoint{S,I,NL,R}) where {S,I,NL,R} = layer_sizes(zero(NL), I, NL)
@inline layer_sizes(z, I::Tuple, NL) = (lsize(z, I[1], I[2]), layer_sizes(I, NL)...)
@inline layer_sizes(I::Tuple, NL) = (lsize(I[1], I[2], I[3]), layer_sizes(tail(I), NL)...)
@inline layer_sizes(I::Tuple{X,X}, NL) where X = (lsize(I[1], I[2], NL),)

@inline lsize(a, b, c) = (c + b)  / 2 - (b + a) / 2

@inline layer_bounds(env::MicroclimPoint{S,I,NL,R}) where {S,I,NL,R} = layer_bounds(I, NL)
@inline layer_bounds(inc::Tuple{X,Vararg}, next) where X = (bound(inc[1], inc[2]), layer_bounds(tail(inc), next)...)
@inline layer_bounds(inc::Tuple{X}, next) where X = (bound(inc[1], next),)

@inline bound(a::Number, b::Number) = (a + b) / 2

@inline max_height(env) = sum(layer_sizes(env))
@inline layer_props(env) = layer_sizes(env) ./ max_height(env)

" Calculate current interpolation layers and fraction from NicheMapR data"
@inline increment_interpolator(env, height) = begin
    increments = get_increments(env)
    for (i, upper_height) in enumerate(increments)
        if upper_height > height
            lower_height = increments[i - 1]
            frac = (height - lower_height)/(upper_height - lower_height)
            return LinearLayerInterpolator(i, frac)
        end
    end
    # Otherwise it's taller/deeper than we have data, so use the largest we have.
    LinearLayerInterpolator(lastindex(increments), 1.0)
end

@inline range_interpolator(env, height) = begin
    lower, upper = get_range(env)
    h = min(max(height, lower), upper)
    frac = (h - lower)/(upper - lower)
    LinearLayerInterpolator(2, frac)
end

" Interpolate between two layers of environmental data. "
@inline interp_layer(layers, interp, i) =
    lin_interp(layers, i, interp.layer - 1) * (oneunit(interp.frac) - interp.frac) +
    lin_interp(layers, i, interp.layer) * interp.frac

" Linear interpolation "
@inline lin_interp(vector::Vector, i) = begin
    int = floor(Int64, i)
    frac::Float64 = i - int
    vector[int+1] * (1 - frac) + vector[int] * frac
end
@inline lin_interp(matrix::Matrix, i, j) = begin
    int = floor(Int64, i)
    frac::Float64 = i - int
    matrix[int+1, j] * (1 - frac) + matrix[int, j] * frac
end
@inline lin_interp(vector::Vector, i::Int) = vector[i]
@inline lin_interp(matrix::Matrix, i::Int, j) = matrix[i, j]

@inline weightedmean(env, layers, inpterp::LinearLayerInterpolators, i) =
    weightedmean(env, layers, inpterp.height, i)
weightedmean(env, layers, height::Number, i) = begin
    wmean = zero(layers[1][1])*m
    lbounds = layer_bounds(env)
    lsizes = layer_sizes(env)
    h = min(height, max_height(env))
    # sum all layers < height, after size adjustment
    for l = 1:length(lbounds)
        if h > lbounds[l]
            # add the whole layer
            wmean += lin_interp(layers, i, l) * lsizes[l]
        else
            # add the last fractional layer and break the loop
            frac = 1 - ((lbounds[l] - h) / lsizes[l])
            wmean += lin_interp(layers, i, l) * lsizes[l] * frac
            break
        end
    end
    wmean / h
end

@inline layermax(layers, interp, i) = begin
    val = layers[1][i]
    for l = 2:interp.layer
        val = max(val, layers[l][i])
    end
    val
end


text:
using CounterfactualExplanations
using CounterfactualExplanations.Models
using Random
using LinearAlgebra
using Flux

@testset "Exceptions" begin
    @testset "LogisticModel" begin
        w = [1,2] # vector instead of matrix
        b = 0 # scalar instead of array
        @test_throws MethodError LogisticModel(w, b)
    end

    @testset "BayesianLogisticModel" begin

        μ = [0, 1.0, -2.0] # vector instead of matrix
        Σ = Symmetric(reshape(randn(9),3,3).*0.1 + UniformScaling(1.0)) # MAP covariance matrix
        @test_throws MethodError BayesianLogisticModel(μ, Σ)

        # Dimensions not matching:
        μ = [0 1.0]
        Σ = Symmetric(reshape(randn(9),3,3).*0.1 + UniformScaling(1.0))
        @test_throws DimensionMismatch BayesianLogisticModel(μ, Σ)

    end

    @testset "logits" begin
        M(x) = [2 1] * x # model not declared as subtype of AbstractFittedModel
        x = [1,1]
        @test_throws MethodError logits(M, x)
    end

    @testset "probs" begin
        M(x) = [2 1] * x # model not declared as subtype of AbstractFittedModel
        x = [1,1]
        @test_throws MethodError probs(M, x)
    end
end

@testset "Predictions" begin

    @testset "LogisticModel" begin
        M = LogisticModel([1 1],[0])
        x = [1,1]
        @test logits(M, x)[1] == 2
        @test probs(M, x)[1] == σ(2)
    end

    @testset "BayesianLogisticModel" begin

        # MLE:
        μ = [0 1.0 1.0] # vector instead of matrix
        Σ = zeros(3,3) # MAP covariance matrix
        M = BayesianLogisticModel(μ, Σ)
        x = [1,1]
        @test logits(M, x)[1] == 2
        @test probs(M, x)[1] == σ(2)

        # Not MLE:
        μ = [0 1.0 1.0] # vector instead of matrix
        Σ = zeros(3,3) + UniformScaling(1) # MAP covariance matrix
        M = BayesianLogisticModel(μ, Σ)
        x = [1,1]
        @test logits(M, x)[1] == 2
        @test probs(M, x)[1] != σ(2) # posterior predictive using probit link function

    end

end


text:
include("convective_bl_model.jl")
function main()

    # TODO: this will move to the future namelist functionality
    cbl_args = ArgParseSettings(autofix_names = true)
    add_arg_group!(cbl_args, "ConvectiveBoundaryLayer")
    @add_arg_table! cbl_args begin
        "--surface-flux"
        help = "specify surface flux for energy and moisture"
        metavar = "prescribed|bulk"
        arg_type = String
        default = "bulk"
    end

    cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = cbl_args)

    surface_flux = cl_args["surface_flux"]

    FT = Float64
    config_type = AtmosLESConfigType

    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(80)
    Δv = FT(80)

    resolution = (Δh, Δh, Δv)

    # Prescribe domain parameters
    xmax = FT(4800)
    ymax = FT(4800)
    zmax = FT(3200)

    t0 = FT(0)

    # Full simulation requires 16+ hours of simulated time
    timeend = FT(3600 * 0.1)
    CFLmax = FT(0.4)

    # Choose default Explicit solver
    ode_solver_type = ClimateMachine.ExplicitSolverType()

    model = convective_bl_model(FT, config_type, zmax, surface_flux)
    ics = model.problem.init_state_prognostic

    # Assemble configuration
    driver_config = ClimateMachine.AtmosLESConfiguration(
        "ConvectiveBoundaryLayer",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        ics,
        solver_type = ode_solver_type,
        model = model,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
        CFL_direction = HorizontalDirection(),
    )
    dgn_config = config_diagnostics(driver_config)

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            ("moisture.ρq_tot",),
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    check_cons = (
        ClimateMachine.ConservationCheck("ρ", "1mins", FT(0.0001)),
        ClimateMachine.ConservationCheck("ρe", "1mins", FT(0.0025)),
    )

    result = ClimateMachine.invoke!(
        solver_config;
        user_callbacks = (cbtmarfilter,),
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        check_euclidean_distance = true,
    )
end

main()


text:
module NuFlux

using DocStringExtensions
using DelimitedFiles
using StatsBase
using Statistics
using Corpuscles
using Interpolations
using LRUCache

abstract type Flux end

struct FluxTable <: Flux
    coszentihbinedges::Vector{Float64}
    azimuthbinedges::Vector{Float64}
    energies::Vector{Float64}
    particle::Particle
    flux::Array{Float64, 3}
end

struct FluxBinInfo
    CosZenithBin::Tuple{Rational, Rational}
    AzimutBin::Tuple{Rational, Rational}
end

FluxBinInfo(cosZmin, cosZmax, azimutmin, azimutmax) = FluxBinInfo((cosZmin, cosZmax), (azimutmin, azimutmax))

function parsefluxbininfo(line)
    regex = eachmatch(r"-?\d+\.?\d*", line)
    numbers = collect(rationalize(parse(Float64, m.match)) for m in regex)
    FluxBinInfo(numbers...)
end

function _getbinmids(binedges)
    Float64.(mean.(binedges[2:end] .+ binedges[1:end-1]) ./ 2)
end

function _makerange(x)
    diffs = x[2:end] - x[1:end-1]
    std(diffs) > 1e-4 && error("value spacing not constant")
    minimum(x):mean(diffs):maximum(x)
end

function readfluxfile(io)
    lines = readlines(io)
    idx = findall(l->occursin("average",l), lines)
    fluxes = Dict()
    coszedges = Set(Vector{Rational}())
    azimuthedges = Set(Vector{Rational}())
    energies = nothing
    for (i,j) in collect(zip(idx[1:end], [idx[2:end]..., lastindex(lines)+1]))
        bininfo = parsefluxbininfo(lines[i])
        coszedges = union(coszedges, Set(bininfo.CosZenithBin))
        azimuthedges = union(azimuthedges, Set(bininfo.AzimutBin))
        bfr = IOBuffer()
        for l in lines[i+2:j-1]
            write(bfr, strip(l))
            write(bfr, "\n")
        end
        seekstart(bfr)
        tmp = readdlm(bfr, ' ','\n')
        fluxes[bininfo] = tmp[:,2:end]
        if isnothing(energies)
            energies = tmp[:,1]
        end
    end
    coszedges = sort(collect(coszedges))
    azimuthedges = sort(collect(azimuthedges))
    data = zeros(Float64, length(coszedges)-1, length(azimuthedges)-1, length(energies), 4)
    for (k,v) in fluxes
        coszbin = searchsortedfirst(coszedges, mean(k.CosZenithBin))-1
        azimutbin = searchsortedfirst(azimuthedges, mean(k.AzimutBin))-1
        data[coszbin, azimutbin, :, :] = v
    end
    retval = Vector{FluxTable}()
    for (i,p) in enumerate([Particle(14), Particle(-14), Particle(12), Particle(-12)])
        tmp = data[:,:,:,i]
        # itp = interpolate(tmp, BSpline(Cubic(Line(OnGrid()))))
        # sitp = scale(itp, _makerange(_getbinmids(coszedges)), _makerange(_getbinmids(azimuthedges)), _makerange(log10.(energies)))
        # etp = extrapolate(sitp, Flat())
        push!(retval, FluxTable(coszedges, azimuthedges, energies, p, tmp))
    end
    retval
end

function readfluxfile(filepath::String)
    io = open(filepath, "r")
    ret = readfluxfile(io)
    close(io)
    ret
end

function _getbinindex(binedges, x)
    if x ≈ first(binedges)
        return 1
    else
        return searchsortedfirst(binedges, x)-1
    end
end

const lru_energy_etp = LRU{NuFlux.FluxTable, Interpolations.Extrapolation}(maxsize=20000)

"""
$(SIGNATURES)

# Arguments
- `flux`:       Flux data
- `energy`:     Energy in GeV
- `interpol`:   Interpolate the data
"""
function flux(f::FluxTable, energy::S; interpol::Bool=false) where {S <: Real}
    if interpol
        etp = get!(lru_energy_etp, f) do
            tmp = mean(f.flux, dims=(1,2))[1,1,:]
            itp = interpolate(tmp, BSpline(Cubic(Line(OnGrid()))))
            sitp = scale(itp, _makerange(log10.(f.energies)))
            extrapolate(sitp, Flat())
        end
        return etp(log10(energy))
    else
        idx_energy = findmin(abs.(f.energies .- energy))[2]
        return mean(f.flux[:, :, idx_energy])
    end
end

const lru_zenith_energy_etp = LRU{NuFlux.FluxTable, Interpolations.Extrapolation}(maxsize=20000)

"""
$(SIGNATURES)

# Arguments
- `flux`:       Flux data
- `energy`:     Energy in GeV
- `cosθ`:       Cosine of the zenith angle
- `interpol`:   Interpolate the data
"""
function flux(f::FluxTable, energy::S, cosθ::T; interpol::Bool=false) where {S,T <: Real}
    if interpol
        etp = get!(lru_zenith_energy_etp, f) do
            tmp = mean(f.flux, dims=2)[:,1,:]
            itp = interpolate(tmp, BSpline(Cubic(Line(OnGrid()))))
            sitp = scale(itp, _makerange(_getbinmids(f.coszentihbinedges)), _makerange(log10.(f.energies)))
            extrapolate(sitp, Flat())
        end
        return etp(cosθ, log10(energy))
    else
        idx_energy = findmin(abs.(f.energies .- energy))[2]
        idx_coszenith = _getbinindex(f.coszentihbinedges, cosθ)
        return mean(f.flux[idx_coszenith, :, idx_energy])
    end
end

const lru_zenith_azimuth_energy_etp = LRU{NuFlux.FluxTable, Interpolations.Extrapolation}(maxsize=20000)

"""
$(SIGNATURES)

# Arguments
- `flux`:       Flux data
- `energy`:     Energy in GeV
- `cosθ`:       Cosine of the zenith angle
- `ϕ`:          Azimuth angle
- `interpol`:   Interpolate the data
"""
function flux(f::FluxTable, energy::S, cosθ::T, ϕ::U; interpol::Bool=false) where {S,T,U <: Real}
    if interpol
        etp = get!(lru_zenith_azimuth_energy_etp, f) do
            itp = interpolate(f.flux, BSpline(Cubic(Line(OnGrid()))))
            sitp = scale(itp, _makerange(_getbinmids(f.coszentihbinedges)), _makerange(_getbinmids(f.azimuthbinedges)), _makerange(log10.(f.energies)))
            extrapolate(sitp, Flat())
        end
        return etp(cosθ, ϕ, log10(energy))
    else
        idx_energy = findmin(abs.(f.energies .- energy))[2]
        idx_azimuth = _getbinindex(f.azimuthbinedges, ϕ)
        idx_coszenith = _getbinindex(f.coszentihbinedges, cosθ)
        return f.flux[idx_coszenith, idx_azimuth, idx_energy]
    end
end

end # module


text:
"""
mutable struct STL{T} <: AbstractGeometry{T, 3}

Initial attempt at a complex, mesh based volume
"""

using MeshIO
using GeometryTypes
using FileIO
using NPZ
using LinearAlgebra
using Distances

struct mesh_surface{T} <: AbstractVolumePrimitive{T, 3}
    x::Array{Float32}
    y::Array{Float32}
    z::Array{Float32}
    vertices::Array{Int32}
    tri_vertices::Array{Float32}
    translate::Union{CartesianVector{T}, Missing}
end

function in(pt::CartesianPoint{T}, g::mesh_surface{T})::Bool where {T}

    use_nthreads::Int = Base.Threads.nthreads()

    triangles = permutedims(g.tri_vertices, [2, 3, 1])

    point = [pt[1] pt[2] pt[3]; pt[1] pt[2] pt[3]; pt[1] pt[2] pt[3]]

    M = triangles .- point
    M2 = (M).^2
    M3 = dropdims(sum(M2, dims=2), dims=2)
    M4 = sqrt.(M3)

    winding_sum = 0.0
    for (index_val, M4val) in enumerate(eachcol(M4)) #mapslices might work here?
        pA = det(M[:,:,index_val])
        pB = M4val[1] * M4val[2] * M4val[3]

        pC = M4val[3] * dot(M[1,:,index_val], M[2,:,index_val])
        pD = M4val[2] * dot(M[3,:,index_val], M[1,:,index_val])
        pE = M4val[1] * dot(M[2,:,index_val], M[3,:,index_val])

        pF = pC + pD + pE + pB

        winding_angle = atan(pA,pF)
        winding_sum = winding_sum + winding_angle

    end

    truth = ((1.9 * pi) < winding_sum < (2.1* pi))

#	open("truefalse.txt", "a") do io
#	   println(io, string(pt[1]), " ", string(pt[2]), " ", string(pt[3]), " ", string(truth))
           #write(io, pt[1], " ", pt[2], " ", pt[3], " ", truth)
#   end

    return truth

end

@inline in(pt::CylindricalPoint, g::mesh_surface)::Bool = in(CartesianPoint(pt), g)

function mesh_surface{T}(dict::Union{Dict{Any, Any}, Dict{String, Any}}, inputunit_dict::Dict{String,Unitful.Units})::mesh_surface{T} where {T <: SSDFloat}

    meshname = dict["filepath"]
    println("Reading mesh: ", meshname )
    mesh = load(meshname)
    println("STL mesh read")
    scale = [0.001, 0.001, 0.001] # should probably be user defined in a dict?
    stl_vertices = GeometryTypes.vertices(mesh)
    total_points = length(stl_vertices)
    stl_faces = GeometryTypes.faces(mesh)
    total_faces = length(stl_faces)

    stl_tri_points = zeros(Float32, total_faces, 3, 3)
    stl_face_points = zeros(Int32, total_faces, 3)
    x = zeros(Float32, total_points)
    y = zeros(Float32, total_points)
    z = zeros(Float32, total_points)

    index_counter = 1
    for vertex in stl_vertices
        x[index_counter] = vertex[1]*scale[1]
        y[index_counter] = vertex[2]*scale[2]
        z[index_counter] = vertex[3]*scale[3]
        index_counter = index_counter + 1
    end

    index_counter = 1
    for stl_face in stl_faces
        stl_tri_points[index_counter, :, :] =  [stl_vertices[stl_face[1]][1]*scale[1] stl_vertices[stl_face[1]][2]*scale[2] stl_vertices[stl_face[1]][3]*scale[3]; stl_vertices[stl_face[2]][1]*scale[1] stl_vertices[stl_face[2]][2]*scale[2] stl_vertices[stl_face[2]][3]*scale[3]; stl_vertices[stl_face[3]][1]*scale[1] stl_vertices[stl_face[3]][2]*scale[2] stl_vertices[stl_face[3]][3]*scale[3]]
        stl_face_points[index_counter, :] =  [stl_face[1], stl_face[2], stl_face[3]]

        index_counter = index_counter + 1

    end

    if haskey(dict, "translate")
        translate = CartesianVector{T}(
            haskey(dict["translate"],"x") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["x"]) * inputunit_dict["length"] ))) : 0.0,
            haskey(dict["translate"],"y") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["y"]) * inputunit_dict["length"] ))) : 0.0,
            haskey(dict["translate"],"z") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["z"]) * inputunit_dict["length"] ))) : 0.0)
    else
        translate = CartesianVector{T}(0.0, 0.0, 0.0)
    end

    return mesh_surface{T}(
        x,
        y,
        z,
        stl_face_points,
        stl_tri_points,
        translate)
end

function Geometry(T::DataType, t::Val{:mesh_surface}, dict::Dict{Any, Any}, inputunit_dict::Dict{String,Unitful.Units})
    return mesh_surface{T}(dict, inputunit_dict)
end

function get_important_points(g::mesh_surface{T})::NTuple{3, Vector{T}} where {T <: SSDFloat}
    v1::Vector{T} = T[g.x[1], g.x[2]] #[g.x[1], g.x[2]]
    v2::Vector{T} = T[g.y[1], g.y[2]]
    v3::Vector{T} = T[g.z[1], g.z[2]]
    return v1, v2, v3
end

function get_important_points(g::mesh_surface{T}, ::Val{:r})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T.(abs.([g.x[1], g.x[2], g.y[1], g.y[2]])))
end
function get_important_points(g::mesh_surface{T}, ::Val{:φ})::Vector{T} where {T <: SSDFloat}
    return T[]
end
function get_important_points(g::mesh_surface{T}, ::Val{:z})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[g.z[1], g.z[2]])
end
function get_important_points(g::mesh_surface{T}, ::Val{:x})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[g.x[1], g.x[2]])
end
function get_important_points(g::mesh_surface{T}, ::Val{:y})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[g.y[1], g.y[2]])
end

function vertices(cb::mesh_surface{T})::Vector{CartesianPoint{T}} where {T <: SSDFloat}

    v = CartesianPoint{T}[]

    #println("Making vertices")
    index_counter = 1
    for x_val in cb.x
        push!(v, CartesianPoint{T}(cb.x[index_counter], cb.y[index_counter], cb.z[index_counter]))

        index_counter = index_counter + 1
    end

    ismissing(cb.translate) ? nothing : v = map(x -> x + cb.translate, v)
    #println("Vertices produced")
    return v
end

function LineSegments(cb::mesh_surface{T})::Vector{LineSegment{T, 3, :cartesian}} where {T}
    v::Vector{CartesianPoint{T}} = vertices(cb)

    ls = LineSegment{T, 3, :cartesian}[]

    for face_indices in eachrow(cb.vertices)
        push!(ls, LineSegment(v[face_indices[1]], v[face_indices[2]]))
        push!(ls, LineSegment(v[face_indices[2]], v[face_indices[3]]))
        push!(ls, LineSegment(v[face_indices[3]], v[face_indices[1]]))
    end

    return ls
end

@recipe function f(cb::mesh_surface{T}) where {T <: SSDFloat}
    label-->"mesh_surface"
    ls = LineSegments(cb)
    @series begin
        ls
    end
end

function sample(cb::mesh_surface{T}, stepsize::Vector{T})  where {T <: SSDFloat}
    samples  = CartesianPoint{T}[]
    for x in cb.x[1] : stepsize[1] : cb.x[2]
        for y in cb.y[1] :stepsize[2] : cb.y[2]
            for z in cb.z[1] : stepsize[3] : cb.z[2]
                push!(samples,CartesianPoint{T}(x, y, z))
            end
        end
    end
    ismissing(cb.translate) ? nothing : samples = map(x -> x + cb.translate, samples)
    return samples
end

function (+)(b::mesh_surface{T}, translate::Union{CartesianVector{T},Missing})::mesh_surface{T} where {T <: SSDFloat}
    if ismissing(translate)
        return b
    elseif ismissing(b.translate)
        return mesh_surface(b.x, b.y, b.z, b.vertices, translate)
    else
        return mesh_surface(b.x, b.y, b.z, b.vertices, b.tri_vertices, b.translate + translate)
    end
 end


text:
using Automata
using Images

function exec()
    println("Starting...")

    # @time automaton(ElementaryRule(30), 10000, reshape(Bool[ 1 ], 1, 1))
    # save("30.png", Gray.(1 .- Float32.(automaton(ElementaryReversibleRule(90), 100))))
    # save("30R.png", Gray.(1 .- Float32.(automaton(ElementaryReversibleRule(30), 100))))
    # save("2530_2.png", Gray.(1 .- Float32.(automaton(TotalisticRule(2530, 4), 100)) ./ 3))
    # save("d3C.png", Gray.(1 .- Float32.(automaton(ContinuousRule(@inline (x -> (x * (3/2)) % 1)), 100))))
    # automaton(ElementaryRule(30), 10)
    # @time automaton(ElementaryRule(30), 10000)
    # @time automaton(ElementaryRule(30), 1000)
    # @time elementary(30, 1000)
end

exec()


text:
using SpikingNN
using Plots

# simulation parameters
T = 1000

# neuron parameters
vᵣ = 0
τᵣ = 1.0
vth = 1.0

# create population
# neuron 1, 2 excite neuron 3
# neuron 3 inhibits neuron 1, 2
weights = [ 0  0  1;
            0  0  1;
           -5 -5  0]
pop = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
                          synapse = Synapse.Alpha,
                          threshold = () -> Threshold.Ideal(vth))

# create input currents
low = ConstantRate(0.1)
high = ConstantRate(0.99)
switch(t; dt = 1) = (t < Int(T/2)) ? low(t; dt = dt) : high(t; dt = dt)
n1synapse = QueuedSynapse(Synapse.Alpha())
n2synapse = QueuedSynapse(Synapse.Alpha())
excite!(n1synapse, filter(x -> x != 0, [low(t) for t = 1:T]))
excite!(n2synapse, filter(x -> x != 0, [switch(t) for t = 1:T]))

# simulate
voltages = Dict([(i, Float64[]) for i in 1:3])
cb = () -> begin
    for id in 1:size(pop)
        push!(voltages[id], getvoltage(pop[id]))
    end
end
@time outputs = simulate!(pop, T; cb = cb, inputs = [n1synapse, n2synapse, (t; dt) -> 0])

rasterplot(outputs, label = ["Input 1", "Input 2", "Inhibitor"])
title!("Raster Plot")
xlabel!("Time (sec)")


text:
using Plots

export true_positives
export true_negatives
export false_positives
export false_negatives
export precision_classification
export recall
export accuracy
export f1_score
export confusion_matrix
export predict_from_threshold
export best_f1_search

"true positives. If rate: sensitivity (true positives / actual positives)"
true_positives(actual::Array{Int64,1}, pred::Array{Int64,1}; rate = false) =
    sum(actual .== pred .== 1) / (rate ? sum(actual .== 1) : 1)

"true negatives. If rate: specificity (true negatives / actual negatives)"
true_negatives(actual::Array{Int64,1}, pred::Array{Int64,1}; rate = false) =
    sum(actual .== pred .== 0) / (rate ? sum(actual .== 1) : 1)

"false positives. If rate: fall-out (false positives / actual negatives)"
false_positives(actual::Array{Int64,1}, pred::Array{Int64,1}; rate = false) =
    sum((actual .== 0) .& (pred .== 1)) / (rate ? sum(actual .== 0) : 1)

"true negatives. If rate: (false negatives / actual positives)"
false_negatives(actual::Array{Int64,1}, pred::Array{Int64,1}; rate = false) =
    sum((actual .== 1) .& (pred .== 0)) / (rate ? sum(actual .== 1) : 1)

"precision (false positives / (true positives + false positives))"
function precision_classification(actual::Array{Int64,1}, pred::Array{Int64,1})
    tp = true_positives(actual, pred)
    fp = false_positives(actual, pred)
    return tp / (tp + fp)
end

"recall (true positives / (true positives + false negatives))"
function recall(actual::Array{Int64,1}, pred::Array{Int64,1})
    tp = true_negatives(actual, pred)
    fn = false_negatives(actual, pred)
    return tp / (tp + fn)
end

"proportion of classifications made correctly"
accuracy(actual::Array{Int64,1}, pred::Array{Int64,1}) = sum(actual .== pred) / length(actual)

"F1 score, harmonic mean of precision and recall, 2 * ((precision * recall) / (precision + recall))"
function f1_score(actual::Array{Int64,1}, pred::Array{Int64,1})
    prec = precision_classification(actual, pred)
    rec = recall(actual, pred)
    return 2 * ((prec * rec) / (prec + rec))
end

"""
_given array of 1s and 0s, actuals and predictions, returns confusion matrix_

#### parameters:
    actual : Array{Int}
        actual labels, ∈ [0,1]
    pred : Array{Float64}
        an array of prediction values 0 <= x <= 1
    ratio : Boolean
        if true, return matrix as ratios, sum = 1, else return absolute numbers

#### returns: Array{Float64,2}
    confusion matrix of form:
    [tp fp]
    [fn tn]
    X axis = actual
    y axis = predicted
"""
function confusion_matrix(
    actual::Array{Int64,1},
    pred::Array{Int64,1};
    ratio = false,
)
    [
     true_positives(actual, pred) false_positives(actual, pred)
     false_negatives(actual, pred) true_negatives(actual, pred)
    ] ./ (ratio ? (length(actual) |> x -> round.(x, digits = 3)) : 1)
end

"given an array of predicted probabilities, return 1 or 0 depending on threshold"
predict_from_threshold(pred::Array{Float64}, threshold::Float64) = [perc > threshold ? 1 : 0 for perc in pred]

"""
_given array of true labels and % predictions, return best threshold and F1 score, plot of F1 score a>_

#### parameters:
    actual : Array{Int}
        actual labels, ∈ [0,1]
    pred : Array{Float64}
        an array of prediction values 0 <= x <= 1

#### returns: Tuple{Float64, Float64, Plots.Plot{Plots.GRBackend}}
    tuple of best threshold, best F1 score, plot of F1 scores
"""
function best_f1_search(actual::Array{Int}, pred::Array{Float64})
    x = 0:0.01:1
    accuracies = [sum(predict_from_threshold(pred, threshold) .== actual) / length(actual) for threshold in x]
    f1s = [f1_score(actual, predict_from_threshold(pred, threshold)) for threshold in x]

    p1 = plot(
        x,
        accuracies,
        label = "Accuracy",
        xlabel = "threshold",
        ylabel = "F1 Score",
        title = "F1 Score",
        legend = :outerbottomright,
        ylim = (0, 1),
    )
    plot!(x, f1s, label = "F1")

    best_f1, best_threshold_index = findmax([isnan(f1) ? 0 : f1 for f1 in f1s])
    return x[best_threshold_index], best_f1, p1
end


text:
# Converting composite objects to Expr nodes for pretty printing.

struct QuoteContext
    limit::Bool
    vars::IdDict{Any, Symbol}
    colors::Vector{Symbol}

    QuoteContext(;
                 limit = false,
                 vars = IdDict{Any, Symbol}(),
                 colors = [:normal]) =
        new(limit, vars, colors)
end


text:
# Evaluate the equality constraint function
function eval_ceq!(gE::Base.ReshapedArray,gE_N::SubArray,solver::Solver,res::DircolResults)
    N = solver.N
    pE = size(gE,1)
    pE_N = size(gE_N,1)

    # Custom constraints
    X,U = res.X,res.U
    if pE > 0
        for k = 1:N-1
            gE[:,k] = obj.cE(X[:,k],U[:,k])
        end
    end
    if pE_N > 0
        gE_N .= obj.cE(X[:,N])
    end
end

function eval_c!(gI::Base.ReshapedArray,gI_N::SubArray,solver::Solver,res::DircolResults)
    N = solver.N
    pI = size(gI,1)
    pI_N = size(gI_N,1)

    # Custom constraints
    X,U = res.X,res.U
    if pI > 0
        for k = 1:N-1
            gE[:,k] = obj.cI(X[:,k],U[:,k])
        end
    end
    if pI_N > 0
        gE_N .= obj.cI(X[:,N])
    end
end

"""
$(SIGNATURES)
Generate the custom function to be passed into SNOPT, as well as `eval_f` and
`eval_g` used to calculate the objective and constraint functions.

# Arguments
* model: TrajectoryOptimization.Model for the dynamics
* obj: TrajectoryOptimization.ConstrainedObjective to describe the cost function
* dt: time step
* pack: tuple of important sizes (n,m,N) -> (num states, num controls, num knot points)
* method: Collocation method. Either :trapezoid or :hermite_simpson_separated
* grads: Specifies the gradient information provided to SNOPT.
    :none - returns no gradient information to SNOPT
    :auto - uses ForwardDiff to calculate gradients
    :quadratic - uses functions exploiting quadratic cost functions
"""
function gen_usrfun(solver::Solver, results::DircolResults, method::Symbol; grads=:none)::Function
    n,N = size(results.X)
    m,N_ = size(results.U_)
    NN = (n+m)N  # Number of decision variables
    dt = solver.dt

    # Count constraints
    pI = 0
    pE = (N-1)*n # Collocation constraints
    p_colloc = pE

    pI_obj, pE_obj = count_constraints(solver.obj)
    pI_c,   pE_c   = pI_obj[2], pE_obj[2]  # Number of custom stage constraints
    pI_N_c, pE_N_c = pI_obj[4], pE_obj[4]  # Number of custom terminal constraints

    pI += pI_c*(N-1) + pI_N_c  # Add custom inequality constraints
    pE += pE_c*(N-1) + pE_N_c  # Add custom equality constraints

    # Allocate Arrays
    c = zeros(pI)
    ceq = zeros(pE)
    grad_f = zeros((n+m)N)
    jacob_ceq = spzeros((N-1)*n,N*(n+m))

    g_colloc = view(ceq,1:p_colloc)
    gE = reshape(view(ceq,p_colloc .+ (1:pI_c*(N-1))),pI_c,N-1)
    gE_N = view(ceq,p_colloc+pI_c*(N-1) .+ (1:pE_N_c))

    gI = reshape(view(c,1:pI_c*(N-1)),pI_c,N-1)
    gI_N = view(c,pE_c .+ (1:pI_N_c))

    # User defined function passed to SNOPT (via Snopt.jl)
    function usrfun(Z::Vector{Float64})
        results.Z .= Z

        # Update Derivatives
        update_derivatives!(solver,results,method)

        # Get points used in integration
        get_traj_points!(solver,results,method)
        get_traj_points_derivatives!(solver,results,method)

        # Objective Function (Cost)
        J = cost(solver,results)

        # Collocation Constraints
        collocation_constraints!(g_colloc, solver, results, method)
        eval_ceq!(gE,gE_N,solver,results) # No inequality constraints for now
        eval_c!(gI,gI_N,solver,results)

        fail = false

        if grads == :none
            return J, c, ceq, fail
        else
            if pE_c > 0 || pE_N_c > 0 || pI_c > 0 || pI_N_c > 0
                error("Gradients not defined for custom constraints yets")
            end

            # Calculate dynamics jacobians at integration points
            update_jacobians!(solver,results,method)

            # Gradient of Objective
            cost_gradient!(grad_f, solver,results,method)

            # Constraint Jacobian
            # constraint_jacobian!(jacob_ceq, solver,results,method)
            jacob_ceq = constraint_jacobian(solver, results, method::Symbol)
            jacob_c = Float64[]

            return J, c, ceq, grad_f, jacob_c, jacob_ceq, fail
        end

    end

    function usrfun(Z::Vector{Float64},mode::Symbol)
        if mode == :ceq
            row,col = constraint_jacobian_sparsity(solver,method)
        else
            row,col = [],[]
        end
        return row,col
    end

    # Specify sparsity structure
    # gceq = constraint_jacobian_sparsity(solver,results,method)
    # gceq_c = spzeros(pE_c,NN)
    # gceq_c[:,1:NN-(n+m)] = 1
    # gceq_N = spzeros(pE_N_c,NN)
    # gceq_N[:,NN-(n+m)+1:end] = 1
    # gceq = [gceq; gceq_c; gceq_N]
    #
    # gc_c = spzeros(pI_c,NN)
    # gc_c[:,1:NN-(n+m)] = 1
    # gc_N = spzeros(pI_N_c,NN)
    # gc_N[:,NN-(n+m)+1:end] = 1
    # gc = [gc_c; gc_N]
    #
    # function usrfun(s::Symbol)
    #     if s == :gceq
    #         return gceq
    #     elseif s == :gc
    #         return gc
    #     end
    # end

    return usrfun
end

function constraint_jacobian_sparsity(solver::Solver, method::Symbol)
    sparsity = collocation_constraint_jacobian_sparsity(solver, method)
    r,c,inds = findnz(sparsity)
    v = sortperm(inds)
    rows = r[v]
    cols = c[v]
    return rows, cols
end

"""
$(SIGNATURES)
Solve a trajectory optimization problem with direct collocation

# Arguments
* model: TrajectoryOptimization.Model for the dynamics
* obj: TrajectoryOptimization.ConstrainedObjective to describe the cost function
* dt: time step. Used to determine the number of knot points. May be modified
    by the solver in order to achieve an integer number of knot points.
* method: Collocation method.
    :midpoint - Zero order interpolation on states and controls.
    :trapezoidal - First order interpolation on states and zero order on control.
    :hermite_simpson_separated - Hermite Simpson collocation with the midpoints
        included as additional decision variables with constraints
    :hermite_simpson - condensed verision of Hermite Simpson. Currently only
        supports grads=:auto or grads=:none (recommended)
* grads: Specifies the gradient information provided to SNOPT.
    :none - returns no gradient information to SNOPT
    :auto - uses ForwardDiff to calculate gradients
    :quadratic - uses functions exploiting quadratic cost functions
"""
function solve_snopt(solver::Solver,X0::Matrix,U0::Matrix;
        method::Symbol=:auto, grads::Symbol=:quadratic, start=:cold)

    N,N_ = get_N(solver,method)
    n,m = get_sizes(solver)

    # Create results structure
    results = DircolResults(n,m,solver.N,method)
    var0 = DircolVars(X0,U0)
    Z0 = var0.Z

    # Generate the objective/constraint function and its gradients
    usrfun = gen_usrfun(solver, results, method, grads=grads)

    # Set up the problem
    lb,ub = get_bounds(solver,method)

    # Set options
    options = Dict{String, Any}()
    options["Derivative option"] = 0
    options["Verify level"] = 1
    options["Major feasibility tol"] = solver.opts.constraint_tolerance
    # options["Minor optimality  tol"] = solver.opts.eps_intermediate
    options["Major optimality  tol"] = solver.opts.cost_tolerance
    sumfile = "logs/snopt-summary.out"
    printfile = "logs/snopt-print.out"

    # Solve the problem
    if solver.opts.verbose
        println("DIRCOL with $method")
        println("Passing Problem to SNOPT...")
    end
    row,col = constraint_jacobian_sparsity(solver,method)
    # prob = Snopt.createProblem(usrfun, Z0, lb, ub, iE=row, jE=col)
    # prob = Snopt.createProblem(usrfun, Z0, lb, ub)
    # prob.x = Z0
    t_eval = @elapsed z_opt, fopt, info = snopt(usrfun, Z0, lb, ub, options,
        printfile=printfile, sumfile=sumfile)
    stats = parse_snopt_summary()
    stats["info"] = info
    stats["runtime"] = t_eval
    # @time snopt(usrfun, Z0, lb, ub, options, start=start)
    # xopt, fopt, info = Z0, Inf, "Nothing"
    sol = DircolVars(z_opt,n,m,N)

    if solver.opts.verbose
        println(info)
    end
    return sol, stats
end

"""
$(SIGNATURES)
MESH REFINEMENT:
Solve by warm starting with a coarse time step and warm-starting the solver
with the previous solution

# Arguments
* mesh: vector of step sizes for refinement prior to step size specified by solver.
"""
function solve_dircol(solver::Solver, X0::Matrix, U0::Matrix, mesh::Vector;
        method::Symbol=:auto, grads::Symbol=:quadratic, start=:cold)
    x_opt, u_opt = X0, U0
    stats = Dict{String,Any}("iterations"=>0, "major iterations"=>0, "objective calls"=>0, "info"=>Vector{String}())
    tic()
    for dt in mesh
        solver.opts.verbose ? println("Refining mesh at dt=$dt") : nothing
        solver_mod = Solver(solver,dt=dt)
        x_int, u_int = interp_traj(solver_mod.N, solver.obj.tf, x_opt, u_opt)
        x_opt, u_opt, f_opt, stats_run = solve_dircol(solver_mod, x_int, u_int, method=method, grads=grads, start=start)
        for key in keys(stats)
            if key == "info"
                push!(stats["info"], stats_run["info"])
            else
                stats[key] += stats_run[key]
            end
        end
        start = :warm  # Use warm starts after the first one
    end

    # Run the original time step
    x_int, u_int = interp_traj(solver.N, solver.obj.tf, x_opt, u_opt)
    x_opt, u_opt, f_opt, stats_run = solve_dircol(solver, x_int, u_int, method=method, grads=grads, start=:warm)
    for key in keys(stats)
        if key == "info"
            push!(stats["info"], stats_run["info"])
        else
            stats[key] += stats_run[key]
        end
    end
    stats["runtime"] = toq()
    return x_opt, u_opt, f_opt, stats
end

function solve_dircol(solver::Solver, mesh::Vector;
        method::Symbol=:auto, grads::Symbol=:quadratic, start=:cold)
    # Constants
    N = solver.N
    N = convert_N(N,method)

    X0, U0 = get_initial_state(solver.obj,N)
    solve_dircol(solver, X0, U0, mesh, method=method, grads=grads, start=start)
end

"""
$(SIGNATURES)
Extract important information from the SNOPT output file(s)
"""
function parse_snopt_summary(file="logs/snopt-summary.out")
    props = Dict()

    function stash_prop(ln::String,prop::String,prop_name::String=prop,vartype::Type=Float64)
        if occursin(prop,ln)
            loc = findfirst(prop,ln)
            if vartype <: Real
                val = split(ln[loc[end]+1:end])[1]
                val = convert(vartype,parse(Float64,val))
                props[prop_name] = val
                return true
            end
        end
        return false
    end

    open(file) do f
        for ln in eachline(f)
            stash_prop(ln,"No. of iterations","iterations",Int64)
            stash_prop(ln,"No. of major iterations","major iterations",Int64)
            stash_prop(ln,"No. of calls to funobj","objective calls",Int64)
            stash_prop(ln,"Nonlinear constraint violn","c_max")
        end
    end
    return props
end


text:
using Test
using MPI
using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.Mesh.DSS
using KernelAbstractions

DA = ClimateMachine.array_type()

function test_dss()
    FT = Float64
    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    @assert csize == 1

    N = (4, 4, 5)
    brickrange = (0:2, 5:6, 0:1)
    periodicity = (false, false, false)
    nvars = 1

    topl = StackedBrickTopology(
        comm,
        brickrange,
        periodicity = periodicity,
        boundary = ((1, 2), (3, 4), (5, 6)),
        connectivity = :full,
    )
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = DA,
        polynomialorder = N,
    )

    Nq = N .+ 1
    Np = prod(Nq)
    Q = MPIStateArray{FT}(
        comm,
        DA,
        Np,
        nvars,
        length(topl.elems),
        realelems = topl.realelems,
        ghostelems = topl.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topl.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
    )
    realelems = Q.realelems
    ghostelems = Q.ghostelems

    ldof = DA{FT, 1}(reshape(1:Np, Np))

    for ivar in 1:nvars
        for ielem in realelems
            Q.data[:, ivar, ielem] .= ldof .+ FT((ielem - 1) * Np)
        end
        Q.data[:, ivar, ghostelems] .= FT(0)
    end

    pre_dss = Array(Q.data)

    dss!(Q, grid)
    #---------Tests-------------------------
    nodes = reshape(1:Np, Nq)
    interior = nodes[2:(Nq[1] - 1), 2:(Nq[2] - 1), 2:(Nq[3] - 1)][:]
    vertmap = Array(grid.vertmap)
    edgemap = Array(grid.edgemap)
    facemap = Array(grid.facemap)
    post_dss = Array(Q.data)

    ne1r = 1:max(Nq[1] - 2, 0)
    ne2r = 1:max(Nq[2] - 2, 0)
    ne3r = 1:max(Nq[3] - 2, 0)

    nf1r = 1:max((Nq[2] - 2) * (Nq[3] - 2), 0)
    nf2r = 1:max((Nq[1] - 2) * (Nq[3] - 2), 0)
    nf3r = 1:max((Nq[1] - 2) * (Nq[2] - 2), 0)

    compare(el1, i1, el2, i2, efmap, rng, post, pre) =
        post[efmap[rng, i1, 1], 1, el1] == pre[efmap[rng, i2, 1], 1, el2]

    compare(el1, i1, el2, i2, el3, i3, efmap, rng, post, pre) =
        post[efmap[rng, i1, 1], 1, el1] ==
        pre[efmap[rng, i2, 1], 1, el2] .+ pre[efmap[rng, i3, 1], 1, el3]

    el1, el2 = 1, 2
    # Element # 1 --------------------------------------------------------------------------
    # interior dof should not be affected
    @test pre_dss[interior, 1, 1] == post_dss[interior, 1, 1]
    # vertex check
    @test post_dss[vertmap, 1, 1] == [
        pre_dss[vertmap[1], 1, 1],
        pre_dss[vertmap[2], 1, 1] + pre_dss[vertmap[1], 1, 2],
        pre_dss[vertmap[3], 1, 1],
        pre_dss[vertmap[4], 1, 1] + pre_dss[vertmap[3], 1, 2],
        pre_dss[vertmap[5], 1, 1],
        pre_dss[vertmap[6], 1, 1] + pre_dss[vertmap[5], 1, 2],
        pre_dss[vertmap[7], 1, 1],
        pre_dss[vertmap[8], 1, 1] + pre_dss[vertmap[7], 1, 2],
    ]
    # edge check
    @test compare(el1, 1, el1, 1, edgemap, ne1r, post_dss, pre_dss)          # edge  1 (unaffected)
    @test compare(el1, 2, el1, 2, edgemap, ne1r, post_dss, pre_dss)          # edge  2 (unaffected)
    @test compare(el1, 3, el1, 3, edgemap, ne1r, post_dss, pre_dss)          # edge  3 (unaffected)
    @test compare(el1, 4, el1, 4, edgemap, ne1r, post_dss, pre_dss)          # edge  4 (unaffected)

    @test compare(el1, 5, el1, 5, edgemap, ne2r, post_dss, pre_dss)          # edge  5 (unaffected)
    @test compare(el1, 6, el1, 6, el2, 5, edgemap, ne2r, post_dss, pre_dss) # edge  6 (shared edge)
    @test compare(el1, 7, el1, 7, edgemap, ne2r, post_dss, pre_dss)          # edge  7 (unaffected)
    @test compare(el1, 8, el1, 8, el2, 7, edgemap, ne2r, post_dss, pre_dss) # edge  8 (shared edge)

    @test compare(el1, 9, el1, 9, edgemap, ne3r, post_dss, pre_dss)          # edge  9 (unaffected)
    @test compare(el1, 10, el1, 10, el2, 9, edgemap, ne3r, post_dss, pre_dss) # edge 10 (shared edge)
    @test compare(el1, 11, el1, 11, edgemap, ne3r, post_dss, pre_dss)          # edge 11 (unaffected)
    @test compare(el1, 12, el1, 12, el2, 11, edgemap, ne3r, post_dss, pre_dss) # edge 12 (shared edge)

    # face check
    @test compare(el1, 1, el1, 1, facemap, nf1r, post_dss, pre_dss)         # face 1 (unaffected)
    @test compare(el1, 2, el1, 2, el2, 1, facemap, nf1r, post_dss, pre_dss) # face 2 (shared face)

    @test compare(el1, 3, el1, 3, facemap, nf2r, post_dss, pre_dss)         # face 3 (unaffected)
    @test compare(el1, 4, el1, 4, facemap, nf2r, post_dss, pre_dss)         # face 4 (unaffected)

    @test compare(el1, 5, el1, 5, facemap, nf3r, post_dss, pre_dss)         # face 5 (unaffected)
    @test compare(el1, 6, el1, 6, facemap, nf3r, post_dss, pre_dss)         # face 6 (unaffected)

    # Element # 2 --------------------------------------------------------------------------
    # interior dof should not be affected
    @test pre_dss[interior, 1, 2] == post_dss[interior, 1, 2]
    # vertex check
    @test post_dss[vertmap, 1, 2] == [
        pre_dss[vertmap[1], 1, 2] + pre_dss[vertmap[2], 1, 1],
        pre_dss[vertmap[2], 1, 2],
        pre_dss[vertmap[3], 1, 2] + pre_dss[vertmap[4], 1, 1],
        pre_dss[vertmap[4], 1, 2],
        pre_dss[vertmap[5], 1, 2] + pre_dss[vertmap[6], 1, 1],
        pre_dss[vertmap[6], 1, 2],
        pre_dss[vertmap[7], 1, 2] + pre_dss[vertmap[8], 1, 1],
        pre_dss[vertmap[8], 1, 2],
    ]
    # edge check
    @test compare(el2, 1, el2, 1, edgemap, ne1r, post_dss, pre_dss)          # edge  1 (unaffected)
    @test compare(el2, 2, el2, 2, edgemap, ne1r, post_dss, pre_dss)          # edge  2 (unaffected)
    @test compare(el2, 3, el2, 3, edgemap, ne1r, post_dss, pre_dss)          # edge  3 (unaffected)
    @test compare(el2, 4, el2, 4, edgemap, ne1r, post_dss, pre_dss)          # edge  4 (unaffected)

    @test compare(el2, 5, el2, 5, el1, 6, edgemap, ne2r, post_dss, pre_dss) # edge  5 (shared edge)
    @test compare(el2, 6, el2, 6, edgemap, ne2r, post_dss, pre_dss)          # edge  6 (unaffected)
    @test compare(el2, 7, el2, 7, el1, 8, edgemap, ne2r, post_dss, pre_dss) # edge  7 (shared edge)
    @test compare(el2, 8, el2, 8, edgemap, ne2r, post_dss, pre_dss)          # edge  8 (unaffected)

    @test compare(el2, 9, el2, 9, el1, 10, edgemap, ne3r, post_dss, pre_dss) # edge  9 (shared edge)
    @test compare(el2, 10, el2, 10, edgemap, ne3r, post_dss, pre_dss)          # edge 10 (unaffected)
    @test compare(el2, 11, el2, 11, el1, 12, edgemap, ne3r, post_dss, pre_dss) # edge 11 (shared edge)
    @test compare(el2, 12, el2, 12, edgemap, ne3r, post_dss, pre_dss)          # edge 12 (unaffected)

    # face check
    @test compare(el2, 1, el2, 1, el1, 2, facemap, nf1r, post_dss, pre_dss) # face 1 (shared face)
    @test compare(el2, 2, el2, 2, facemap, nf1r, post_dss, pre_dss)         # face 2 (unaffected)

    @test compare(el2, 3, el2, 3, facemap, nf2r, post_dss, pre_dss)         # face 3 (unaffected)
    @test compare(el2, 4, el2, 4, facemap, nf2r, post_dss, pre_dss)         # face 4 (unaffected)

    @test compare(el2, 5, el2, 5, facemap, nf3r, post_dss, pre_dss)         # face 5 (unaffected)
    @test compare(el2, 6, el2, 6, facemap, nf3r, post_dss, pre_dss)         # face 6 (unaffected)

end

test_dss()