text:
using Distributed

nprocs() == 1 && addprocs(1)

@everywhere module Sharding

using Discord

function main()
    c = Client(ENV["DISCORD_TOKEN"])

    add_handler!(c, AbstractEvent) do c, e
        println("Shard $(c.shard) received $(typeof(e))")
    end

    open(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    @everywhere begin
        c = Sharding.main()
        wait(c)
    end
end


text:
struct WithLabelArg{T, U} <: Distribution{T}
    base :: Distribution{U}
end

function logpdf(d::WithLabelArg{T, U}, x::T, collection, base_args...) where {T, U}
    if is_discrete(d.base)
        # Accumulate
        logprobs::Array{Float64, 1} = []
        for p in pairs(collection)
            (index, item) = (p.first, p.second)
            if item == x
                push!(logprobs, logpdf(d.base, index, base_args...))
            end
        end
        isempty(logprobs) ? -Inf : logsumexp(logprobs)
    else
        error("Cannot relabel a continuous distribution")
    end
end

function logpdf_grad(d::WithLabelArg{T, U}, x::T, collection, base_args...) where {T, U}
    base_arg_grads = fill(nothing, length(base_args))

    for p in pairs(collection)
        (index, item) = (p.first, p.second)
        if item == x
            new_grads = logpdf_grad(d.base, index, base_args...)
            for (arg_idx, grad) in enumerate(new_grads)
                if base_arg_grads[arg_idx] === nothing
                    base_arg_grads[arg_idx] = grad
                elseif grad !== nothing
                    base_arg_grads[arg_idx] += grad
                end
            end
        end
    end
    (nothing, nothing, base_arg_grads...)
end

function random(d::WithLabelArg{T, U}, collection, base_args...)::T where {T, U}
    collection[random(d.base, base_args...)]
end

is_discrete(d::WithLabelArg{T, U}) where {T, U} = true

(d::WithLabelArg{T, U})(collection, base_args...) where {T, U} = random(d, collection, base_args...)

function has_output_grad(d::WithLabelArg{T, U}) where {T, U}
    false
end

has_argument_grads(d::WithLabelArg{T, U}) where {T, U} = (false, has_argument_grads(d.base)...)

struct RelabeledDistribution{T, U} <: Distribution{T}
    base :: Distribution{U}
    collection::Union{AbstractArray{T}, AbstractDict{U, T}}
end

function logpdf(d::RelabeledDistribution{T, U}, x::T, base_args...) where {T, U}
    if is_discrete(d.base)
        # Accumulate
        logprobs::Array{Float64, 1} = []
        for p in pairs(d.collection)
            (index, item) = (p.first, p.second)
            if item == x
                push!(logprobs, logpdf(d.base, index, base_args...))
            end
        end
        isempty(logprobs) ? -Inf : logsumexp(logprobs)
    else
        error("Cannot relabel a continuous distribution")
    end
end

function logpdf_grad(d::RelabeledDistribution{T, U}, x::T, base_args...) where {T, U}
    base_arg_grads = fill(nothing, length(base_args))

    for p in pairs(d.collection)
        (index, item) = (p.first, p.second)
        if item == x
            new_grads = logpdf_grad(d.base, index, base_args...)
            for (arg_idx, grad) in enumerate(new_grads)
                if base_arg_grads[arg_idx] === nothing
                    base_arg_grads[arg_idx] = grad
                elseif grad !== nothing
                    base_arg_grads[arg_idx] += grad
                end
            end
        end
    end
    (nothing, base_arg_grads...)
end

function random(d::RelabeledDistribution{T, U}, base_args...)::T where {T, U}
    d.collection[random(d.base, base_args...)]
end

is_discrete(d::RelabeledDistribution{T, U}) where {T, U} = true

(d::RelabeledDistribution{T, U})(base_args...) where {T, U} = random(d, base_args...)

function has_output_grad(d::RelabeledDistribution{T, U}) where {T, U}
    false
end

has_argument_grads(d::RelabeledDistribution{T, U}) where {T, U} = has_argument_grads(d.base)

ext: jl


text:
module BrokemiaHelperCelesteNetFlagSynchronizer

using ..Ahorn, Maple

@mapdef Entity "BrokemiaHelper/CelesteNetFlagSynchronizer" CelesteNetFlagSynchronizer(x::Integer, y::Integer, flag::String="")

const placements = Ahorn.PlacementDict(
    "CelesteNet Flag Synchronizer (BrokemiaHelper)" => Ahorn.EntityPlacement(
        CelesteNetFlagSynchronizer,
        "point",
        Dict{String, Any}()
    )
)
sprite = "Ahorn/BrokemiaHelper/CelesteNetFlagSynchronizer.png"

function Ahorn.selection(entity::CelesteNetFlagSynchronizer)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CelesteNetFlagSynchronizer, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
#!/usr/bin/env julia

include("read.jl")

function sort(file=ARGS[1])
    lines=read(file)
    readCounts = Dict()
    for ln in lines
        if ismatch(r"^CHROM", ln)
            print()
        else
            cols = split(ln, '\t')
            println(cols[1:5])
            positiveSamples = 0
            for col in cols[5:end]
                read = split(col, ':')[2]
                if read >= 1
                    positiveSamples += 1
                end
            end
        end
    pair = cols[1] * positiveSamples
    push!(readCounts, pair)
    end
end

test = sort()


text:
"""
    ImpulseResponse{TF<:AbstractFloat, VCE<:CovarianceEstimator} <: StatisticalModel

Estimates for an impulse response function.
"""
struct ImpulseResponse{TF<:AbstractFloat, VCE<:CovarianceEstimator} <: StatisticalModel
    B::Vector{TF}
    SE::Vector{TF}
    T::Vector{Int}
    vce::VCE
    yname::VarName
    xname::VarName
    minhorz::Int
end

coef(f::ImpulseResponse) = f.B
stderror(f::ImpulseResponse) = f.SE

function confint(f::ImpulseResponse; level::Real=0.9, horz=Colon())
    scale = criticalvalue.(Ref(f.vce), level, view(f.T, horz))
    se = stderror(f)[horz]
    b = coef(f)[horz]
    return b .- scale .* se, b .+ scale .* se
end

"""
    irf(r::LocalProjectionResult, yname::VarName, xwname::VarName; lag::Int=0)

Return the [`ImpulseResponse`](@ref) of outcome variable `yname`
with respect to variable `xwname` based on the estimation result `r`.
If `lag` is not specified, `xwname` is assumed to be contemporaneous.
"""
function irf(r::LocalProjectionResult, yname::VarName, xwname::VarName; lag::Int=0)
    TF = typeof(r).parameters[4]
    H = length(r.T)
    B = Vector{TF}(undef, H)
    SE = Vector{TF}(undef, H)
    for h in 1:H
        B[h] = coef(r, h, xwname, yname=yname, lag=lag)
        SE[h] = sqrt(vcov(r, h, xwname, yname1=yname, lag1=lag))
    end
    return ImpulseResponse(B, SE, r.T, r.vce, yname, xwname, r.minhorz)
end

function coeftable(f::ImpulseResponse;
        level::Real=0.90, horz::Union{Integer,AbstractVector,Colon}=Colon())
    H = length(f.B)
    if horz !== Colon()
        cf = coef(f)[horz]
        se = stderror(f)[horz]
        cil, ciu = confint(f, level=level, horz=horz)
        cnames = horz .+ f.minhorz .-1
    else
        cf = coef(f)
        se = stderror(f)
        cil, ciu = confint(f, level=level)
        cnames = f.minhorz:f.minhorz+H-1
    end
    ts = cf ./ se
    pv = pvalue.(Ref(f.vce), ts, view(f.T, horz))
    levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)
    return CoefTable(Vector[cf, se, ts, pv, cil, ciu],
        ["Estimate", "Std. Error", "z", "Pr(>|z|)", "Lower $levstr%", "Upper $levstr%"],
        ["$(cnames[i])" for i = 1:length(cf)], 4, 3)
end

show(io::IO, f::ImpulseResponse) = print(io, typeof(f).name.name)

function show(io::IO, ::MIME"text/plain", f::ImpulseResponse)
    H = length(f.B)
    nr, nc = displaysize(io)
    print(io, typeof(f).name.name, " with $H horizon")
    println(io, H > 1 ? "s:" : ":")
    if H > nr-8
        show(io, coeftable(f, horz=1:nr-9))
        printstyled(io, "\n", H-nr+9, " rows of estimates omitted"; color=:cyan)
    else
        show(io, coeftable(f))
    end
end


text:
using Origin
using Test

@testset "Test for plain" begin
@origin (a=>0, b=>2) function test()
    a = collect(0:10)
    b = collect(1:5)
    @test a[0] == 0
    @test a[10] == 10
    @test b[2] == 1
    @test b[4] == 3
    @test b[end] == 5
    @test b[2:end] == [1, 2, 3, 4, 5]
    @test b[a[5]] == 4
end
test()
end

@testset "Test for nesting" begin
@origin (a=>0, b=>2) function test()
    a = collect(0:10)
    b = collect(1:5)
    @test a[0] == 0
    @test b[2] == 1
    @origin a=>100 begin
        @test a[101] == 1
    end
end
test()
end

@testset "Test for eachindex" begin
@origin (a=>0, b=>2) function test()
    a = collect(0:10)
    for i = eachindex(a)
        a[i] = i
    end
end
@test_throws BoundsError test()
end

@testset "Test for eachindex2" begin
@origin (b=>2) function test()
    a = collect(0:10)
    for i = eachindex(a)
        a[i] = i
    end
    a
end
a = test()
for i = eachindex(a)
    @test a[i] == i
end
end


text:
using Documenter, NetalignUtils

makedocs(
           format = :html,
           sitename = "NetalignUtils",
           pages = [
                    "index.md"
           ],
    modules = [NetalignUtils]
       )

deploydocs(
           repo = "github.com/vvjn/NetalignUtils.jl.git",
           target = "build",
           deps   = nothing,
           make   = nothing,
           julia = "0.6"
)


text:
"""
# module Fortran

# Examples

```jldoctest
julia>
```
"""
module Fortran

end


text:
using CharFuncPricing, Test
using TaylorSeries

@testset "SZLogContinuity" begin
 v0 = 0.08
 ρ=-0.8
 σ=0.25
 θ=0.1
 κ=3.0
 #params from Cui et al for Heston.
 params = SchobelZhuParams{Float64}(sqrt(v0), κ, sqrt(θ), ρ, σ)
 cf = DefaultCharFunc(params)
 #with the original formula from SZ, there is a discontinuity in the log.
 sz = CharFuncPricing.evaluateLogCharFunc( cf, 2.0+0.0im, 15.0)
 lk = CharFuncPricing.evaluateLogCharFuncLK( cf, 2.0+0.0im, 15.0)
 @test isapprox(lk,sz, atol=1e-13)
 for u = 0.1:0.1:4.0
      lk = CharFuncPricing.evaluateLogCharFuncLK(cf, u+0.0im, 15.0)
      sz=CharFuncPricing.evaluateLogCharFunc( cf, u+0.0im, 15.0)
      @test isapprox(lk,sz, atol=1e-13)
  end

  v0 = 0.010201
  ρ  = -0.7
  σ = 0.61
  κ = 6.21
  θ = 0.019
  α = 3.35861
  params = SchobelZhuParams{Float64}(sqrt(v0), κ, sqrt(θ), ρ, σ)
  cf = DefaultCharFunc(params)
   #with the original formula from SZ, there is a discontinuity in the log.
  sz = CharFuncPricing.evaluateLogCharFunc(cf, 10.0+α*1im, 10.0)
  lk = CharFuncPricing.evaluateLogCharFuncLK(cf, 10.0+α*1im, 10.0)
  @test isapprox(lk,sz, atol=1e-13)

end

@testset "SZCumulantsJoshi" begin
    κ = 4.0
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.1
    τ = 0.01
    params = SchobelZhuParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc{SchobelZhuParams{Float64}, Taylor1{Complex}}(params)
    t = Taylor1(Float64, 4)
    cft = evaluateLogCharFunc(cf, t, τ)
    c1, c2, c4 = computeCumulants(params, τ)
    println(c1, " ", c2 / 2, " ", c4 / (2 * 3 * 4), " ", cft)
    @test isapprox(imag(cft.coeffs[2]), c1, atol = 1e-12)
    @test isapprox(-real(cft.coeffs[3]) * 2, c2, atol = 1e-12)
    @test isapprox(cft.coeffs[5] * 2 * 3 * 4, c4, atol = 1e-16)
end

@testset "SchobelZhuPaper" begin
    #reference values from "Stochastic Volatility with an Ornstein-Uhlenbeck Process: An Extension" (1998)
    κ = 4.0
    θ = 0.2
    σ = 0.1
    ρ = -0.5
    v0 = 0.2
    τ = 0.5
    spot = 100.0
    r = 0.0953
    forward = spot * exp(r * τ)
    df = exp(-r * τ)
    params = SchobelZhuParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    m = 128
    l = 6
    pricer =
        makeCosCharFuncPricer(cf, τ, m, l)
    strikes = [90.0, 95.0, 100.0, 105.0, 110.0, 115.0, 120.0]
    refPrices = [15.29, 11.50, 8.24, 5.60, 3.58, 2.16, 1.22]
    for (strike, refPrice) in zip(strikes, refPrices)
        price = priceEuropean(pricer, true, strike, forward, τ, df)
        println(strike, " ", price)
        @test isapprox(refPrice, price, atol = 0.01)
    end
    θ = 0.0001
    v0 = 0.15
    params = SchobelZhuParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)

    pricer =
        makeCosCharFuncPricer(cf, τ, m, l)
    refPrices = [14.22, 9.60, 5.37, 2.15, 0.50, 0.06, 0.004]
    for (strike, refPrice) in zip(strikes, refPrices)
        price = priceEuropean(pricer, true, strike, forward,τ, df)
        println(strike, " ", price)
        @test isapprox(refPrice, price, atol = 0.01)
    end

end


text:
abstract type AbstractLShapedSolver <: AbstractStructuredModel end

nscenarios(lshaped::AbstractLShapedSolver) = lshaped.nscenarios
ncuts(lshaped::AbstractLShapedSolver) = lshaped.data.ncuts
niterations(lshaped::AbstractLShapedSolver) = lshaped.data.iterations
tolerance(lshaped::AbstractLShapedSolver) = lshaped.parameters.τ

# Initialization #
# ======================================================================== #
function initialize!(lshaped::AbstractLShapedSolver, subsolver::SubSolver)
    # Initialize progress meter
    lshaped.progress.thresh = lshaped.parameters.τ
    # Initialize subproblems
    initialize_subproblems!(lshaped, scenarioproblems(lshaped.stochasticprogram), decision(lshaped), subsolver)
    # Prepare the master optimization problem
    prepare_master!(lshaped)
    # Initialize feasibility handler
    initialize_feasibility!(lshaped)
    # Initialize regularization policy
    initialize_regularization!(lshaped)
    # Finish initialization
    finish_initilization!(lshaped)
    return nothing
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function set_params!(lshaped::AbstractLShapedSolver; kwargs...)
    for (k,v) in kwargs
        setfield!(lshaped.parameters,k,v)
    end
    return nothing
end

function prepare_master!(lshaped::AbstractLShapedSolver)
    # θs
    for i = 1:nthetas(lshaped)
        MPB.addvar!(lshaped.mastersolver.lqmodel ,-1e10, Inf, initial_theta_coefficient(lshaped.feasibility))
        push!(lshaped.mastervector, -1e10)
    end
    return nothing
end

function update_solution!(lshaped::AbstractLShapedSolver, solver::LQSolver)
    ncols = decision_length(lshaped.stochasticprogram)
    nb = nthetas(lshaped)
    x = copy(getsolution(solver))
    lshaped.mastervector .= x[1:ncols+nb]
    lshaped.x[1:ncols] .= x[1:ncols]
    set_model_objectives(lshaped, x[ncols+1:ncols+nb])
    return nothing
end

function get_obj(lshaped::AbstractLShapedSolver)
    first_stage = StochasticPrograms.get_stage_one(lshaped.stochasticprogram)
    c = JuMP.prepAffObjective(first_stage)
    return first_stage.objSense == :Min ? c : -c
end

function current_objective_value(lshaped::AbstractLShapedSolver)
    return get_obj(lshaped)⋅current_decision(lshaped) + sum(subobjectives(lshaped))
end

function calculate_estimate(lshaped::AbstractLShapedSolver)
    return get_obj(lshaped)⋅lshaped.x + sum(model_objectives(lshaped))
end

function log!(lshaped::AbstractLShapedSolver; optimal = false)
    @unpack Q,θ,iterations = lshaped.data
    @unpack keep, offset, indent = lshaped.parameters
    push!(lshaped.Q_history, Q)
    push!(lshaped.θ_history, θ)
    lshaped.data.iterations += 1

    log_regularization!(lshaped)
    if lshaped.parameters.log
        current_gap = optimal ? 0.0 : gap(lshaped)
        ProgressMeter.update!(lshaped.progress,current_gap,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", objective(lshaped)),
                                  ("$(indentstr(indent))Gap", current_gap),
                                  ("$(indentstr(indent))Number of cuts", ncuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
    end
    return nothing
end

function log!(lshaped::AbstractLShapedSolver, t::Integer; optimal = false)
    @unpack Q,θ,iterations = lshaped.data
    @unpack keep, offset, indent = lshaped.parameters
    lshaped.Q_history[t] = Q
    lshaped.θ_history[t] = θ

    log_regularization!(lshaped,t)

    if lshaped.parameters.log
        current_gap = optimal ? 0.0 : gap(lshaped)
        ProgressMeter.update!(lshaped.progress, current_gap,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", objective(lshaped)),
                                  ("$(indentstr(indent))Gap", current_gap),
                                  ("$(indentstr(indent))Number of cuts", ncuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
    end
    return nothing
end

function indentstr(n::Integer)
    return repeat(" ", n)
end

function check_optimality(lshaped::AbstractLShapedSolver)
    @unpack τ = lshaped.parameters
    @unpack θ = lshaped.data
    return θ > -Inf && gap(lshaped) <= τ
end
# ======================================================================== #

# Cut functions #
# ======================================================================== #
active(lshaped::AbstractLShapedSolver, hyperplane::AbstractHyperPlane) = active(hyperplane, decision(lshaped), tolerance(lshaped))
active(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = optimal(cut, decision(lshaped), model_objectives(lshaped)[cut.id], tolerance(lshaped))
active(lshaped::AbstractLShapedSolver, cut::AggregatedOptimalityCut) = optimal(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]), tolerance(lshaped))
satisfied(lshaped::AbstractLShapedSolver, hyperplane::AbstractHyperPlane) = satisfied(hyperplane, decision(lshaped), tolerance(lshaped))
satisfied(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = satisfied(cut, decision(lshaped), model_objectives(lshaped)[cut.id], tolerance(lshaped))
satisfied(lshaped::AbstractLShapedSolver, cut::AggregatedOptimalityCut) = satisfied(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]), tolerance(lshaped))
violated(lshaped::AbstractLShapedSolver, hyperplane::AbstractHyperPlane) = !satisfied(lshaped, hyperplane)
gap(lshaped::AbstractLShapedSolver, hyperplane::AbstractHyperPlane) = gap(hyperplane, decision(lshaped))
gap(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = gap(cut, decision(lshaped), model_objectives(lshaped)[cut.id])
gap(lshaped::AbstractLShapedSolver, cut::AggregatedOptimalityCut) = gap(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]))

function add_cut!(lshaped::AbstractLShapedSolver, cut::AbstractHyperPlane; consider_consolidation = true, check = true)
    added = add_cut!(lshaped, cut, model_objectives(lshaped), subobjectives(lshaped), check = check)
    added && update_objective!(lshaped, cut)
    if consider_consolidation
        added && add_cut!(lshaped, lshaped.consolidation, cut)
    end
    return added
end

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}, θs::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    θ = θs[cut.id]
    @unpack τ, cut_scaling = lshaped.parameters
    # Update objective
    subobjectives[cut.id] = Q
    # Check if cut gives new information
    if check && θ > -Inf && (θ + τ >= Q || θ + τ >= cut(lshaped.x))
        # Optimal with respect to this subproblem
        return false
    end
    # Add optimality cut
    process_cut!(lshaped, cut)
    MPB.addconstr!(lshaped.mastersolver.lqmodel, lowlevel(cut, cut_scaling)...)
    lshaped.data.ncuts += 1
    if lshaped.parameters.debug
        push!(lshaped.cuts, cut)
    end
    return true
end
function add_cut!(lshaped::AbstractLShapedSolver, cut::AggregatedOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    θs = θs[cut.ids]
    θ = sum(θs)
    @unpack τ, cut_scaling = lshaped.parameters
    # Update objective
    subobjectives[cut.ids] .= Q/length(cut.ids)
    # Check if cut gives new information
    if check && θ > -Inf && (θ + τ >= Q || θ + τ >= cut(lshaped.x))
        # Optimal with respect to these subproblems
        return false
    end
    # Add optimality cut
    process_cut!(lshaped, cut)
    MPB.addconstr!(lshaped.mastersolver.lqmodel, lowlevel(cut, cut_scaling)...)
    lshaped.data.ncuts += 1
    if lshaped.parameters.debug
        push!(lshaped.cuts, cut)
    end
    return true
end
add_cut!(lshaped::AbstractLShapedSolver, cut::AnyOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector, x::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, cut(x); check = check)
add_cut!(lshaped::AbstractLShapedSolver, cut::AnyOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, current_decision(lshaped); check = check)

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{FeasibilityCut}, ::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    # Ensure that there is no false convergence
    subobjectives[cut.id] = Q
    # Add feasibility cut
    process_cut!(lshaped, cut)
    MPB.addconstr!(lshaped.mastersolver.lqmodel, lowlevel(cut)...)
    lshaped.data.ncuts += 1
    if lshaped.parameters.debug
        push!(lshaped.feasibility.cuts, cut)
    end
    return true
end
add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{FeasibilityCut}, θs::AbstractVector, subobjectives::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, Inf)

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{Infeasible}, ::AbstractVector, subobjectives::AbstractVector; check = true)
    subobjectives[cut.id] = Inf
    return true
end

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{Unbounded}, ::AbstractVector, subobjectives::AbstractVector; check = true)
    subobjectives[cut.id] = -Inf
    return true
end


text:
struct EvalFunc{F} <: Function
  f::F
end
(f::EvalFunc)(args...) = f.f(args...)

NO_TSPAN_PROBS = Union{AbstractLinearProblem,AbstractNonlinearProblem,
                       AbstractIntegralProblem,AbstractSteadyStateProblem,
                       AbstractJumpProblem}

has_kwargs(_prob::DEProblem) = has_kwargs(typeof(_prob))
Base.@pure has_kwargs(::Type{T}) where {T} = :kwargs ∈ fieldnames(T)

const allowedkeywords = (
  :dense,
  :saveat,
  :save_idxs,
  :tstops,
  :d_discontinuities,
  :save_everystep,
  :save_on,
  :save_start,
  :save_end,
  :initialize_save,
  :adaptive,
  :abstol,
  :reltol,
  :dt,
  :dtmax,
  :dtmin,
  :force_dtmin,
  :internalnorm,
  :controller,
  :gamma,
  :beta1,
  :beta2,
  :qmax,
  :qmin,
  :qsteady_min,
  :qsteady_max,
  :qoldinit,
  :failfactor,
  :calck,
  :alias_u0,
  :maxiters,
  :callback,
  :isoutofdomain,
  :unstable_check,
  :verbose,
  :merge_callbacks,
  :progress,
  :progress_steps,
  :progress_name,
  :progress_message,
  :timeseries_errors,
  :dense_errors,
  :weak_timeseries_errors,
  :weak_dense_errors,
  :calculate_errors,
  :initializealg,
  :alg,
  :save_noise,
  :delta,
  :seed,
  :alg_hints,
  :kwargshandle,
  :trajectories,
  :batch_size,
  :sensealg,
  :advance_to_tstop,
  :stop_at_next_tstop,
  # These two are from the default algorithm handling
  :default_set,
  :second_time
)

const KWARGWARN_MESSAGE =
"""
Unrecognized keyword arguments found. Future versions will error.
The only allowed keyword arguments to `solve` are:
$allowedkeywords

See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.

Set kwargshandle=KeywordArgError for an error message.
Set kwargshandle=KeywordArgSilent to ignore this message.
"""

const KWARGERROR_MESSAGE =
"""
Unrecognized keyword arguments found.
The only allowed keyword arguments to `solve` are:
$allowedkeywords

See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.
"""

struct CommonKwargError <: Exception
  kwargs
end

function Base.showerror(io::IO, e::CommonKwargError)
  println(io, KWARGERROR_MESSAGE)
  notin = collect(map(x -> x ∉ allowedkeywords, keys(e.kwargs)))
  unrecognized = collect(keys(e.kwargs))[notin]
  print(io, "Unrecognized keyword arguments: ")
  printstyled(io, unrecognized; bold=true, color=:red)
  print(io, "\n\n")
end

@enum KeywordArgError KeywordArgWarn KeywordArgSilent

const INCOMPATIBLE_U0_MESSAGE =
"""
Initial condition incompatible with functional form.
Detected an in-place function with an initial condition of type Number or SArray.
This is incompatible because Numbers cannot be mutated, i.e.
`x = 2.0; y = 2.0; x .= y` will error.

If using a immutable initial condition type, please use the out-of-place form.
I.e. define the function `du=f(u,p,t)` instead of attempting to "mutate" the immutable `du`.

If your differential equation function was defined with multiple dispatches and one is
in-place, then the automatic detection will choose in-place. In this case, override the
choice in the problem constructor, i.e. `ODEProblem{false}(f,u0,tspan,p,kwargs...)`.

For a longer discussion on mutability vs immutability and in-place vs out-of-place, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation
"""

struct IncompatibleInitialConditionError <: Exception end

function Base.showerror(io::IO, e::IncompatibleInitialConditionError)
  print(io, INCOMPATIBLE_U0_MESSAGE)
end

const NO_DEFAULT_ALGORITHM_MESSAGE =
"""
Default algorithm choices require DifferentialEquations.jl.
Please specify an algorithm (e.g., `solve(prob, Tsit5())` for an ODE)
or import DifferentialEquations directly.

You can find the list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/
and its associated pages.
"""

struct NoDefaultAlgorithmError <: Exception end

function Base.showerror(io::IO, e::NoDefaultAlgorithmError)
  print(io, NO_DEFAULT_ALGORITHM_MESSAGE)
end

const NO_TSPAN_MESSAGE =
"""
No tspan is set in the problem or chosen in the init/solve call
"""

struct NoTspanError <: Exception end

function Base.showerror(io::IO, e::NoTspanError)
  print(io, NO_TSPAN_MESSAGE)
end

const NON_SOLVER_MESSAGE =
"""
The arguments to solve are incorrect.
The second argument must be a solver choice, `solve(prob,alg)`
where `alg` is a `<: DEAlgorithm`, i.e. `Tsti5()`.

Please double check the arguments being sent to the solver.
"""

struct NonSolverError <: Exception end

function Base.showerror(io::IO, e::NonSolverError)
  print(io, NON_SOLVER_MESSAGE)
end

const PROBSOLVER_PAIRING_MESSAGE =
"""
Incompatible problem+solver pairing.
For example, this can occur if an ODE solver is passed with an SDEProblem.
Solvers are only capable of handling specific problem types. Please double
check that the chosen pairing is capable for handling the given problems.
"""

struct ProblemSolverPairingError <: Exception
  prob
  alg
end

function Base.showerror(io::IO, e::ProblemSolverPairingError)
  println(io, PROBSOLVER_PAIRING_MESSAGE)
  println(io, "Problem type: $(SciMLBase.__parameterless_type(typeof(e.prob)))")
  println(io, "Solver type: $(SciMLBase.__parameterless_type(typeof(e.alg)))")
  println(io, "Problem types compatible with the chosen solver: $(compatible_problem_types(e.prob,e.alg))")
end

function compatible_problem_types(prob,alg)
  if alg isa AbstractODEAlgorithm
    ODEProblem
  elseif alg isa AbstractSDEAlgorithm
    (SDEProblem, SDDEProblem)
  elseif alg isa AbstractDDEAlgorithm # StochasticDelayDiffEq.jl just uses the SDE alg
    DDEProblem
  elseif alg isa AbstractDAEAlgorithm
    DAEProblem
  elseif alg isa AbstractSteadyStateAlgorithm
    SteadyStateProblem
  end
end

const DIRECT_AUTODIFF_INCOMPATABILITY_MESSAGE =
"""
Incompatible solver + automatic differentiation pairing.
The chosen automatic differentiation algorithm requires the ability
for compiler transforms on the code which is only possible on pure-Julia
solvers such as those from OrdinaryDiffEq.jl. Direct differentiation methods
which require this ability include:

- Direct use of ForwardDiff.jl on the solver
- `ForwardDiffSensitivity`, `ReverseDiffAdjoint`, `TrackerAdjoint`, and `ZygoteAdjoint`
  sensealg choices for adjoint differentiation.

Either switch the choice of solver to a pure Julia method, or change the automatic
differentiation method to one that does not require such transformations.

For more details on automatic differentiation, adjoint, and sensitivity analysis
of differential equations, see the documentation page:

https://diffeq.sciml.ai/stable/analysis/sensitivity/
"""

struct DirectAutodiffError <: Exception end

function Base.showerror(io::IO, e::DirectAutodiffError)
  println(io, DIRECT_AUTODIFF_INCOMPATABILITY_MESSAGE)
end

const NONCONCRETE_ELTYPE_MESSAGE =
"""
Non-concrete element type inside of an `Array` detected.
Arrays with non-concrete element types, such as
`Array{Union{Float32,Float64}}`, are not supported by the
differential equation solvers. Anyways, this is bad for
performance so you don't want to be doing this!

If this was a mistake, promote the element types to be
all the same. If this was intentional, for example,
using Unitful.jl with different unit values, then use
an array type which has fast broadcast support for
heterogeneous values such as the ArrayPartition
from RecursiveArrayTools.jl. For example:

```julia
using RecursiveArrayTools
x = ArrayPartition([1.0,2.0],[1f0,2f0])
y = ArrayPartition([3.0,4.0],[3f0,4f0])
x .+ y # fast, stable, and usable as u0 into DiffEq!
```

Element type:
"""

struct NonConcreteEltypeError <: Exception
  eltype
end

function Base.showerror(io::IO, e::NonConcreteEltypeError)
  print(io, NONCONCRETE_ELTYPE_MESSAGE)
  print(io,e.eltype)
end

const GENERIC_NUMBER_TYPE_ERROR_MESSAGE =
"""
Non-standard number type (i.e. not Float32, Float64,
ComplexF32, or ComplexF64) detected as the element type
for the initial condition or time span. These generic
number types are only compatible with the pure Julia
solvers which support generic programming, such as
OrdinaryDiffEq.jl. The chosen solver does not support
this functionality. Please double check that the initial
condition and time span types are correct, and check that
the chosen solver was correct.
"""

struct GenericNumberTypeError <: Exception
  alg
  uType
  tType
end

function Base.showerror(io::IO, e::GenericNumberTypeError)
  println(io, GENERIC_NUMBER_TYPE_ERROR_MESSAGE)
  println(io, "Solver: $(e.alg)")
  println(io, "u0 type: $(e.uType)")
  print(io,"Timespan type: $(e.tType)")
end

const COMPLEX_SUPPORT_ERROR_MESSAGE =
"""
Complex number type (i.e. ComplexF32, or ComplexF64)
detected as the element type for the initial condition
with an algorithm that does not support complex numbers.
Please check that the initial condition type is correct.
If complex number support is needed, try different solvers
such as those from OrdinaryDiffEq.jl.
"""

struct ComplexSupportError <: Exception
  alg
end

function Base.showerror(io::IO, e::ComplexSupportError)
  println(io, COMPLEX_SUPPORT_ERROR_MESSAGE)
  println(io, "Solver: $(e.alg)")
end

const COMPLEX_TSPAN_ERROR_MESSAGE =
"""
Complex number type (i.e. ComplexF32, or ComplexF64)
detected as the element type for the independent variable
(i.e. time span). Please check that the tspan type is correct.
No solvers support complex time spans. If this is required,
please open an issue.
"""

struct ComplexTspanError <: Exception end

function Base.showerror(io::IO, e::ComplexTspanError)
  println(io, COMPLEX_TSPAN_ERROR_MESSAGE)
end

const TUPLE_STATE_ERROR_MESSAGE =
"""
Tuple type used as a state. Since a tuple does not have vector
properties, it will not work as a state type in equation solvers.
Instead, change your equation from using tuple constructors `()`
to static array constructors `SA[]`. For example, change:

```julia
function ftup((a,b),p,t)
  return b,-a
end
u0 = (1.0,2.0)
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)
```

to:

```julia
using StaticArrays
function fsa(u,p,t)
    SA[u[2],u[1]]
end
u0 = SA[1.0,2.0]
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)
```

This will be safer and fast for small ODEs. For more information, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Further-Optimizations-of-Small-Non-Stiff-ODEs-with-StaticArrays
"""

struct TupleStateError <: Exception end

function Base.showerror(io::IO, e::TupleStateError)
  println(io, TUPLE_STATE_ERROR_MESSAGE)
end

function init_call(_prob, args...; merge_callbacks=true, kwargshandle=KeywordArgWarn, kwargs...)

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  checkkwargs(kwargshandle; kwargs...)

  if hasfield(typeof(_prob), :f) && hasfield(typeof(_prob.f), :f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__init, _prob, args...; kwargs...)#::T
  else
    __init(_prob, args...; kwargs...)#::T
  end
end

function init(prob::DEProblem, args...; kwargs...)
  if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _prob = get_concrete_problem(prob, isadaptive(alg); kwargs...)
    check_prob_alg_pairing(_prob, alg)
    init_call(_prob, alg, args...; kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _prob = get_concrete_problem(prob, isadaptive(alg); kwargs...)
    check_prob_alg_pairing(_prob, alg)
    init_call(_prob, args...; kwargs...)
  else
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); kwargs...)
    init_call(_prob, args...; kwargs...)
  end
end

function solve_call(_prob, args...; merge_callbacks=true, kwargshandle=KeywordArgWarn, kwargs...)
  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  checkkwargs(kwargshandle; kwargs...)
  isdefined(_prob, :u0) && _prob.u0 isa Array &&
                          !isconcretetype(RecursiveArrayTools.recursive_unitless_eltype(_prob.u0)) &&
                          throw(NonConcreteEltypeError(RecursiveArrayTools.recursive_unitless_eltype(_prob.u0)))

  if hasfield(typeof(_prob), :f) && hasfield(typeof(_prob.f), :f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__solve, _prob, args...; kwargs...)#::T
  else
    __solve(_prob, args...; kwargs...)#::T
  end
end

"""
solve(prob::DEProblem, alg::Union{DEAlgorithm,Nothing}; kwargs...)

## Arguments

The only positional argument is `alg` which is optional. By default, `alg = nothing`.
If `alg = nothing`, then `solve` dispatches to the DifferentialEquations.jl automated
algorithm selection (if `using DifferentialEquations` was done, otherwise it will
error with a MethodError`).

## Keyword Arguments

The DifferentialEquations.jl universe has a large set of common arguments available
for the `solve` function. These arguments apply to `solve` on any problem type and
are only limited by limitations of the specific implementations.

Many of the defaults depend on the algorithm or the package the algorithm derives
from. Not all of the interface is provided by every algorithm.
For more detailed information on the defaults and the available options
for specific algorithms / packages, see the manual pages for the solvers of specific
problems. To see whether a specific package is compaible with the use of a
given option, see the [Solver Compatibility Chart](@ref)

### Default Algorithm Hinting

To help choose the default algorithm, the keyword argument `alg_hints` is
provided to `solve`. `alg_hints` is a `Vector{Symbol}` which describe the
problem at a high level to the solver. The options are:

* `:auto` vs `:nonstiff` vs `:stiff` - Denotes the equation as nonstiff/stiff.
  `:auto` allow the default handling algorithm to choose stiffness detection
  algorithms. The default handling defaults to using `:auto`.

Currently unused options include:

* `:interpolant` - Denotes that a high-precision interpolation is important.
* `:memorybound` - Denotes that the solver will be memory bound.

This functionality is derived via the benchmarks in
[SciMLBenchmarks.jl](https://github.com/SciML/SciMLBenchmarks.jl)

#### SDE Specific Alghints

* `:additive` - Denotes that the underlying SDE has additive noise.
* `:stratonovich` - Denotes that the solution should adhere to the Stratonovich
  interpretation.

### Output Control

These arguments control the output behavior of the solvers. It defaults to maximum
output to give the best interactive user experience, but can be reduced all the
way to only saving the solution at the final timepoint.

The following options are all related to output control. See the "Examples"
section at the end of this page for some example usage.

* `dense`: Denotes whether to save the extra pieces required for dense (continuous)
  output. Default is `save_everystep && !isempty(saveat)` for algorithms which have
  the ability to produce dense output, i.e. by default it's `true` unless the user
  has turned off saving on steps or has chosen a `saveat` value. If `dense=false`,
  the solution still acts like a function, and `sol(t)` is a linear interpolation
  between the saved time points.
* `saveat`: Denotes specific times to save the solution at, during the solving
  phase. The solver will save at each of the timepoints in this array in the
  most efficient manner available to the solver. If only `saveat` is given, then
  the arguments `save_everystep` and `dense` are `false` by default.
  If `saveat` is given a number, then it will automatically expand to
  `tspan[1]:saveat:tspan[2]`. For methods where interpolation is not possible,
  `saveat` may be equivalent to `tstops`. The default value is `[]`.
* `save_idxs`: Denotes the indices for the components of the equation to save.
  Defaults to saving all indices. For example, if you are solving a 3-dimensional ODE,
  and given `save_idxs = [1, 3]`, only the first and third components of the
  solution will be outputted.
  Notice that of course in this case the outputed solution will be two-dimensional.
* `tstops`: Denotes *extra* times that the timestepping algorithm must step to.
  This should be used to help the solver deal with discontinuities and
  singularities, since stepping exactly at the time of the discontinuity will
  improve accuracy. If a method cannot change timesteps (fixed timestep
  multistep methods), then `tstops` will use an interpolation,
  matching the behavior of `saveat`. If a method cannot change timesteps and
  also cannot interpolate, then `tstops` must be a multiple of `dt` or else an
  error will be thrown. Default is `[]`.
* `d_discontinuities:` Denotes locations of discontinuities in low order derivatives.
  This will force FSAL algorithms which assume derivative continuity to re-evaluate
  the derivatives at the point of discontinuity. The default is `[]`.
* `save_everystep`: Saves the result at every step.
  Default is true if `isempty(saveat)`.
* `save_on`: Denotes whether intermediate solutions are saved. This overrides the
  settings of `dense`, `saveat` and `save_everystep` and is used by some applicatioins
  to manually turn off saving temporarily. Everyday use of the solvers should leave
  this unchanged. Defaults to `true`.
* `save_start`: Denotes whether the initial condition should be included in
  the solution type as the first timepoint. Defaults to `true`.
* `save_end`: Denotes whether the final timepoint is forced to be saved,
  regardless of the other saving settings. Defaults to `true`.
* `initialize_save`: Denotes whether to save after the callback initialization
  phase (when `u_modified=true`). Defaults to `true`.

Note that `dense` requires `save_everystep=true` and `saveat=false`. If you need
additional saving while keeping dense output, see
[the SavingCallback in the Callback Library](@ref saving_callback).

### Stepsize Control

These arguments control the timestepping routines.

#### Basic Stepsize Control

These are the standard options for controlling stepping behavior. Error estimates
do the comparison

```math
err_{scaled} = err/(abstol + max(uprev,u)*reltol)
```

The scaled error is guaranteed to be `<1` for a given local error estimate
(note: error estimates are local unless the method specifies otherwise). `abstol`
controls the non-scaling error and thus can be thought of as the error around zero.
`reltol` scales with the size of the dependent variables and so one can interpret
`reltol=1e-3` as roughly being (locally) correct to 3 digits. Note tolerances can
be specified element-wise by passing a vector whose size matches `u0`.

* `adaptive`: Turns on adaptive timestepping for appropriate methods. Default
  is true.
* `abstol`: Absolute tolerance in adaptive timestepping. This is the tolerance
  on local error estimates, not necessarily the global error (though these quantities
  are related). Defaults to `1e-6` on deterministic equations (ODEs/DDEs/DAEs) and `1e-2`
  on stochastic equations (SDEs/RODEs).
* `reltol`: Relative tolerance in adaptive timestepping.  This is the tolerance
  on local error estimates, not necessarily the global error (though these quantities
  are related). Defaults to `1e-3` on deterministic equations (ODEs/DDEs/DAEs) and `1e-2`
  on stochastic equations (SDEs/RODEs).
* `dt`: Sets the initial stepsize. This is also the stepsize for fixed
  timestep methods. Defaults to an automatic choice if the method is adaptive.
* `dtmax`: Maximum dt for adaptive timestepping. Defaults are
  package-dependent.
* `dtmin`: Minimum dt for adaptive timestepping. Defaults are
  package-dependent.
* `force_dtmin`: Declares whether to continue, forcing the minimum `dt` usage.
  Default is `false`, which has the solver throw a warning and exit early when
  encountering the minimum `dt`. Setting this true allows the solver to continue,
  never letting `dt` go below `dtmin` (and ignoring error tolerances in those
  cases). Note that `true` is not compatible with most interop packages.

#### Fixed Stepsize Usage

Note that if a method does not have adaptivity, the following rules apply:

* If `dt` is set, then the algorithm will step with size `dt` each iteration.
* If `tstops` and `dt` are both set, then the algorithm will step with either a
  size `dt`, or use a smaller step to hit the `tstops` point.
* If `tstops` is set without `dt`, then the algorithm will step directly to
  each value in `tstops`
* If neither `dt` nor `tstops` are set, the solver will throw an error.

#### [Advanced Adaptive Stepsize Control](@id advanced_adaptive_stepsize_control)

These arguments control more advanced parts of the internals of adaptive timestepping
and are mostly used to make it more efficient on specific problems. For detained
explanations of the timestepping algorithms, see the
[timestepping descriptions](@ref timestepping)

* `internalnorm`: The norm function `internalnorm(u,t)` which error estimates
  are calculated. Required are two dispatches: one dispatch for the state variable
  and the other on the elements of the state variable (scalar norm).
  Defaults are package-dependent.
* `controller`: Possible examples are [`IController`](@ref),
  [`PIController`](@ref), [`PIDController`](@ref), [`PredictiveController`](@ref).
  Default is algorithm-dependent.
* `gamma`: The risk-factor γ in the q equation for adaptive timestepping
  of the controllers using it.
  Default is algorithm-dependent.
* `beta1`: The Lund stabilization α parameter.
  Default is algorithm-dependent.
* `beta2`: The Lund stabilization β parameter.
  Default is algorithm-dependent.
* `qmax`: Defines the maximum value possible for the adaptive q.
  Default is algorithm-dependent.
* `qmin`: Defines the minimum value possible for the adaptive q.
  Default is algorithm-dependent.
* `qsteady_min`: Defines the minimum for the range around 1 where the timestep
  is held constant. Default is algorithm-dependent.
* `qsteady_max`: Defines the maximum for the range around 1 where the timestep
  is held constant. Default is algorithm-dependent.
* `qoldinit`: The initial `qold` in stabilization stepping.
  Default is algorithm-dependent.
* `failfactor`: The amount to decrease the timestep by if the Newton iterations
  of an implicit method fail. Default is 2.

### Memory Optimizations

* `calck`: Turns on and off the internal ability for intermediate
  interpolations (also known as intermediate density). Not the same as `dense`, which is post-solution interpolation.
  This defaults to `dense || !isempty(saveat) ||  "no custom callback is given"`.
  This can be used to turn off interpolations
  (to save memory) if one isn't using interpolations when a custom callback is
  used. Another case where this may be used is to turn on interpolations for
  usage in the integrator interface even when interpolations are used nowhere else.
  Note that this is only required if the algorithm doesn't have
  a free or lazy interpolation (`DP8()`). If `calck = false`, `saveat` cannot be used.
  The rare keyword `calck` can be useful in event handling.
* `alias_u0`: allows the solver to alias the initial condition array that is contained
  in the problem struct. Defaults to false.

### Miscellaneous

* `maxiters`: Maximum number of iterations before stopping. Defaults to 1e5.
* `callback`: Specifies a callback. Defaults to a callback function which
  performs the saving routine. For more information, see the
  [Event Handling and Callback Functions manual page](@ref callbacks).
* `isoutofdomain`: Specifies a function `isoutofdomain(u,p,t)` where, when it
  returns true, it will reject the timestep. Disabled by default.
* `unstable_check`: Specifies a function `unstable_check(dt,u,p,t)` where, when
  it returns true, it will cause the solver to exit and throw a warning. Defaults
  to `any(isnan,u)`, i.e. checking if any value is a NaN.
* `verbose`: Toggles whether warnings are thrown when the solver exits early.
  Defaults to true.
* `merge_callbacks`: Toggles whether to merge `prob.callback` with the `solve` keyword
  argument `callback`. Defaults to `true`.

### Progress Monitoring

These arguments control the usage of the progressbar in the Juno IDE.

* `progress`: Turns on/off the Juno progressbar. Default is false.
* `progress_steps`: Numbers of steps between updates of the progress bar.
  Default is 1000.
* `progress_name`: Controls the name of the progressbar. Default is the name
  of the problem type.
* `progress_message`: Controls the message with the progressbar. Defaults to
  showing `dt`, `t`, the maximum of `u`.

### Error Calculations

If you are using the test problems (ex: `ODETestProblem`), then the following
options control the errors which are calculated:

* `timeseries_errors`: Turns on and off the calculation of errors at the steps
  which were taken, such as the `l2` error. Default is true.
* `dense_errors`: Turns on and off the calculation of errors at the steps which
  require dense output and calculate the error at 100 evenly-spaced points
  throughout `tspan`. An example is the `L2` error. Default is false.

### Sensitivity Algorithms (`sensealg`)

`sensealg` is used for

## Examples

The following lines are examples of how one could use the configuration of
`solve()`. For these examples a 3-dimensional ODE problem is assumed, however
the extention to other types is straightforward.

1. `solve(prob, AlgorithmName())` : The "default" setting, with a user-specified
  algorithm (given by `AlgorithmName()`). All parameters get their default values.
  This means that the solution is saved at the steps the Algorithm stops internally
  and dense output is enabled if the chosen algorithm allows for it.

  All other integration parameters (e.g. stepsize) are chosen automatically.
2. `solve(prob, saveat = 0.01, abstol = 1e-9, reltol = 1e-9)` : Standard setting
  for accurate output at specified (and equidistant) time intervals, used for
  e.g. Fourier Transform. The solution is given every 0.01 time units,
  starting from `tspan[1]`. The solver used is `Tsit5()` since no keyword
  `alg_hits` is given.

3. `solve(prob, maxiters = 1e7, progress = true, save_idxs = [1])` : Using longer
  maximum number of solver iterations can be useful when a given `tspan` is very
  long. This example only saves the first of the variables of the system, either
  to save size or because the user does not care about the others. Finally, with
  `progress = true` you are enabling the progress bar.
"""
function solve(prob::DEProblem, args...; sensealg=nothing,
  u0=nothing, p=nothing, kwargs...)
  u0 = u0 !== nothing ? u0 : prob.u0
  p = p !== nothing ? p : prob.p
  if sensealg === nothing && haskey(prob.kwargs, :sensealg)
    sensealg = prob.kwargs[:sensealg]
  end
  solve_up(prob, sensealg, u0, p, args...; kwargs...)
end

function solve_up(prob::DEProblem, sensealg, u0, p, args...; kwargs...)

  if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _alg = prepare_alg(alg, u0, p, prob)
    _prob = get_concrete_problem(prob, isadaptive(_alg); u0=u0, p=p, kwargs...)
    check_prob_alg_pairing(_prob, alg) # use alg for improved inference
    solve_call(_prob, _alg, args...; kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _alg = prepare_alg(alg, u0, p, prob)
    _prob = get_concrete_problem(prob, isadaptive(_alg); u0=u0, p=p, kwargs...)
    check_prob_alg_pairing(_prob, alg)
    solve_call(_prob, _alg, Base.tail(args)...; kwargs...)
  elseif isempty(args) # Default algorithm handling
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
    solve_call(_prob, args...; kwargs...)
  else
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
    solve_call(_prob, args...; kwargs...)
  end
end

function solve(prob::EnsembleProblem, args...; kwargs...)
  if isempty(args) || length(args) == 1 && typeof(args[1]) <: EnsembleAlgorithm
    __solve(prob, nothing, args...; kwargs...)
  else
    __solve(prob, args...; kwargs...)
  end
end

function solve(prob::AbstractNoiseProblem, args...; kwargs...)
  __solve(prob, args...; kwargs...)
end

function solve(prob::AbstractJumpProblem, args...; kwargs...)
  __solve(prob, args...; kwargs...)
end

function checkkwargs(kwargshandle; kwargs...)
  if any(x -> x ∉ allowedkeywords, keys(kwargs))
    if kwargshandle == KeywordArgError
      throw(CommonKwargError(kwargs))
    elseif kwargshandle == KeywordArgWarn
      @warn KWARGWARN_MESSAGE
      unrecognized = setdiff(keys(kwargs), allowedkeywords)
      print("Unrecognized keyword arguments: ")
      printstyled(unrecognized; bold=true, color=:red)
      print("\n\n")
    else
      @assert kwargshandle == KeywordArgSilent
    end
  end
end

@non_differentiable checkkwargs(kwargshandle)

function get_concrete_problem(prob::AbstractJumpProblem, isadapt; kwargs...)
  prob
end

function get_concrete_problem(prob::SteadyStateProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, Inf, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0=u0)
end

function get_concrete_problem(prob::NonlinearProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, nothing, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0=u0)
end

function get_concrete_problem(prob::AbstractEnsembleProblem, isadapt; kwargs...)
  prob
end

function solve(prob::PDEProblem, alg::DiffEqBase.DEAlgorithm, args...;
  kwargs...)
  solve(prob.prob, alg, args...; kwargs...)
end

function init(prob::PDEProblem, alg::DiffEqBase.DEAlgorithm, args...;
  kwargs...)
  init(prob.prob, alg, args...; kwargs...)
end

function get_concrete_problem(prob, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  u0_promote = promote_u0(u0, p, tspan[1])
  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
     prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
     p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f=f_promote, u0=u0_promote, p=p, tspan=tspan_promote)
  end
end

function get_concrete_problem(prob::DAEProblem, isadapt; kwargs...)

  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  du0 = get_concrete_du0(prob, isadapt, tspan[1], kwargs)

  u0_promote = promote_u0(u0, p, tspan[1])
  du0_promote = promote_u0(du0, p, tspan[1])

  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
     isconcretedu0(prob, tspan[1], kwargs) && typeof(du0_promote) === typeof(prob.du0) &&
     prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
     p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f=f_promote, du0=du0_promote, u0=u0_promote, p=p, tspan=tspan_promote)
  end
end

function get_concrete_problem(prob::DDEProblem, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)

  if prob.constant_lags isa Function
    constant_lags = prob.constant_lags(p)
  else
    constant_lags = prob.constant_lags
  end

  u0 = promote_u0(u0, p, tspan[1])
  tspan = promote_tspan(u0, p, tspan, prob, kwargs)

  remake(prob; u0=u0, tspan=tspan, p=p, constant_lags=constant_lags)
end

function promote_f(f::F, u0) where {F}
  # Ensure our jacobian will be of the same type as u0
  uElType = u0 === nothing ? Float64 : eltype(u0)
  if isdefined(f, :jac_prototype) && f.jac_prototype isa AbstractArray
    f = @set f.jac_prototype = similar(f.jac_prototype, uElType)
  end
  return f
end

promote_f(f::SplitFunction, u0) = typeof(f.cache) === typeof(u0) && isinplace(f) ? f : remake(f, cache=zero(u0))
prepare_alg(alg, u0, p, f) = alg

function get_concrete_tspan(prob, isadapt, kwargs, p)
  if prob.tspan isa Function
    tspan = prob.tspan(p)
  elseif haskey(kwargs, :tspan)
    tspan = kwargs[:tspan]
  elseif prob.tspan === (nothing, nothing)
    throw(NoTspanError())
  else
    tspan = prob.tspan
  end

  isadapt && eltype(tspan) <: Integer && (tspan = float.(tspan))

  tspan
end

function isconcreteu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.u0 !== nothing && !isdistribution(prob.u0)
end

function isconcretedu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.du0 !== nothing && !isdistribution(prob.du0)
end

function get_concrete_u0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif haskey(kwargs, :u0)
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  isadapt && eltype(u0) <: Integer && (u0 = float.(u0))

  _u0 = handle_distribution_u0(u0)

  if isinplace(prob) && (_u0 isa Number || _u0 isa SArray)
    throw(IncompatibleInitialConditionError())
  end

  if _u0 isa Tuple
    throw(TupleStateError())
  end

  _u0
end

function get_concrete_du0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.du0)
    du0 = prob.du0(prob.p, t0)
  elseif haskey(kwargs, :du0)
    du0 = kwargs[:du0]
  else
    du0 = prob.du0
  end

  isadapt && eltype(du0) <: Integer && (du0 = float.(du0))

  _du0 = handle_distribution_u0(du0)

  if isinplace(prob) && (_du0 isa Number || _du0 isa SArray)
    throw(IncompatibleInitialConditionError())
  end

  _du0
end

function get_concrete_p(prob, kwargs)
  if haskey(kwargs, :p)
    p = kwargs[:p]
  else
    p = prob.p
  end
end

handle_distribution_u0(_u0) = _u0
handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
isdistribution(_u0::Distributions.Sampleable) = true

eval_u0(u0::Function) = true
eval_u0(u0) = false

function __solve(prob::DEProblem, args...; default_set=false, second_time=false, kwargs...)
  if second_time
    throw(NoDefaultAlgorithmError())
  elseif length(args) > 0 && !(typeof(args[1]) <: Union{Nothing,DEAlgorithm})
    throw(NonSolverError())
  else
    __solve(prob::DEProblem, nothing, args...; default_set=false, second_time=true, kwargs...)
  end
end

function check_prob_alg_pairing(prob, alg)
  if prob isa ODEProblem && !(alg isa AbstractODEAlgorithm) ||
     prob isa SDEProblem && !(alg isa AbstractSDEAlgorithm) ||
     prob isa SDDEProblem && !(alg isa AbstractSDEAlgorithm) ||
     prob isa DDEProblem && !(alg isa AbstractDDEAlgorithm) ||
     prob isa DAEProblem && !(alg isa AbstractDAEAlgorithm) ||
     prob isa SteadyStateProblem && !(alg isa AbstractSteadyStateAlgorithm)

     throw(ProblemSolverPairingError(prob, alg))
  end

  if isdefined(prob, :u0) && eltype(prob.u0) <: ForwardDiff.Dual && !SciMLBase.isautodifferentiable(alg)
    throw(DirectAutodiffError())
  end

  # Complex number support comes before arbitrary number support for a more direct
  # error message.
  if !SciMLBase.allowscomplex(alg)
    if isdefined(prob, :u0) && RecursiveArrayTools.recursive_unitless_eltype(prob.u0) <: Complex
        throw(ComplexSupportError(alg))
    end
  end

  if isdefined(prob, :tspan) && eltype(prob.tspan) <: Complex
    throw(ComplexTspanError())
  end

  # Check for concrete element type so that the non-concrete case throws a better error
  if !SciMLBase.allows_arbitrary_number_types(alg)
    if isdefined(prob, :u0)
      uType = RecursiveArrayTools.recursive_unitless_eltype(prob.u0)
      if Base.isconcretetype(uType) &&
        !(uType <: Union{Float32,Float64,ComplexF32,ComplexF64})

        throw(GenericNumberTypeError(alg, isdefined(prob, :u0) ? prob.u0 : nothing,
        isdefined(prob, :tspan) ? prob.tspan : nothing))
      end
    end

    if isdefined(prob, :tspan)
      tType = RecursiveArrayTools.recursive_unitless_eltype(prob.tspan)
      if Base.isconcretetype(tType) &&
        !(tType <: Union{Float32,Float64,ComplexF32,ComplexF64})

        throw(GenericNumberTypeError(alg, isdefined(prob, :u0) ? prob.u0 : nothing,
        isdefined(prob, :tspan) ? prob.tspan : nothing))
      end
    end
  end

end

################### Differentiation

"""
Ignores all adjoint definitions (i.e. `sensealg`) and proceeds to do standard
AD through the `solve` functions. Generally only used internally for implementing
discrete sensitivity algorithms.
"""
struct SensitivityADPassThrough <: SciMLBase.DEAlgorithm end

function ChainRulesCore.frule(::typeof(solve_up), prob,
  sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
  u0, p, args...;
  kwargs...)
  _solve_forward(prob, sensealg, u0, p, args...; kwargs...)
end

function ChainRulesCore.rrule(::typeof(solve_up), prob::SciMLBase.DEProblem,
  sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
  u0, p, args...;
  kwargs...)
  _solve_adjoint(prob, sensealg, u0, p, args...; kwargs...)
end

###
### Legacy Dispatches to be Non-Breaking
###

@deprecate concrete_solve(prob::SciMLBase.DEProblem, alg::Union{SciMLBase.DEAlgorithm,Nothing},
  u0=prob.u0, p=prob.p, args...; kwargs...) solve(prob, alg, args...; u0=u0, p=p, kwargs...)

function _solve_adjoint(prob, sensealg, u0, p, args...; merge_callbacks=true, kwargs...)

  _prob = if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  else
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_adjoint(_prob, nothing, sensealg, u0, p; kwargs...)
  else
    _concrete_solve_adjoint(_prob, args[1], sensealg, u0, p, Base.tail(args)...; kwargs...)
  end
end

function _solve_forward(prob, sensealg, u0, p, args...; merge_callbacks=true, kwargs...)

  _prob = if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  else
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_forward(prob, nothing, sensealg, u0, p; kwargs...)
  else
    _concrete_solve_forward(prob, args[1], sensealg, u0, p, Base.tail(args)...; kwargs...)
  end
end

function _concrete_solve_adjoint(args...; kwargs...)
  error("No adjoint rules exist. Check that you added `using DiffEqSensitivity`")
end

function _concrete_solve_forward(args...; kwargs...)
  error("No sensitivity rules exist. Check that you added `using DiffEqSensitivity`")
end


text:
abstract type Question end

mutable struct Stringq <: Question
    re::Regex
    label
    hint
    explanation
    placeholder
end

"""
    stringq(re::Regex; label="", hint="", explanation="", placeholder="")

Match string answer with regular expression

Arguments:

* `re`: a regular expression for grading

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

* `placeholder`: text shown when input widget is initially drawn

## Example

```
re = Regex("^abc")
stringq(re, label="First 3 letters...")
```

"""
stringq(re::Regex; label="", hint="", explanation="",  placeholder=nothing) =
    Stringq(re, label, hint, explanation, placeholder)

##
mutable struct Numericq <: Question
    val
    atol
    units
    label
    hint
    explanation
    placeholder
end

"""
    numericq(value, atol=1e-3; label="", hint="", units="", explanation="", placeholder=nothing)

Match a numeric answer

Arguments:

* `value`: the numeric answer

* `atol`: ``|answer - value| \\le atol`` is used to determine correctness

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `units`: a string indicating expected units

* `placeholder`: text shown when input widget is initially drawn

"""
function numericq(val, atol=1e-3, args...;
                  label="",
                  hint="", units::AbstractString="",
                  explanation="",
                  placeholder=nothing)

    Numericq(val, atol,  units, label, hint, explanation, placeholder)
end

numericq(val::Int; kwargs...) = numericq(val, 0; kwargs...)

##
mutable struct Radioq <: Question
    choices
    answer
    values
    labels
    label
    hint
    explanation
    inline
end

"""
    radioq(choices, answer; label="", hint="", explanation="", keep_order=false)

Multiple choice question (one of several)

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answer::Int`: index of correct choice

* `keep_order::Boolean`: if `true` keeps display order of choices, otherwise they are shuffled

* `inline::Bool`: hint to render inline (or not) if supported

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

## Example:

```
choices = ["beta", raw"``\\beta``", "`beta`"]
answer = 2
radioq(choices, answer; hint="Which is the Greek symbol?")
```

"""
function radioq(choices, answer::Integer;
                label="", hint="", explanation="",
                inline::Bool=(hint!=""),
                keep_order::Bool=false)
    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices
    !keep_order && shuffle!(inds)

    Radioq(choices[inds], findfirst(isequal(answer), inds),
           values, labels[inds], label, hint, explanation, inline)
end

##
mutable struct Buttonq <: Question
    choices
    answer
    label
    hint
    explanation
    colors
end

"""
    buttonq(choices, answer; label="", hint="", [colors])

Use buttons for multiple choice (one of many). Show answer after first click.

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answer::Int`: index of correct choice

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

## Example:

```
choices = ["beta", raw"``\\beta``", "`beta`"]
answer = 2
explanation = "The other two answers are not a symbol, but a name"
buttonq(choices, answer; label="Which is the Greek symbol?",
        explanation=explanation)
```
"""
function buttonq(choices, answer::Integer;
                 label="", hint="", explanation="",
                 colors=(GREEN="#FF0000AA",
                         RED = "#00AA33AA",
                         BLUE = nothing) #"#0033CC11",
                 )

    answers = [i == answer ? "correct" : "incorrect" for i ∈ eachindex(choices)]
    Buttonq(choices, answers, label, hint, explanation, colors)
end

##
mutable struct Multiq <: Question
    choices
    answer
    values
    labels
    label
    hint
    explanation
    inline
end

"""
    multiq(choices, answers; label="", hint="", explanation="", keep_order=false)

Multiple choice question (one *or more* of several)

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answers::Vector{Int}`: index of correct choice(s)

* `keep_order::Boolean`: if `true` keeps display order of choices, otherwise they are shuffled

* `inline::Bool`: hint to render inline (or not) if supported

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

## Example

```
choices = ["pear", "tomato", "banana"]
answers = [1,3]
multiplecq(choices, answers; label="yellow foods", hint="not the red one!")
```

"""
function multiq(choices, answers;
                label="", hint="", explanation="",
                inline::Bool=(hint!=""),
                keep_order::Bool=false)
    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices
    !keep_order && shuffle!(inds)

    Multiq(choices[inds], findall(in(answers), inds),
           values, labels[inds], label, hint, explanation, inline)
end

## ----
mutable struct MultiButtonq <: Question
    choices
    answer
    values
    labels
    label
    hint
    explanation
    inline
end

"""
    multibuttonq(choices, answers; label="", hint="", explanation="", keep_order=false)

Multiple choice question (zero, one *or more* of several) using buttons and a "done" button to reveal  the answers and whether the user input is correct.

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answers::Vector{Int}`: index of correct choice(s)

* `keep_order::Boolean`: if `true` keeps display order of choices, otherwise they are shuffled

* `inline::Bool`: hint to render inline (or not) if supported

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

## Example

```
choices = ["pear", "tomato", "banana"]
answers = [1,3]
multibuttonq(choices, answers; label="yellow foods", hint="not the red one!")
```

"""
function multibuttonq(choices, answers;
                label="", hint="", explanation="",
                inline::Bool=(hint!=""),
                      keep_order::Bool=false)

    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices
    !keep_order && shuffle!(inds)

    MultiButtonq(choices[inds], findall(in(answers), inds),
           values, labels[inds], label, hint, explanation, inline)
end

##
mutable struct Matchq <: Question
    questions
    choices
    answer
    label
    hint
    explanation
end

"""
    matchq(questions, choices, answers; label="", hint="", explanation="")
    matchq(d::Dictionary; label="", hint="", explanation="")

Use a drop down to select the right match for each question.

Arguments:

* `questions`: Indexable collection of questions

* `choices`: indexable collection of choices for each question. As seen in the example, choices can be formatted with markdown.

* `answers`: for each question, the index from `choices` of the correct answer

* `d`: As an alternative, a dictionary of questions and answers can be specified. The choices will be taken from the values then randomized, the answers will be computed

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

* `explanation`: text to display on a wrong selection

## Examples

```
questions = ("Select a Volvo", "Select a Mercedes", "Select an Audi")
choices = ("XC90", "A4", "GLE 350", "X1") # may be more than questions
answer = (1,3,2) # indices of correct
matchq(questions, choices, answer)
```

This example uses a dictionary to specify the questions and choices:

```
d = Dict("Select a Volvo" => "XC90", "Select a Mercedes" => "GLE 250")
matchq(d)
```

"""
function matchq(questions, choices, answers;
                label="", hint="", explanation="")

    @assert length(questions) == length(answers)

    Matchq(questions, choices, answers,
           label, hint, explanation)
end

function matchq(d::AbstractDict; kwargs...)

    questions = collect(keys(d))
    choices = collect(values(d))
    n = length(questions)
    inds = shuffle(1:n)

    matchq(questions, choices[inds], sortperm(inds); kwargs...)
end

"""
    booleanq(ans; [label, hint, explanation])
True of false questions:

## Example:

```
booleanq(true; label="Does it hurt...")
```

"""
function booleanq(ans::Bool;
                  labels::Vector=["true", "false"],
                  label="", hint::AbstractString="", explanation="",
                  inline::Bool=true)
    choices = labels[1:2]
    ans = 2 - ans
    radioq(choices, ans;
           label=label, hint=hint, explanation=explanation,
           inline=inline, keep_order=true)
end

"""
    yesnoq(ans; [label, hint, explanation])

Boolean question with `yes` or `no` labels.

Examples:

```
yesnoq("yes")
yesnoq(true)
```

"""
yesnoq(ans::AbstractString, args...; kwargs...) = radioq(["Yes", "No"], ans == "yes" ? 1 : 2, args...; keep_order=true, kwargs...)
yesnoq(ans::Bool, args...; kwargs...) = yesnoq(ans ? "yes" : "no", args...;kwargs...)

## --------

abstract type FillBlankQ <: Question end

mutable struct FillBlankChoiceQ <: FillBlankQ
    question
    choices
    answer
    label
    hint
    explanation
end

mutable struct FillBlankStringQ <: FillBlankQ
    question
    re::Regex
    label
    hint
    explanation
    placeholder
end

mutable struct FillBlankNumericQ <: FillBlankQ
    question
    val
    atol
    label
    hint
    explanation
    placeholder
end

"""
    fillblankq(question answer::Regex; placeholder=nothing, [label, hint, explanation])
    fillblankq(question, choices, answer; keep_order=false,[label, hint, explanation])
    fillblankq(question, val, atol=0; placeholder=nothing, [label, hint, explanation])

Present a fill-in-the-blank question where the blank can be a selection, a number, or a string graded by a regular expression.

* `question`: A string. Use `____` (4 or more under scores) to indicate the blank.

Other rguments from `stringq`, `radioq`, and `numericq`

## Examples
```
question = "The quick brown fox jumped over the ____ dog"
fillblankq(question, r"lazy")
fillblankq(question, ("lazy", "brown", "sleeping"), 1)
fillblankq("____ ``+ 2  = 4``", 2)
```
"""
fillblankq(question, answer::Regex; label="", hint="", explanation="",
           placeholder=nothing) =
    FillBlankStringQ(question, answer, label, hint, explanation, placeholder)

function fillblankq(question, choices, answer; label="", hint="", explanation="",
                    keep_order=false)

    if !keep_order
        inds = collect(1:length(choices))
        values = copy(inds)
        labels = choices
        !keep_order && shuffle!(inds)
        choices = choices[inds]
        answer = findfirst(isequal(answer), inds)
    end

    FillBlankChoiceQ(question, choices, answer, label, hint, explanation)
end

fillblankq(question, val::Real, atol=0; label="", hint="", explanation="",
           placeholder=nothing) =
    FillBlankNumericQ(question, val, atol, label, hint, explanation, placeholder)

## ------

mutable struct HotspotQ <: Question
    imgfile
    xs
    ys
    label
    hint
    explanation
    correct_answer
end

"""
    hotspotq(imagefile, xs, ys=(0,1); label="", hint="", explanation="",
        correct_answer=nothing)

Question type to check if user clicks in a specified region of an image.

* `imgfile`: File of an image. The images will be encoded and embedded in the web page.

* `xs`: iterable specifying `(xmin, xmax)` with `0 <= xmin <= xmax <= 1`

* `ys`: iterable specifying `(ymin, ymax)`

* `correct_answer`: A text snippet of JavaScript which can be specified to add more complicated logic to test if an answer is correct. It must use `x` and `y` for the coordinates of the click.

## Examples

```
using Plots
p1 = plot(x -> x^2, axis=nothing, legend=false)
p2 = plot(x -> x^3, axis=nothing, legend=false)
p3 = plot(x -> -x^2, axis=nothing, legend=false)
p4 = plot(x -> -x^3, axis=nothing, legend=false)
l = @layout [a b; c d]
p = plot(p1, p2, p3, p4, layout=l)
imgfile = tempname() * ".png"
savefig(p, imgfile)
hotspotq(imgfile, (0,1/2), (0, 1/2), label="What best matches the graph of ``f(x) = -x^4``?")
```

!!! note
    The display of the image is not adjusted by this question type and must
    be managed separately.

"""
function hotspotq(imgfile, xs, ys=(0,1);
                  label = "", hint="", explanation="",
                  correct_answer=nothing)
    HotspotQ(imgfile, xs, ys, label, hint, explanation, correct_answer)
end

## ----
mutable struct PlotlyLightQ <: Question
    p
    xs
    ys
    label
    hint
    explanation
    correct_answer
end

"""
    plotlylightq(p, xs=(-Inf, Inf), ys=(-Inf, Inf);
                 label="", hint="", explanation="",
                 correct_answer=nothing)

From a plotly graph of a function present a question about the `x-y` point on a graph shown on hover. (For figures with multiple graphs, the hover of the first one is taken. For parameterized plots, the hover may be computed in an unexpected manner.)

By default, correct answers select a value on the graph with `x` in the range specified by `xs` and `y` in the range specified by `ys`.

* `xs`: specifies interval for selected point `[x₀,x₁]`, defaults to `(-Inf,Inf)`

* `ys`: range `[y₀,y₁]`

* `correct_answer`: When specified, allows more advanced notions of correct. This is a JavaScript code snippet with `x` and `y` representing the hovered point on the graph that is highlighted on clicking.

## Examples

```
using PlotlyLight
xs = range(0, 2pi, length=100)
ys = sin.(xs)
p = Plot(Config(x=xs, y=ys))
plotlylightq(p, (3,Inf); label="Click a value with ``x>3``")
```

An example where the default grading script needs modification. Note
also, the `x`, `y` values refer to the *first* graph. (One could
modify their definition in `correct answer`; they are found through
`x=e.points[0].x`, `y=e.points[0].y`.)

```
xs = range(0, 2pi, length=100)
ys = sin.(xs)
p = Plot(Config(x=xs, y=ys));
push!(p.data, Config(x=xs, y=cos.(xs)));  # add layer
question = "Click a value where `sin(x)` is increasing"
# evalute pi/2 as no pi in JavaScript, also Inf -> Infinity
correct_answer = "((x >= 0 && x <= $(pi/2)) || (x >= $(3pi/2) && x <= $(2pi)))"
plotlylightq(p; label=question, correct_answer=correct_answer)
```

!!! note
    The use of `PlotlyLight` graphics works with `Weave` and `Pluto`, but is unusable from `Documenter`.

"""
function plotlylightq(p, xs=(-Inf, Inf), ys=(-Inf,Inf);
                  label = "", hint="", explanation="",
                  correct_answer=nothing)
    PlotlyLightQ(p, xs, ys, label, hint, explanation, correct_answer)
end


text:
module Todict

export todict
"todict(123, :asd) -> Dict(:asd => 123)"
todict(x, k::Symbol) = Dict( k=>x )

"123 |> todict(:asd) -> Dict(:asd => 123)"
todict(k::Symbol) = x->todict(x,k)

"todict(123, :a, :b=>x->x+1) -> Dict(:a => 123, :b => 124)"
function todict(x, k::Symbol, p::Pair)
 Dict( k=>x, p.first=>p.second(x) )
end

"123 |> todict(:a, :b=>x->x+1) -> Dict(:a => 123, :b => 124)"
todict(k::Symbol, p::Pair) = x->todict(x,k,p)

"""
todict(123, :a, :b=>x->x+1, :c=>x->x+2) -> Dict(:a=>123,:b=>124,:c=>125)
"""
function todict(x, k::Symbol, p::Pair, p2::Pair)
 Dict( k=>x, p.first=>p.second(x), p2.first=>p2.second(x) )
end

"123|>todict(:a, :b=>x->x+1, :c=>x->x+2) -> Dict(:a=>123,:b=>124,:c=>125)"
todict(k::Symbol, p::Pair, p2::Pair) =
    x->todict(x,k,p,p2)

"""
todict(123, :a, :b=>x->x+1, :c=>x->x+2, :d=>x->x+3) ->
    Dict(:a=>123,:b=>124,:c=>125,:d=>126)
"""
function todict(x, k::Symbol, p::Pair, p2::Pair, p3::Pair)
 Dict( k=>x, p[1]=>p[2](x), p2[1]=>p2[2](x), p3[1]=>p3[2](x) )
end

"""
123|>todict(:a, :b=>x->x+1, :c=>x->x+2, :d=>x->x+3) ->
    Dict(:a=>123,:b=>124,:c=>125,:d=>126)
"""
todict(k::Symbol, p::Pair, p2::Pair, p3::Pair) =
    x->todict(x,k,p,p2,p3)

"""
todict(123, :a, :b=>x->x+1, :c=>x->x+2, :d=>x->x+3, :e=>x->x+4) -> Dict(...)

todict(123, :a, :b=>x->x+1, :c=>x->x+2, :d=>x->x+3, :e=>x->x+4, :f=>x->x+5)
"""
function todict(
       x, k::Symbol, p::Pair, p2::Pair, p3::Pair, p4::Pair, pp::Pair...)
 dict = Dict( k=>x, p[1]=>p[2](x),
           p2[1]=>p2[2](x), p3[1]=>p3[2](x), p4[1]=>p4[2](x) )
 for (k,v) in pp
    get!(dict, k, v(x))
 end
 dict
end

"""
123|>todict(:a, :b=>x->x+1, :c=>x->x+2, :d=>x->x+3, :e=>x->x+4) -> Dict(...)
"""
todict(
    k::Symbol, p::Pair, p2::Pair, p3::Pair, p4::Pair, pp::Pair...) =
       x->todict(x, k, p, p2, p3, p4, pp...)

"todict( Dict(:a=>123), :a=>:b, x->x+1 ) -> Dict(:a=>123,:b=>124)"
function todict{P<:Pair{Symbol,Symbol}}(d::Dict, p::P, f::Function)
    get!(d, p.second, f(d[p.first]) )
    d
end

"Dict(:a=>123) |> todict( :a=>:b, x->x+1 ) -> Dict(:a=>123,:b=>124)"
todict{P<:Pair{Symbol,Symbol}}(p::P, f::Function) = d->todict(d,p,f)

"""
todict( Dict(:a=>123), :a=>:b, x->x+1, :b=>:c, x->x+1 ) ->
    Dict(:a=>123,:b=>124,:c=>125)
"""
todict{P<:Pair{Symbol,Symbol},F<:Function}(
    d::Dict, p1::P, f1::F, p2::P, f2::F, other...) =
       todict( todict(d,p1,f1), p2,f2, other... )

"Dict(:a=>123) |> todict( :a=>:b, x->x+1, :b=>:c, x->x+1 )"
todict{P<:Pair{Symbol,Symbol},F<:Function}(
    p1::P,f1::F, p2::P,f2::F, other...) =
   d->todict(d,p1,f1,p2,f2,other...)

end # module

max_stars_repo_name: closescreen/Todict

max_issues_repo_name: closescreen/Todict

max_forks_repo_name: closescreen/Todict


text:
printstyled("  diff_x!, int_x!\n", color=:light_green)

# Check that double-integration and double-differentiation are reversible
T = Float64
for i = 1:100
  fs = rand(T[10, 20, 40, 50, 100])
  δ = one(T)/fs
  N = 10^rand(1:5)
  N2 = div(N,2)
  x = randn(T, N)
  y = deepcopy(x)
  gaps = [1, rand(2:N2), rand(N2+1:N-1), N]

  diff_x!(x, gaps, fs)
  int_x!(x, gaps, δ)
  diff_x!(x, gaps, fs)
  int_x!(x, gaps, δ)
  if isapprox(x,y) == false
    println("test failed (part 1) after i = ", i, " trials (length(x) = ", length(x), ")")
  end
  @test isapprox(x,y)

  diff_x!(x, gaps, fs)
  diff_x!(x, gaps, fs)
  int_x!(x, gaps, δ)
  int_x!(x, gaps, δ)
  if isapprox(x,y) == false
    println("test failed (part 2) after i = ", i, " trials (length(x) = ", length(x), ")")
  end
  @test isapprox(x,y)

  int_x!(x, gaps, δ)
  int_x!(x, gaps, δ)
  diff_x!(x, gaps, fs)
  diff_x!(x, gaps, fs)
  if isapprox(x,y) == false
    println("test failed (part 3) after i = ", i, " trials (length(x) = ", length(x), ")")
  end
  @test isapprox(x,y)
end


text:
function measure(data, pattern)
    start = time()
    count = length(collect(eachmatch(pattern, data)))
    elapsed = time() - start
    elapsed_ms = 1000 * elapsed

    println(string(elapsed_ms) * " - " * string(count))
end

if length(ARGS) < 1
    println("Usage: julia benchmark.jl <filename>")
    exit(1)
end

data = open(f->read(f, String), ARGS[1])

measure(data, r"Twain");
measure(data, r"(?i)Twain");
measure(data, r"[a-z]shing");
measure(data, r"Huck[a-zA-Z]+|Saw[a-zA-Z]+");
measure(data, r"\b\w+nn\b");
measure(data, r"[a-q][^u-z]{13}x");
measure(data, r"Tom|Sawyer|Huckleberry|Finn");
measure(data, r"(?i)Tom|Sawyer|Huckleberry|Finn");
measure(data, r".{0,2}(Tom|Sawyer|Huckleberry|Finn)");
measure(data, r".{2,4}(Tom|Sawyer|Huckleberry|Finn)");
measure(data, r"Tom.{10,25}river|river.{10,25}Tom");
measure(data, r"[a-zA-Z]+ing");
measure(data, r"\s[a-zA-Z]{0,12}ing\s");
measure(data, r"([A-Za-z]awyer|[A-Za-z]inn)\s");
measure(data, r"[\"'][^\"']{0,30}[?!\\.][\"']");

max_stars_repo_name: olsaarik/regex-benchmark

max_issues_repo_name: olsaarik/regex-benchmark

max_forks_repo_name: olsaarik/regex-benchmark


text:
#xor-ml.jl - set up to test the xor function in julia.

using SigmoidNumbers
using GenML

import SigmoidNumbers.Sigmoid

#assign sigmoid to the pseudologistic
GenML.TF.sigmoid(x::Posit) = SigmoidNumbers.pseudologistic(x)
GenML.TF.dasysigmoid(x::Posit) = SigmoidNumbers.delta_psl(x)
GenML.CF.logloss(x::Posit) = SigmoidNumbers.pseudologcost(x)

macro widen(x)
  :(reinterpret(Posit{12}, $x))
end
macro noclamp(x)
  esc(:(reinterpret(F, $x)))
end
macro clamp(x)
  esc(:(SigmoidNumbers.__round(reinterpret(F, $x))))
end

function GenML.FCL.scaledouterproductfma{F<:Posit}(matrix::AbstractMatrix{F}, output_deltas::AbstractVector{F}, input_array::AbstractVector, alpha::F, i, o)
  #widen everything to posit16

  for idx = 1:o
    for jdx = 1:i
      if typeof(input_array) <: AbstractVector{Bool}
        scaledouterproductfma_prod = @widen(alpha) * @widen(output_deltas[idx]) * input_array[jdx]
        matrix[idx,jdx] = @noclamp(@widen(matrix[idx, jdx]) - scaledouterproductfma_prod) # + lambda * mlp.transitions[l_idx - 1])
      else
        scaledouterproductfma_prod = @widen(alpha) * (@widen(output_deltas[idx]) * @widen(input_array[jdx]))
        matrix[idx,jdx] = @noclamp(@widen(matrix[idx, jdx]) - scaledouterproductfma_prod) # + lambda * mlp.transitions[l_idx - 1])
      end
      @clamp(matrix[idx,jdx])
    end
  end
end

function xortrain(PType)

  println("working on unreliable xor data set with backpropagation")

  input_matrix = rand(Bool, 10, 500)
  training_results = Array{Bool,2}(1,500)
  training_results[:] = [input_matrix[1, col] $ input_matrix[2, col] $ (rand() < 0.05) for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{PType,(10,2,1)}(randn)

  for rounds = 1:25
  for idx = 1:500
    input_column = input_matrix[:, idx]
    training_column = training_results[:, idx]

    GenML.Optimizers.backpropagationoptimize(xornet, input_column, training_column, GenML.CF.crossentropy)
  end
  end

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > PType(0.5)) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for noisy xor: $wrongcount / 50 => $(2 * (50 - wrongcount)) % correct")

  return wrongcount
end

Base.randn{PType <: Sigmoid}(::Type{PType}, dims::Integer...) = map(PType, randn(dims...))

const trials = 1000
resdict=Dict{Int64, Float64}()
for nindex = 7:-1:6
  results = zeros(Int, trials)

  PType = Posit{nindex}

  for idx = 1:trials
    println("posit $(PType): round $idx")
    results[idx] = xortrain(PType)
  end

  println("results: ", sum(results) / trials)
  resdict[nindex] = sum(results) / trials

  #writedlm(string("./results/posit-$nindex.csv"), results, ',')

end

println(resdict)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
abstract type InfiniteOrderKernel <: ContinuousUnivariateDistribution end

function Base.show(io::IO, d::InfiniteOrderKernel)
    print(io, nameof(typeof(d)))
    print(io, " | bandwidth = ")
    print(io, d.h)
end

"""
    SincKernel(h) <: InfiniteOrderKernel

Implements the `SincKernel` with bandwidth `h` to be used for kernel density estimation
through the `KernelDensity.jl` package. The sinc kernel is defined as follows:
```math
K_{\\text{sinc}}(x) = \\frac{\\sin(x)}{\\pi x}
```
It is not typically used for kernel density estimation, because this kernel is not
a density itself. However, it is particularly well suited to deconvolution problems
and estimation of very smooth densities because its Fourier transform is the following:
```math
K^*_{\\text{sinc}}(t) = \\mathbf 1( t \\in [-1,1])
```
"""
struct SincKernel{H} <: InfiniteOrderKernel
    h::H
end
SincKernel() = SincKernel(DataBasedDefault())

function Empirikos._set_defaults(
    kernel::InfiniteOrderKernel,
    Zs::AbstractVector{<:Empirikos.AbstractNormalSample};
    kwargs...,
)
    Empirikos.skedasticity(Zs) == Empirikos.Homoskedastic() ||
        throw("Only implemented for Homoskedastic Gaussian data.")
    σ = std(Zs[1])
    n = length(Zs)
    h = σ / sqrt(log(n))
    @set kernel.h = h
end

Distributions.cf(a::SincKernel, t) = one(Float64) * (-1 / a.h <= t <= 1 / a.h)

function Distributions.pdf(a::SincKernel, t::Real)
    if t == zero(Float64)
        return (one(Float64) / pi / a.h)
    else
        return (sin(t / a.h) / pi / t)
    end
end

"""
    DeLaValleePoussinKernel(h) <: InfiniteOrderKernel

Implements the `DeLaValleePoussinKernel` with bandwidth `h` to be used for kernel density estimation
through the `KernelDensity.jl` package. The De La Vallée-Poussin kernel is defined as follows:
```math
K_V(x) = \\frac{\\cos(x)-\\cos(2x)}{\\pi x^2}
```
Its use case is similar to the [`SincKernel`](@ref), however it has the advantage of being integrable
(in the Lebesgue sense) and having bounded total variation. Its Fourier transform is the following:
```math
K^*_V(t) = \\begin{cases}
 1, & \\text{ if } |t|\\leq 1 \\\\
 0, &\\text{ if } |t| \\geq 2 \\\\
 2-|t|,& \\text{ if } |t| \\in [1,2]
 \\end{cases}
```
"""
struct DeLaValleePoussinKernel{H} <: InfiniteOrderKernel
    h::H
end

DeLaValleePoussinKernel() = DeLaValleePoussinKernel(DataBasedDefault())

function Distributions.cf(a::DeLaValleePoussinKernel, t)
    if abs(t * a.h) <= 1
        return (one(Float64))
    elseif abs(t * a.h) <= 2
        return (2 * one(Float64) - abs(t * a.h))
    else
        return (zero(Float64))
    end
end

function Distributions.pdf(a::DeLaValleePoussinKernel, t::Real)
    if t == zero(Float64)
        return (3 * one(Float64) / 2 / pi / a.h)
    else
        return (a.h * (cos(t / a.h) - cos(2 * t / a.h)) / pi / t^2)
    end
end

"""
    FlatTopKernel(h) < InfiniteOrderKernel

Implements the `FlatTopKernel` with bandwidth `h` to be used for kernel density estimation
through the `KernelDensity.jl` package. The flat-top kernel is defined as follows:
```math
K(x) = \\frac{\\sin^2(1.1x/2)-\\sin^2(x/2)}{\\pi x^2/ 20}.
```
Its use case is similar to the [`SincKernel`](@ref), however it has the advantage of being integrable
(in the Lebesgue sense) and having bounded total variation. Its Fourier transform is the following:
```math
K^*(t) = \\begin{cases}
 1, & \\text{ if } t|\\leq 1 \\\\
 0, &\\text{ if } |t| \\geq 1.1 \\\\
 11-10|t|,& \\text{ if } |t| \\in [1,1.1]
 \\end{cases}
```

```jldoctest
julia> Empirikos.FlatTopKernel(0.1)
FlatTopKernel | bandwidth = 0.1
```
"""
struct FlatTopKernel{H} <: InfiniteOrderKernel
    h::H
end

FlatTopKernel() = FlatTopKernel(DataBasedDefault())

function Distributions.cf(a::FlatTopKernel, t)
    if abs(t * a.h) <= 1
        return (one(Float64))
    elseif abs(t * a.h) <= 1.1
        return (one(Float64) - 10 * (abs(t * a.h) - one(Float64)))
    else
        return (zero(Float64))
    end
end

# TODO: This is *not* numerically stable.
# But presumably do not need it to be (since the KDE code just calls cf)
function Distributions.pdf(a::FlatTopKernel, t::Real)
    h = a.h
    th = t / h
    if t == zero(Float64)
        return ((0.55^2 - 0.5^2) * 20 / pi / h)
    else
        return ((abs2(sin(0.55 * th)) - abs2(sin(0.5 * th))) / pi / th^2 * 20 / h)
    end
end

"""
    InfinityNormDensityBand(;a_min,
                             a_max,
                             kernel  =  Empirikos.FlatTopKernel(),
                             bootstrap = :Multinomial,
                             nboot = 1000,
                             α = 0.05,
                             rng = Random.MersenneTwister(1)
                        )  <: FLocalization

This struct contains hyperparameters that will be used for constructing a neighborhood
of the marginal density. The steps of the method (and corresponding hyperparameter meanings)
are as follows
* First a kernel density estimate ``\\bar{f}`` with `kernel` is fit to the data.
* Second, a `bootstrap` (options: `:Multinomial` or `Poisson`) with `nboot` bootstrap replicates will be used to estimate ``c_n``, such that:
```math
\\liminf_{n \\to \\infty}\\mathbb{P}\\left[\\sup_{x \\in [a_{\\text{min}} , a_{\\text{max}}]} | \\bar{f}(x) - f(x)| \\leq c_ n\\right] \\geq 1-\\alpha
```
Note that the bound is valid from `a_min` to `a_max`. `α` is the nominal level and finally
`rng` sets the seed for the bootstrap samples.
"""
Base.@kwdef struct InfinityNormDensityBand <: FLocalization
    a_min = DataBasedDefault()
    a_max = DataBasedDefault()
    npoints::Integer = 1024
    kernel = FlatTopKernel()
    bootstrap = :Multinomial
    nboot::Integer = 1000
    α::Float64 = 0.05
    rng = Random.MersenneTwister(1)
end

Empirikos.vexity(::InfinityNormDensityBand) = Empirikos.LinearVexity()
function Empirikos._set_defaults(
    method::InfinityNormDensityBand,
    Zs::AbstractVector{<:Empirikos.AbstractNormalSample{<:Number}};
    hints...,
)
    Empirikos.skedasticity(Zs) == Empirikos.Homoskedastic() ||
        throw("Only implemented for Homoskedastic Gaussian data.")

    q = get(hints, :quantile, 0.1)
    a_min, a_max = quantile(response.(Zs), (q, 1 - q))
    if isa(method.a_min, DataBasedDefault)
        method = @set method.a_min = a_min
    end
    if isa(method.a_max, DataBasedDefault)
        method = @set method.a_max = a_max
    end
    method
end

"""
    FittedInfinityNormDensityBand

The result of running
```julia
StatsBase.fit(opt::InfinityNormDensityBand, Zs)```
Here `opt` is an instance
of [`InfinityNormDensityBand`](@ref) and `Zs` is a vector of [`AbstractNormalSample`](@ref)s
distributed according to a density ``f``..

## Fields:
* `a_min`,`a_max`, `kernel`: These are the same as the fields in `opt::InfinityNormDensityBand`.
* `C∞`: The half-width of the L∞ band.
* `fitted_kde`: The fitted `KernelDensity` object.
"""
Base.@kwdef struct FittedInfinityNormDensityBand{T<:Real,S,K} <: FittedFLocalization
    C∞::T
    a_min::T
    a_max::T
    fitted_kde::Any
    interp_kde::Any
    midpoints::S
    estimated_density::K
    boot_samples::Any
    method = nothing
end

Empirikos.vexity(method::FittedInfinityNormDensityBand) = Empirikos.LinearVexity()

function Empirikos.nominal_alpha(inftyband::FittedInfinityNormDensityBand)
    Empirikos.nominal_alpha(inftyband.method)
end

function StatsBase.fit(
    opt::InfinityNormDensityBand,
    Zs::AbstractVector{<:Empirikos.AbstractNormalSample{<:Number}};
    kwargs...,
)
    Empirikos.skedasticity(Zs) == Empirikos.Homoskedastic() ||
        throw("Only implemented for Homoskedastic Gaussian data.")
    opt = Empirikos.set_defaults(opt, Zs; kwargs...)

    @unpack a_min, a_max, npoints, kernel, nboot, α, bootstrap, rng = opt

    # deepcopying below to make sure RNG status for sampling here remains the same.
    rng = deepcopy(rng)

    res = certainty_banded_KDE(
        response.(Zs),
        a_min,
        a_max;
        npoints = npoints,
        rng = rng,
        kernel = kernel,
        bootstrap = bootstrap,
        nboot = nboot,
        α = α,
    )
    res = @set res.method = opt
    Z = Zs[1]
    normal_midpoints = [@set Z.Z = mdpt for mdpt in res.midpoints]
    res = @set res.midpoints = normal_midpoints
    res
end

function certainty_banded_KDE(
    Xs,
    a_min,
    a_max;
    kernel,
    rng = Random._GLOBAL_RNG,
    npoints = 4096,
    nboot = 1_000,
    bootstrap = :Multinomial,
    α = 0.5,
)

    h = kernel.h
    m = length(Xs)

    lo, hi = extrema(Xs)
    # avoid FFT wrap-around numerical difficulties
    lo_kde, hi_kde = min(a_min, lo - 6 * h), max(a_max, hi + 6 * h)
    midpts = range(lo_kde, hi_kde; length = npoints)

    fitted_kde = kde(Xs, KernelDensity.UniformWeights(m), midpts, kernel)
    interp_kde = InterpKDE(fitted_kde)

    # quantities only at points of interest
    midpts_idx = findall((midpts .>= a_min) .& (midpts .<= a_max))
    midpoints = fitted_kde.x[midpts_idx]
    estimated_density = fitted_kde.density[midpts_idx]

    C∞_boot = Vector{Float64}(undef, nboot)

    if bootstrap === :Multinomial
        multi = Multinomial(m, fill(1 / m, m))
    end
    for k = 1:nboot
        if bootstrap === :Poisson # Poisson bootstrap to estimate certainty band
            Z_pois = rand(rng, Poisson(1), m)
            ws = Weights(Z_pois / sum(Z_pois))
        elseif bootstrap === :Multinomial # Efron bootstrap
            ws = Weights(vec(rand(rng, multi, 1)))
        else
            throw(error("Only :Multinomial and :Poisson supported."))
        end
        f_kde_pois = kde(Xs, ws, midpts, kernel)
        C∞_boot[k] = maximum(abs.(estimated_density .- f_kde_pois.density[midpts_idx]))
    end

    C∞ = quantile(C∞_boot, 1 - α)

    FittedInfinityNormDensityBand(
        C∞ = C∞,
        a_min = a_min,
        a_max = a_max,
        fitted_kde = fitted_kde,
        interp_kde = interp_kde,
        midpoints = midpoints,
        estimated_density = estimated_density,
        boot_samples = C∞_boot,
    )
end

@recipe function f(ctband::FittedInfinityNormDensityBand; subsample = 300)
    y_all = ctband.estimated_density
    x_all = response.(ctband.midpoints)

    n = length(y_all)
    _step = div(n - 2, subsample)
    idxs = [1; 2:_step:(n-1); n]
    x = x_all[idxs]
    y = y_all[idxs]

    xlims --> (a_min, a_max)
    ylims --> (0, maximum(y) * 1.1)
    yguide --> "Density"

    α = nominal_alpha(ctband)
    ci_level = 100 * (1 - α)

    cis_ribbon = ctband.C∞
    fillalpha --> 0.36
    seriescolor --> "#018AC4"
    ribbon --> cis_ribbon
    framestyle --> :box
    legend --> :topleft
    label --> "$ci_level% CI"

    x, y
end

function Empirikos.flocalization_constraint!(
    model,
    ctband::FittedInfinityNormDensityBand,
    prior::Empirikos.PriorVariable,
)

    C∞ = ctband.C∞
    _midpoints = ctband.midpoints
    _density_values = ctband.estimated_density

    # TODO: perhaps move this evaluation below
    _min, _max = extrema(Empirikos.MarginalDensity(_midpoints[1]))
    for (index, Z) in enumerate(_midpoints)
        _density_hat = _density_values[index]
        marginal_pdf = pdf(prior, Z::EBayesSample)

        if _density_hat + C∞ < _max
            @constraint(model, marginal_pdf <= _density_hat + C∞)
        end
        if _density_hat - C∞ > _min
            @constraint(model, marginal_pdf >= _density_hat - C∞)
        end
    end
    model
end


text:
"""
**Nadam constructor (Nesterov Adam)**

```julia
     Nadam(;η::Real=0.001, β₁::Real=0.9, β₂::Real=0.999, ϵ::Real=10e-8)
```

Algorithm :

```math
\\begin{align*}
m_t =& \\beta_1 m_{t-1} + (1-\\beta_1)g_t\\\\
v_t =& \\beta_2 v_{t-1} + (1-\\beta_2)g_t^2\\\\
\\hat{m}_t =& \\frac{m_t}{1-\\beta_1^t}\\\\
\\hat{v}_t =& \\frac{v_t}{1-\\beta_2^t}\\\\
\\Delta x_t =&\\frac{\\eta}{\\sqrt{\\hat{v}_t+\\epsilon}}(\\beta_1\\hat{m}_t+\\frac{(1-\\beta_1)g_t}{1-\\beta_1^t})
\\end{align*}
```

[Algorithm Reference](http://cs229.stanford.edu/proj2015/054_report.pdf) and [other reference](http://www.cs.toronto.edu/~fritz/absps/momentum.pdf)
"""
mutable struct Nadam <: Optimizer
    opt_type::String
    t::Int64
    ϵ::Float64
    η::Float64
    β₁::Float64
    β₂::Float64
    m_t::AbstractArray
    v_t::AbstractArray
end

"Construct Nadam optimizer"
function Nadam(;η::Real=0.001, β₁::Real=0.9, β₂::Real=0.999, ϵ::Real=10e-8)
    @assert η > 0.0 "η must be greater than 0"
    @assert β₁ > 0.0 "β₁ must be greater than 0"
    @assert β₂ > 0.0 "β₂ must be greater than 0"
    @assert ϵ > 0.0 "ϵ must be greater than 0"

    Nadam("Nadam", 0, ϵ, η, β₁, β₂, [], [])
end

params(opt::Nadam) = "ϵ=$(opt.ϵ), η=$(opt.η), β₁=$(opt.β₁), β₂=$(opt.β₂)"

function update(opt::Nadam, g_t::AbstractArray{T}) where {T<:Real}
    # resize biased moment estimates if first iteration
    if opt.t == 0
        opt.m_t = zero(g_t)
        opt.v_t = zero(g_t)
    end

    # update timestep
    opt.t += 1

    # update biased first moment estimate
    opt.m_t = opt.β₁ * opt.m_t + (one(T) - opt.β₁) * g_t

    # update biased second raw moment estimate
    opt.v_t = opt.β₂ * opt.v_t + (one(T) - opt.β₂) * ((g_t) .^2)

    # compute bias corrected first moment estimate
    m̂_t = opt.m_t / (one(T) - opt.β₁^opt.t)

    # compute bias corrected second raw moment estimate
    v̂_t = opt.v_t / (one(T) - opt.β₂^opt.t)

    # apply update
    ρ = opt.η ./ (sqrt.(v̂_t .+ opt.ϵ))
    ρ .*= (opt.β₁ * m̂_t + (one(T) - opt.β₁) * g_t / (one(T) - opt.β₁^opt.t))

    return ρ
end