text:
"Hamiltonian Monte Carlo Sampling"
struct HMCAlg <: SamplingAlgorithm end

"Hamiltonian Monte Carlo Sampling"
const HMC = HMCAlg()
isapproximate(::HMCAlg) = true
defΩ(::HMCAlg) = SimpleΩ{ID, Float64}

"Hamiltonian monte carlo with leapfrog integration: https://arxiv.org/pdf/1206.1901.pdf"
function hmc(U, ∇U, nsteps, stepsize, current_q::Vector)
  q = transform.(current_q)
  p = randn(length(q))
  current_p = p

  # Make a half step for momentum at beginning
  # Rejects proposals outside domain TODO: Something smarter
  # any(notunit, q) && return (current_q, false)
  invq = inv_transform.(q)
  p = p - stepsize * ∇U(invq) .* jacobian.(q) / 2.0

  for i = 1:nsteps
    # Half step for the position and momentum
    q = q .+ stepsize .* p
    if i != nsteps
      invq = inv_transform.(q)
      p = p - stepsize * ∇U(invq) .* jacobian.(q) ./ 2.0
    end
  end

  # Make half a step for momentum at the end
  # any(notunit, q) && return current_q, false
  invq = inv_transform.(q)
  p = p .- stepsize .* ∇U(invq) .* jacobian.(q) ./ 2.0

  # Evaluate the potential and kinetic energies at start and end
  current_U = U(current_q)
  current_K =  sum(current_p.^2) / 2.0
  proposed_U = U(invq)
  proposed_K = sum(p.^2) / 2.0

  if log(rand()) < current_U - proposed_U + current_K - proposed_K
    return (proposed_U, invq, true) # accept ω
  else
    return (current_U, current_q, false)  # reject ω
  end
end

"Sample from `x | y == true` with Hamiltonian Monte Carlo"
function Base.rand(y::RandVar,
                   n::Integer,
                   alg::HMCAlg,
                   ΩT::Type{OT};
                   nsteps = 10,
                   stepsize = 0.001) where {OT <: Ω}
  ω = ΩT()
  indomainₛ(y, ω) # Initialize omega
  ωvec = linearize(ω)

  ωsamples = ΩT[]
  U(ω) = -logerr(indomainₛ(y, ω))
  U(ωvec::Vector) = U(unlinearize(ωvec, ω))
  ∇U(ωvec) = gradient(y, ω, ωvec)

  accepted = 0
  for i = 1:n
    p_, ωvec, wasaccepted = hmc(U, ∇U, nsteps, stepsize, ωvec)
    ω_ = unlinearize(ωvec, ω)
    push!(ωsamples, unlinearize(ωvec, ω))
    if wasaccepted
      accepted += 1
    end
    lens(Loop, (ω = ω_, accepted = accepted, p = p_, i = i))
  end
  [applynotrackerr(y, ω_) for ω_ in ωsamples]
end


text:
module CutCellDG

using LinearAlgebra
using SparseArrays
using NearestNeighbors
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
import ImplicitDomainQuadrature: extend

include("utilities.jl")
include("cell_map.jl")
include("cg_mesh.jl")
include("dg_mesh.jl")
include("periodic_bc.jl")
include("levelset.jl")
include("cut_mesh.jl")

include("cell_quadratures.jl")
include("face_quadratures.jl")
include("interface_quadratures.jl")
include("cell_merging.jl")

include("hooke_stiffness.jl")

include("weak_form.jl")
include("interface_conditions.jl")
include("transformation_strain.jl")

include("seed_zero_levelset.jl")
include("levelset_closest_points.jl")
# include("levelset_reinitialize.jl")
include("levelset_propagate.jl")

include("assembly.jl")
include("assemble_hermite_levelset_initializers.jl")
include("assemble_displacement_bilinear_forms.jl")
include("assemble_interelement_condition.jl")
include("assemble_interface_condition.jl")
include("assemble_penalty_displacement_bc.jl")
include("assemble_penalty_displacement_component_bc.jl")
include("assemble_body_force.jl")
include("assemble_traction_force_linear_form.jl")
include("assemble_traction_force_component_linear_form.jl")
include("assemble_transformation_strain.jl")
include("postprocess.jl")

end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Roots of a quadratic function

function quad_roots{S<:Real,T<:Real,U<:Real}(x::S, y::T, z::U)
    (a, b, c) = promote(float(x), y, z)
    2eps(a) < abs(a) || return [-c/b]
    disc = b^2 - 4a*c
    2eps(b^2) < abs(disc) || return [-sqrt(c/a)]
    if disc < 0
        disc += 0.0im
    end
    d = sqrt(disc)
    if b < 0
        d += -b
        return [d/2a, 2c/d]
    else
        d = -b - d
        return [2c/d, d/2a]
    end
end

a = {1, 1, 1.0, 10}
b = {10, 2, -10.0^9, 1}
c = {1, 1, 1, 1}

for i in 1:length(a)
    pstr = @sprintf "%fx^2 + %fx + %f" a[i] b[i] c[i]
    println("The roots of ", pstr, " are: ")
    println("    ", join(quad_roots(a[i], b[i], c[i]), " and "), "\n")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module TestMLJUnit

using Test

import MLJTuning

using TreeParzen

struct DummyModel
    init::Bool # dummy field to initialise because we can't do DummyModel()
end
DummyModel(;kwargs...) = DummyModel(true) # we literally don't care, we just need a constructor

function setup(;n_simultaneous=1, n_startup=20, space=Dict{Symbol}, suggest=Dict{Symbol}[])

    tuning = MLJTreeParzen.MLJTreeParzenTuning(;max_simultaneous_draws=n_simultaneous, random_trials=n_startup)
    state = MLJTreeParzen.MLJTreeParzenSpace(space, suggest)

    return tuning, state

end

function complete_trials(array_of_stoof)

    trials = last.(array_of_stoof)

    for trial in trials
        TreeParzen.tell!(trial, rand(Float64))
    end

    return [(DummyModel(), (trial_object=trial,)) for trial in trials]

end

@testset "MLJ: default_n" begin

    obj = MLJTreeParzen.MLJTreeParzenTuning()
    @test MLJTuning.default_n(obj, nothing) == 50

end

@testset "MLJ: setup" begin

    tuning = MLJTreeParzen.MLJTreeParzenTuning()
    model = nothing # supposed to be a model but nobody is checking and we dont need it so lets prove it
    verbosity = nothing

    space_dict = Dict(:x => HP.Uniform(:x, 0.0, 1.0))
    suggestions = Dict(:x => 0.222)
    suggestions_list = [Dict(:x => 0.222), Dict(:x => 0.444)]
    preconstructed_space = MLJTreeParzen.MLJTreeParzenSpace(space_dict)
    preconstruted_suggest = MLJTreeParzen.MLJTreeParzenSpace(space_dict, suggestions)
    preconstruted_suggestions = MLJTreeParzen.MLJTreeParzenSpace(space_dict, suggestions_list)

    # test with just a dict space
    @test isa(MLJTuning.setup(tuning, model, space_dict, verbosity), MLJTreeParzen.MLJTreeParzenSpace)

    # test with preconstructed space and that we get same thing
    @test isa(MLJTuning.setup(tuning, model, preconstructed_space, verbosity), MLJTreeParzen.MLJTreeParzenSpace)
    @test MLJTuning.setup(tuning, model, preconstructed_space, verbosity) == preconstructed_space

    # test with preconstructed space with a sugggestion and we get the same thing
    @test isa(MLJTuning.setup(tuning, model, preconstruted_suggest, verbosity), MLJTreeParzen.MLJTreeParzenSpace)
    @test MLJTuning.setup(tuning, model, preconstruted_suggest, verbosity) == preconstruted_suggest

    # test with preconstructed space with multiple suggestions and we get the same thing
    @test isa(MLJTuning.setup(tuning, model, preconstruted_suggestions, verbosity), MLJTreeParzen.MLJTreeParzenSpace)
    @test MLJTuning.setup(tuning, model, preconstruted_suggestions, verbosity) == preconstruted_suggestions

end

@testset "MLJ: models!" begin

    # scenarios to cover (or properties to check):
    # 1) with suggestions < num startup (return num startup followed by num simultaenous draw)
    # 2) with suggestions > num startup (throw exception)
    # 3) with suggestions == num startup (return num startup followed by num simultaenous draw)
    # 4) no suggestions (return num startup followed by num simultaenous draw)
    # 5) round 1 always returns num startup and all subsequent rounds returns num simultaneous draw

    space = Dict(:x => HP.Uniform(:x, -5., 5.))
    testmodel = DummyModel()
    fakehist = complete_trials([(nothing, TreeParzen.ask(space)) for i in 1:100]) # do too many, we can cut it down

    @testset "no suggestions" begin

        suggestions = Dict{Symbol}[]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output = MLJTuning.models!(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3

    end

    @testset "less suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output = MLJTuning.models!(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3
        @test getindex.(getproperty.(last.(output), :hyperparams), :x)[1:2] == [1, 2]

    end

        @testset "exact suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2), Dict(:x => 3)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output = MLJTuning.models!(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3
        @test getindex.(getproperty.(last.(output), :hyperparams), :x) == [1, 2, 3]

    end

        @testset "too many suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2), Dict(:x => 3), Dict(:x => 4)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        @test_throws ArgumentError MLJTuning.models!(tuning, testmodel, nothing, state, 0, 0)

    end

    @testset "num draws" begin

        for sim_value in 1:10

            suggestions = Dict{Symbol}[]
            tuning, state = setup(;n_simultaneous=sim_value,space=space, suggest=suggestions)

            output = MLJTuning.models!(tuning, testmodel, nothing, state, 0, 0)

            @test output isa Vector{Tuple{Any, TreeParzen.Trials.Trial}}
            @test length(output) == 20

            for histlen in 21:30

                output = MLJTuning.models!(tuning, testmodel, fakehist[1:histlen], state, 0, 0)
                @test length(output) == sim_value

            end

        end

    end

end

@testset "MLJ: result" begin

    # this will change but is easy, just different signature
    sample_eval = (measure="auc", measurement=1.0, donkeyfield=22, blabla="shorts")
    trial = TreeParzen.ask(Dict(:x => 0.222))
    complete_trial = tell!(deepcopy(trial), 0.52)
    tuning = MLJTreeParzen.MLJTreeParzenTuning()
    history = (measure="auc", measurement=5.2, trial_object=complete_trial)

    expected_trial_obj = deepcopy(trial)
    tell!(expected_trial_obj, sample_eval.measurement)
    expected_result = (measure=sample_eval.measure, measurement=sample_eval.measurement, trial_object=expected_trial_obj)

    result = MLJTuning.result(tuning, history, nothing, sample_eval, trial)

    @test result.measurement == expected_result.measurement
    @test result.measure == expected_result.measure
    @test result.trial_object.loss == expected_result.trial_object.loss

    # again but this time with empty history too, shouldnt make a difference though
    expected_trial_obj = deepcopy(trial)
    tell!(expected_trial_obj, sample_eval.measurement)
    expected_result = (measure=sample_eval.measure, measurement=sample_eval.measurement, trial_object=expected_trial_obj)

    result = MLJTuning.result(tuning, nothing, nothing, sample_eval, trial)

    @test result.measurement == expected_result.measurement
    @test result.measure == expected_result.measure
    @test result.trial_object.loss == expected_result.trial_object.loss

end

end # module
true


text:
module IoTest

using Avro.Io
using Avro.Schemas
using Test

const BOOLEAN_EXAMPLES =
    [
        (true, Schemas.BOOLEAN, [0x01]),
        (false, Schemas.BOOLEAN, [0x00]),
    ]

const INT_EXAMPLES =
    [
        (Int32(0), Schemas.INT, [0x00]),
        (Int32(-1), Schemas.INT, [0x01]),
        (Int32(1), Schemas.INT, [0x02]),
        (Int32(-2), Schemas.INT, [0x03]),
        (Int32(2), Schemas.INT, [0x04]),
        (Int32(-64), Schemas.INT, [0x7f]),
        (Int32(64), Schemas.INT, [0x80, 0x01]),
        (Int32(8192), Schemas.INT, [0x80, 0x80, 0x01]),
        (Int32(-8193), Schemas.INT, [0x81, 0x80, 0x01]),
    ]

const LONG_EXAMPLES =
    [
        (Int64(0), Schemas.LONG, [0x00]),
        (Int64(-1), Schemas.LONG, [0x01]),
        (Int64(1), Schemas.LONG, [0x02]),
        (Int64(-2), Schemas.LONG, [0x03]),
        (Int64(2), Schemas.LONG, [0x04]),
        (Int64(-64), Schemas.LONG, [0x7f]),
        (Int64(64), Schemas.LONG, [0x80, 0x01]),
        (Int64(8192), Schemas.LONG, [0x80, 0x80, 0x01]),
        (Int64(-8193), Schemas.LONG, [0x81, 0x80, 0x01]),
    ]

const FLOAT_EXAMPLES =
    [
        (Float32(0), Schemas.FLOAT, zeros(UInt8, 4)),
        (Float32(0.1), Schemas.FLOAT, [0xCD, 0xCC, 0xCC, 0x3D]),
        (Float32(0.2), Schemas.FLOAT, [0xCD, 0xCC, 0x4C, 0x3E]),
        (Float32(1.0), Schemas.FLOAT, [0x00, 0x00, 0x80, 0x3f]),
        (Inf32, Schemas.FLOAT, [0x00, 0x00, 0x80, 0x7F]),
        (-Inf32, Schemas.FLOAT, [0x00, 0x00, 0x80, 0xFF]),
        (NaN32, Schemas.FLOAT, [0x00, 0x00, 0xC0, 0x7F]),
    ]

const DOUBLE_EXAMPLES =
    [
        (Float64(0), Schemas.DOUBLE, zeros(UInt8, 8)),
        (Float64(0.1), Schemas.DOUBLE, [0x9A, 0x99, 0x99, 0x99, 0x99, 0x99, 0xb9, 0x3F]),
        (Float64(0.2), Schemas.DOUBLE, [0x9A, 0x99, 0x99, 0x99, 0x99, 0x99, 0xc9, 0x3F]),
        (Float64(1.0), Schemas.DOUBLE, [0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xF0, 0x3F]),
        (Inf64, Schemas.DOUBLE, [0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xF0, 0x7F]),
        (-Inf64, Schemas.DOUBLE, [0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xF0, 0xFF]),
        (NaN64, Schemas.DOUBLE, [0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xF8, 0x7F]),
    ]

const STRING_EXAMPLES =
    [
        ("foo", Schemas.STRING, [0x06, 0x66, 0x6f, 0x6f])
        ("", Schemas.STRING, [0x00])
    ]

buffer = IOBuffer()
encoder = BinaryEncoder(buffer)
decoder = BinaryDecoder(buffer)

@testset "Encoding" begin
    @testset "Integer" for (input, schema, expected) in INT_EXAMPLES
        # Encode the datum
        bytes_written = encode_int(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_int(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written

        take!(buffer)
    end

    @testset "Long" for (input, schema, expected) in LONG_EXAMPLES
        # Encode the datum
        bytes_written = encode_long(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_long(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Float" for (input, schema, expected) in FLOAT_EXAMPLES
        # Encode the datum
        bytes_written = encode_float(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_float(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        if !isnan(input)
            @test input == output
        end
        @test length(expected) == bytes_written
    end

    @testset "Double" for (input, schema, expected) in DOUBLE_EXAMPLES
        # Encode the datum
        bytes_written = encode_double(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_double(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        if !isnan(input)
            @test input == output
        end
        @test length(expected) == bytes_written
    end

    @testset "Boolean" for (input, schema, expected) in BOOLEAN_EXAMPLES
        # Encode the datum
        bytes_written = encode_boolean(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_boolean(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "String" for (input, schema, expected) in STRING_EXAMPLES
        # Encode the datum
        bytes_written = encode_string(encoder, input)

        # Decode the datum
        seekstart(buffer)
        output = decode_string(decoder)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test output == input
        @test length(expected) == bytes_written
    end
end

end


text:
using RDatasets, GLM, Statistics

df = dataset("MASS", "cpus")
df.Freq = map( x->10^9/x , df.CycT)

model = lm(@formula(Perf ~ MMax + Cach + ChMax + Freq), df)
pred(x) = round(coef(model)'*vcat(1,x),digits = 3)

println("n = ", size(df)[1])
println("(Avg,Std) of observed performance: ", (mean(df.Perf),std(df.Perf)))
println(model)
println("Estimated performance for computer A: ", pred([32000, 32, 32, 4*10^7]))
println("Estimated performance for computer B: ", pred([32000, 16, 32, 6*10^7]))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# # Foundations of Bilevel Programming: Example 5
# This example is from the book _Foundations of Bilevel Programming_ by Stephan
# Dempe, Chapter 8.1, Page 255. [url](https://www.springer.com/gp/book/9781402006319)

# Here, only the second level is described

# Model of the problem
# First level
# ```math
# \min 0,\\
# ```
# Second level
# ```math
# \min x,\\
# \notag s.t.\\
# x+y \leq 2,\\
# x-y \leq 2,\\
# -4x+5y \leq 10,\\
# -4x-5y \leq 10,\\
# ```

using BilevelJuMP
using Ipopt
using JuMP
using Test

model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-9))

# Global variables
atol = 1e-3

# First we need to create all of the variables in the upper and lower problems:

# Upper level variables
@variable(Upper(model), y, start = 0)

#Lower level variables
@variable(Lower(model), x, start = 0)

# Then we can add the objective and constraints of the upper problem:

# Upper level objecive function
@objective(Upper(model), Min, 0 * y + 0)

# Followed by the objective and constraints of the lower problem:

# Lower objective function
@objective(Lower(model), Min, x)

# Lower constraints
@constraint(Lower(model), x + y <= 2)
@constraint(Lower(model), x - y <= 2)
@constraint(Lower(model), -4x + 5y <= 10)
@constraint(Lower(model), -4x - 5y <= 10)

# Initial Starting conditions  #src

# Now we can solve the problem and verify the solution again that reported by
# Dempe.

optimize!(model)
primal_status(model)
termination_status(model)


text:
# OpenCL.Event

abstract type CLEvent <: CLObject end

mutable struct Event <: CLEvent
    id :: CL_event

    function Event(evt_id::CL_event; retain=false)
        if retain
            @check api.clRetainEvent(evt_id)
        end
        evt = new(evt_id)
        finalizer(_finalize, evt)
        return evt
    end
end

# wait for completion before running finalizer
mutable struct NannyEvent <: CLEvent
    id::CL_event
    obj::Any

    function NannyEvent(evt_id::CL_event, obj::Any; retain=false)
        if retain
            @check api.clRetainEvent(evt_id)
        end
        nanny_evt = new(evt_id, obj)
        finalizer(nanny_evt) do x
            x.id != C_NULL && wait(x)
            x.obj = nothing
            _finalize(x)
        end
        nanny_evt
    end
end

function _finalize(evt::CLEvent)
    if evt.id != C_NULL
        @check api.clReleaseEvent(evt.id)
        evt.id = C_NULL
    end
end

NannyEvent(evt::Event, obj::Any; retain=false) = NannyEvent(evt.id, obj, retain=retain)

Base.pointer(evt::CLEvent) = evt.id

function Base.show(io::IO, evt::Event)
    ptr_val = convert(UInt, Base.pointer(evt))
    ptr_address = "0x$(string(ptr_val, base = 16, pad = Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Event(@$ptr_address)")
end

Base.getindex(evt::CLEvent, evt_info::Symbol) = info(evt, evt_info)

@ocl_v1_1_only begin

    mutable struct UserEvent <: CLEvent
        id :: CL_event

        function UserEvent(evt_id::CL_event, retain=false)
            if retain
                @check api.clRetainEvent(evt_id)
            end
            evt = new(evt_id)
            finalizer(_finalize, evt)
            return evt
        end
    end

    function UserEvent(ctx::Context; retain=false)
        status = Ref{CL_int}()
        evt_id = api.clCreateUserEvent(ctx.id, status)
        if status[] != CL_SUCCESS
            throw(CLError(status[]))
        end
        try
            return UserEvent(evt_id, retain)
        catch err
            @check api.clReleaseEvent(evt_id)
            throw(err)
        end
    end

    function Base.show(io::IO, evt::UserEvent)
        ptr_val = convert(UInt, Base.pointer(evt))
        ptr_address = "0x$(string(ptr_val, base = 16, pad = Sys.WORD_SIZE>>2))"
        print(io, "OpenCL.UserEvent(@$ptr_address)")
    end

    function complete(evt::UserEvent)
        @check api.clSetUserEventStatus(evt.id, CL_COMPLETE)
        return evt
    end
end

struct _EventCB
    handle :: Ptr{Nothing}
    evt_id :: CL_event
    status :: CL_int
end

function event_notify(evt_id::CL_event, status::CL_int, payload::Ptr{Nothing})
    ptr = convert(Ptr{_EventCB}, payload)
    handle = unsafe_load(ptr).handle

    val = _EventCB(handle, evt_id, status)
    unsafe_store!(ptr, val)

    # Use uv_async_send to notify the main thread
    ccall(:uv_async_send, Nothing, (Ptr{Nothing},), handle)
    nothing
end

function add_callback(evt::CLEvent, callback::Function)
    event_notify_ptr = @cfunction(event_notify, Nothing,
                                  (CL_event, CL_int, Ptr{Cvoid}))

    # The uv_callback is going to notify a task that,
    # then executes the real callback.
    cb = Base.AsyncCondition()
    GC.@preserve cb begin

        # Storing the results of our c_callback needs to be
        # isbits && isimmutable
        r_ecb = Ref(_EventCB(Base.unsafe_convert(Ptr{Cvoid}, cb), 0, 0))

        @check api.clSetEventCallback(evt.id, CL_COMPLETE, event_notify_ptr, r_ecb)

        @async begin
           try
             Base.wait(cb)
             ecb = r_ecb[]
             callback(ecb.evt_id, ecb.status)
           catch
             rethrow()
           finally
             Base.close(cb)
           end
        end
    end
end

function wait(evt::CLEvent)
    evt_id = [evt.id]
    @check api.clWaitForEvents(cl_uint(1), pointer(evt_id))
    return evt
end

function wait(evts::Vector{CLEvent})
    evt_ids = [evt.id for evt in evts]
    if !isempty(evt_ids)
        nevents = cl_uint(length(evt_ids))
        @check api.clWaitForEvents(nevents, pointer(evt_ids))
    end
    return evts
end

@ocl_v1_2_only begin
    function enqueue_marker_with_wait_list(q::CmdQueue,
                                           wait_for::Vector{CLEvent})
        n_wait_events = cl_uint(length(wait_for))
        wait_evt_ids = [evt.id for evt in wait_for]
        ret_evt = Ref{CL_event}()
        @check api.clEnqueueMarkerWithWaitList(q.id, n_wait_events,
                                               isempty(wait_evt_ids) ? C_NULL : wait_evt_ids,
                                               ret_evt)
        @return_event ret_evt[]
    end

    function enqueue_barrier_with_wait_list(q::CmdQueue,
                                            wait_for::Vector{CLEvent})
        n_wait_events = cl_uint(length(wait_for))
        wait_evt_ids = [evt.id for evt in wait_for]
        ret_evt = Ref{CL_event}()
        @check api.clEnqueueBarrierWithWaitList(q.id, n_wait_events,
                                                isempty(wait_evt_ids) ? C_NULL : wait_evt_ids,
                                                ret_evt)
        @return_event ret_evt[]
    end
end

function enqueue_marker(q::CmdQueue)
    evt = Ref{CL_event}()
    @check api.clEnqueueMarker(q.id, evt)
    @return_event evt[]
end
@deprecate enqueue_marker enqueue_marker_with_wait_list

function enqueue_wait_for_events(q::CmdQueue, wait_for::Vector{T}) where {T<:CLEvent}
    n_wait_events = cl_uint(length(wait_for))
    wait_evt_ids = [evt.id for evt in wait_for]
    @check api.clEnqueueWaitForEvents(q.id, n_wait_events,
                                      isempty(wait_evt_ids) ? C_NULL : pointer(wait_evt_ids))
end

function enqueue_wait_for_events(q::CmdQueue, wait_for::CLEvent)
    enqueue_wait_for_events(q, [wait_for])
end

function enqueue_barrier(q::CmdQueue)
    @check api.clEnqueueBarrier(q.id)
    return q
end
@deprecate enqueue_barrier enqueue_barrier_with_wait_list

cl_event_status(s::Symbol) = begin
    if s == :queued
        return CL_QUEUED
    elseif s == :submitted
        return CL_SUBMITTED
    elseif s == :running
        return CL_RUNNING
    elseif s == :complete
        return CL_COMPLETE
    else
        throw(ArgumentError("unrecognized status symbol :$s"))
    end
end

macro profile_info(func, profile_info)
    quote
        function $(esc(func))(evt::CLEvent)
            time = Ref{CL_long}(0)
            err_code = api.clGetEventProfilingInfo(evt.id, $(esc(profile_info)),
                                                   sizeof(CL_ulong), time, C_NULL)
            if err_code != CL_SUCCESS
                if err_code == CL_PROFILING_INFO_NOT_AVAILABLE
                    if evt[:status] != :complete
                        #TODO: evt must have status complete before it can be profiled
                        throw(CLError(err_code))
                    else
                        #TODO: queue must be created with :profile argument
                        throw(CLError(err_code))
                    end
                end
                throw(CLError(err_code))
            end
            return time[]
        end
    end
end

function info(evt::CLEvent, evt_info::Symbol)
    command_queue(evt::CLEvent) = begin
        cmd_q = Ref{CL_command_queue}()
        @check api.clGetEventInfo(evt.id, CL_EVENT_COMMAND_QUEUE,
                                  sizeof(CL_command_queue), cmd_q, C_NULL)
        return CmdQueue(cmd_q[])
    end

    command_type(evt::CLEvent) = begin
        cmd_t = Ref{CL_int}()
        @check api.clGetEventInfo(evt.id, CL_EVENT_COMMAND_TYPE,
                                  sizeof(CL_int), cmd_t, C_NULL)
        return cmd_t[]
    end

    reference_count(evt::CLEvent) = begin
        cnt = Ref{CL_uint}()
        @check api.clGetEventInfo(evt.id, CL_EVENT_REFERENCE_COUNT,
                                  sizeof(CL_uint), cnt, C_NULL)
        return cnt[]
    end

    context(evt::CLEvent) = begin
        ctx = Ref{CL_context}()
        @check api.clGetEventInfo(evt.id, CL_EVENT_CONTEXT,
                                  sizeof(CL_context), CL_context, C_NULL)
        Context(ctx[])
    end

    status(evt::CLEvent) = begin
        st = Ref{CL_int}()
        @check api.clGetEventInfo(evt.id, CL_EVENT_COMMAND_EXECUTION_STATUS,
                                  sizeof(CL_int), st, C_NULL)
        status = st[]
        if status == CL_QUEUED
            return :queued
        elseif status == CL_SUBMITTED
            return :submitted
        elseif status == CL_RUNNING
            return :running
        elseif status == CL_COMPLETE
            return :complete
        else
            throw(ArgumentError("Unknown status value: $status"))
        end
    end

    @profile_info(profile_start,  CL_PROFILING_COMMAND_START)
    @profile_info(profile_end,    CL_PROFILING_COMMAND_END)
    @profile_info(profile_queued, CL_PROFILING_COMMAND_QUEUED)
    @profile_info(profile_submit, CL_PROFILING_COMMAND_SUBMIT)

    profile_duration(evt::Event) = begin
        return evt[:profile_end] - evt[:profile_start]
    end

    info_map = Dict{Symbol, Function}(
        :context => context,
        :command_queue => command_queue,
        :reference_count => reference_count,
        :command_type => command_type,
        :status => status,
        :profile_start => profile_start,
        :profile_end => profile_end,
        :profile_queued => profile_queued,
        :profile_submit => profile_submit,
        :profile_duration => profile_duration,
    )

    try
        func = info_map[evt_info]
        func(evt)
    catch err
        if isa(err, KeyError)
            throw(ArgumentError("OpenCL.Event has no info for: $evt_info"))
        else
            throw(err)
        end
    end
end


text:
module ComradeAHMC

using AbstractDifferentiation
using AbstractMCMC
using Reexport
@reexport using AdvancedHMC
using Comrade
using TupleVectors
using ArgCheck: @argcheck
using Random

export sample, AHMC

Base.@kwdef struct AHMC{S,I,P,T,A,D}
    metric::S
    integrator::I = Leapfrog
    trajectory::P = MultinomialTS
    termination::T = StrictGeneralisedNoUTurn(10, 1000.0)
    adaptor::A = StanHMCAdaptor
    targetacc::Float64 = 0.8
    autodiff::D = AD.ForwardDiffBackend()
end

Comrade.samplertype(::Type{<:AHMC}) = Comrade.IsFlat()

function _initialize_hmc(tpost::Comrade.TransformedPosterior, init_params, nchains)
    isnothing(init_params) && return Comrade.HypercubeTransform.inverse.(Ref(tpost.transform), rand(tpost.lpost.prior, nchains))
    @argcheck length(init_params) == nchains
    return init_params
end

function AbstractMCMC.sample(post::Comrade.Posterior,
    sampler::A, parallel::AbstractMCMC.AbstractMCMCEnsemble,
    nsamples, nchains;
    init_params=nothing, kwargs...
    ) where {A<:AHMC}
    tpost = asflat(post)
    if isnothing(init_params)
        θ0 = prior_sample(tpost, nchains)
    else
        θ0 = Comrade.HypercubeTransform.inverse.(Ref(tpost), init_params)
    end
    return Comrade.sample(tpost,
                            sampler,
                            parallel,
                            nsamples,
                            nchains;
                            init_params=θ0,
                            kwargs...)
end

function AbstractMCMC.sample(tpost::Comrade.TransformedPosterior,
                             sampler::AHMC, parallel::AbstractMCMC.AbstractMCMCEnsemble,
                             nsamples, nchains;
                             init_params=nothing, kwargs...
                             )

    ℓ = logdensityof(tpost)

    ∇ℓ = Comrade.make_pullback(ℓ, sampler.autodiff)
    θ0 = _initialize_hmc(tpost, init_params, nchains)
    model = AdvancedHMC.DifferentiableDensityModel(ℓ, ∇ℓ)
    metric = sampler.metric
    # This is a hack to get a good initial step size
    hamiltonian = Hamiltonian(metric, ℓ, ∇ℓ)
    ϵ0 = find_good_stepsize(hamiltonian, first(θ0))
    integrator = sampler.integrator(ϵ0)

    # form the HMCKernel
    kernel = HMCKernel(Trajectory{sampler.trajectory}(integrator, sampler.termination))
    adaptor = sampler.adaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(sampler.targetacc, integrator))

    res = AbstractMCMC.sample(Random.GLOBAL_RNG, model, kernel, metric, adaptor, parallel, nsamples, nchains; init_params=θ0, chain_type=Array, kwargs...)

    stats = [TupleVector(getproperty.(r, :stat)) for r in res]
    samples = [getproperty.(getproperty.(r, :z), :θ) for r in res]
    chains = [TupleVector(transform.(Ref(tpost), s)) for s in samples]
    return chains, stats

end

function AbstractMCMC.sample(tpost::Comrade.TransformedPosterior, sampler::AHMC, nsamples, args...;
                             init_params=nothing,
                             kwargs...)
    ℓ = logdensityof(tpost)

    ∇ℓ = Comrade.make_pullback(ℓ, sampler.autodiff)

    θ0 = init_params
    if isnothing(init_params)
        @warn "No starting location chosen, picking start from random"
        θ0 = prior_sample(post)
    end
    model = AdvancedHMC.DifferentiableDensityModel(ℓ, ∇ℓ)
    metric = sampler.metric
    # This is a hack to get a good initial step size
    hamiltonian = Hamiltonian(metric, ℓ, ∇ℓ)
    ϵ0 = find_good_stepsize(hamiltonian, θ0)
    integrator = sampler.integrator(ϵ0)

    # form the HMCKernel
    kernel = HMCKernel(Trajectory{sampler.trajectory}(integrator, sampler.termination))
    adaptor = sampler.adaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(sampler.targetacc, integrator))

    res = AbstractMCMC.sample(model, kernel, metric, adaptor, nsamples, args...; init_params=θ0, chain_type=Array, kwargs...)

    stats = TupleVector(getproperty.(res, :stat))
    samples = getproperty.(getproperty.(res, :z), :θ)
    chain = transform.(Ref(tpost), samples)
    return TupleVector(chain), stats
end

end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 44b34e3e-f591-11eb-33c7-0d192fda6f83
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using QuantumOptics
end

# ╔═╡ 07955d5f-1272-4074-a531-5bc868722a7d
begin
	using PyPlot
	ion()
end

# ╔═╡ 9f0d6f82-87cf-405f-9c79-233b2c39ea2e
begin
	const dim = 5

	fbasis = FockBasis(dim)
end

# ╔═╡ f34df059-eb22-4db5-a96b-3396695ec432
begin
	a = destroy(fbasis)
	n = number(fbasis)
end;

# ╔═╡ 37e774dc-ebca-4716-a44e-92bad8c91b7f
coh(α) = coherentstate(fbasis, α)

# ╔═╡ 1dd7ab9a-cc7d-4f2b-a4ce-9ee51a205fc0
dagger(coh(1))*dagger(a)*coh(1)

# ╔═╡ d0991d20-43b3-49b0-b579-6879815bc12b
sqee(ξ) = sparse(exp((1/2) * (conj(ξ) * dense(a)^2 - ξ * dense(dagger(a))^2)))

# ╔═╡ 02163fd0-99b0-492d-b03b-238890284567
diss(α) = sparse(exp(α * dense(dagger(a)) - conj(α) * dense(a)))

# ╔═╡ e86d36fe-28ff-49f3-bd7a-20fef0390695
acoh(ω, λ, γ, α, t) = (α + (λ/ω)) * exp(-γ * t * exp(- 1im * ω / γ)) - (λ/ω)

# ╔═╡ d3c34519-101e-4c6c-ba05-97684fde57b7
tlist = range(0, stop = 5, length = 500)

# ╔═╡ 3e84e3d9-c822-4451-bd8b-96c31e6e9ef6
tacoh = [acoh(4, 1, 5, 4, tlist[τ]) for τ in 1:length(tlist)]

# ╔═╡ 875caefb-fd1a-4c21-8749-03e5c2d8dcaf
begin
	plot(tlist, real.(tacoh))
	plot(tlist, imag.(tacoh))
end

# ╔═╡ Cell order:
# ╠═44b34e3e-f591-11eb-33c7-0d192fda6f83
# ╠═9f0d6f82-87cf-405f-9c79-233b2c39ea2e
# ╠═f34df059-eb22-4db5-a96b-3396695ec432
# ╠═37e774dc-ebca-4716-a44e-92bad8c91b7f
# ╠═1dd7ab9a-cc7d-4f2b-a4ce-9ee51a205fc0
# ╠═d0991d20-43b3-49b0-b579-6879815bc12b
# ╠═02163fd0-99b0-492d-b03b-238890284567
# ╠═e86d36fe-28ff-49f3-bd7a-20fef0390695
# ╠═07955d5f-1272-4074-a531-5bc868722a7d
# ╠═d3c34519-101e-4c6c-ba05-97684fde57b7
# ╠═3e84e3d9-c822-4451-bd8b-96c31e6e9ef6
# ╠═875caefb-fd1a-4c21-8749-03e5c2d8dcaf

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
# samplers
struct LogUniform <: Sampleable{Univariate,Continuous}
    d

    LogUniform(a, b) = new(Uniform(log10(a), log10(b)))
end

Base.rand(rng::AbstractRNG, s::LogUniform) = 10 ^ rand(rng, s.d)

# Batch partitioning with buffer
using Random: AbstractRNG, shuffle!, GLOBAL_RNG, randperm

struct BatchPartition{R<:AbstractRNG}
    indices::Vector{Int}
    n::Int
    batchsize::Int
    epochs::Int
    epochlength::Int
    batches::Int
    lastbatch::Vector{Int}
    shuffle::Bool
    include_buffer::Bool
    buffer
    rng::R
end

function BatchPartition(
    n;
    epochs = 1,
    batchsize = 1,
    lastbatch = rand(1:n, batchsize),
    shuffle = false,
    include_buffer = false,
    buffer = () -> Int[],
    rng = GLOBAL_RNG
)

    batchsize > 0 || throw(ArgumentError("batchsize must be positive integer"))
    batchsize <= n || throw(ArgumentError("number of observations less than batchsize"))

    return BatchPartition(
        collect(1:n),
        n,
        batchsize,
        epochs,
        n ÷ batchsize,
        epochs * (n ÷ batchsize),
        lastbatch,
        shuffle,
        include_buffer,
        buffer,
        rng
    )
end

function Base.show(io::IO, d::BatchPartition)
    println(io, "BatchPartition:")
    println(io, "  ⋅ batch size: ", d.batchsize)
    println(io, "  ⋅ epoch length: ", d.epochlength)
    println(io, "  ⋅ number of epochs: ", d.epochs, " (", d.batches, " batches)")
    println(io, "  ⋅ shuffle: ", d.shuffle)
    print(io,   "  ⋅ include buffer: ", d.include_buffer)
end

Base.length(d::BatchPartition) = d.batches
Base.eltype(::BatchPartition) = Vector{Int}

function Base.iterate(d::BatchPartition, batch_ind = 1)
    if batch_ind > d.batches
        return nothing
    end
    if batch_ind % d.epochlength == 1 && d.shuffle
        shuffle!(d.rng, d.indices)
    end

    inds = collect(1:d.batchsize) .+ rem(batch_ind - 1, d.epochlength)*d.batchsize
    d.include_buffer && update_indices!(d, inds, d.buffer())
    batch = d.indices[inds]
    d.lastbatch .= batch
    return (batch, batch_ind + 1)
end

function update_indices!(::BatchPartition, inds, b_inds)
    @warn("indices provided by buffer muset be of type Int or Vector{Int}")
    return
end

function update_indices!(d::BatchPartition, inds, b_inds::Union{Int, Vector{Int}})
    if !all(0 .< b_inds .<= d.batchsize)
        @warn("indices in buffer are corrupted")
        return
    end
    k = d.lastbatch[b_inds]
    if !all(0 .< k .<= d.n)
        @warn("last batch is corrupted")
        return
    end

    pos = randperm(d.batchsize)[1:length(k)]
    inds[pos] .= k
    return
end


text:
module BigCombinatorics

export Fibonacci
export Factorial, DoubleFactorial, FallingFactorial, RisingFactorial, HyperFactorial
export Binomial, Catalan
export Derangements, MultiChoose, Multinomial
export Bell, Stirling1, Stirling2
export IntPartitions, IntPartitionsDistinct
export Euler, PowerSum, Menage

"""
This is where we cache values already computed
"""
_master_table = Dict{Function,Dict}()

"""
    _save(f::Function, x, val::BigInt)

Save the computed value of `f` with argument(s) `x` and value `val`.
"""
function _save(f::Function, x, val::BigInt)
    d = _master_table[f]
    d[x] = val
    nothing
end

"""
    _has(f::Function, x)

Check if argument `x` for function `f` is already in the `_master_table`.
"""
function _has(f::Function, x)::Bool
    d = _master_table[f]
    return haskey(d, x)
end

"""
    _max_arg(f::Function)

Find the largest argument known for `f`. Note that `f` must be a function
of just a single value.
"""
function _max_arg(f::Function)
    tab = _master_table[f]
    return maximum(keys(tab))
end

"""
    _get(f::Function, x)

Retrieve the function value `f(x)` from the `_master_table`.
"""
function _get(f::Function, x)::BigInt
    d = _master_table[f]
    return d[x]
end

"""
    _make(f::Function, T::Type)

Create an entry in the `_master_table` for the function `f`.
"""
function _make(f::Function, T::Type)
    _master_table[f] = Dict{T,BigInt}()
    nothing
end

"""
`BigCombinatorics.cache_report()` reports how many
entries are saved for each function in the `BigCombinatorics`
module.
"""
function cache_report()
    total = 0
    for func in keys(_master_table)
        parse_func = split(string(func), ".")
        func_name = last(parse_func)
        sz = length(_master_table[func])
        total += sz
        println("$sz\t$func_name")
    end
    println()
    println("$total\tTotal entries")
    nothing
end

"""
    BigCombinatorics.cache_clear()

Clears all cached values.
"""
function cache_clear()
    for f in keys(_master_table)
        f()
    end
    nothing
end

"""
`Fibonacci(n)` returns the `n`-th Fibonacci number.
We begin with `Fibonacci(0)==0` and `Fibonacci(1)==1`.
See https://oeis.org/A000045
"""
function Fibonacci(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonngative"))
    end

    if _has(Fibonacci, n)
        return _get(Fibonacci, n)
    end

    start = _max_arg(Fibonacci) + 1
    for m = start:n
        val = _get(Fibonacci, m - 1) + _get(Fibonacci, m - 2)
        _save(Fibonacci, m, val)
    end

    return _get(Fibonacci, n)
end

function Fibonacci()
    _make(Fibonacci, Integer)
    _save(Fibonacci, 0, big(0))
    _save(Fibonacci, 1, big(1))
end
Fibonacci()

"""
`Factorial(n)` returns `n!` for nonnegative integers `n`.
See https://oeis.org/A000142

`Factorial(n,k)` returns `n!/k!` (to be consistent with Julia's
`factorial`.) Requires `0 <= k <= n`.

See also `FallingFactorial` and `RisingFactorial`.
"""
function Factorial(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonngative"))
    end
    return factorial(big(n))
end

function Factorial(n::Integer, k::Integer)::BigInt
    if k > n
        throw(DomainError((n, k), "$k cannot exceed $n"))
    end
    return div(Factorial(n), Factorial(k))
end

"""
    FallingFactorial(n,k)

returns `n*(n-1)*(n-2)*...*(n-k+1)`
(with a total of `k` factors). Requires `n,k >= 0`.
If `k>n` then `0` is returned.
"""
function FallingFactorial(n::Integer, k::Integer)::BigInt
    if n < 0 || k < 0
        throw(DomainError((n, k), " arguments must be nonnegative"))
    end
    if k > n
        return big(0)
    end
    return Factorial(n, n - k)
end

"""
    RisingFactorial(n,k)

returns `n*(n+1)*(n+2)*...*(n+k-1)`
(with a total of `k` factors). Requires `n,k >= 0`.
"""
function RisingFactorial(n::Integer, k::Integer)::BigInt
    if n < 0 || k < 0
        throw(DomainError((n, k), "arguments must be nonnegative"))
    end
    if k == 0
        return big(1)
    end
    if n == 0
        return big(0)
    end

    return FallingFactorial(n + k - 1, k)
end

"""
    DoubleFactorial(n)

returns `n!!`, i.e., `n*(n-2)*...` with `(-1)!! == 0!! == 1!! == 1`.
See https://oeis.org/A006882
"""
function DoubleFactorial(n::Integer)::BigInt
    if n < -1
        throw(DomainError(n, "argument must be at least -1"))
    end
    if _has(DoubleFactorial, n)
        return _get(DoubleFactorial, n)
    end

    start = _max_arg(DoubleFactorial) + 1
    for m = start:n
        val = m * _get(DoubleFactorial, m - 2)
        _save(DoubleFactorial, m, val)
    end

    return _get(DoubleFactorial, n)
end

function DoubleFactorial()
    _make(DoubleFactorial, Integer)
    _save(DoubleFactorial, -1, big(1))
    _save(DoubleFactorial, 0, big(1))
    _save(DoubleFactorial, 1, big(1))
end
DoubleFactorial()

"""
    HyperFactorial(n)

returns the hyperfactorial of `n`, that is
`1^1 * 2^2 * 3^3 * ... * n^n`.
See https://oeis.org/A002109
"""
function HyperFactorial(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "arument must be nonnegative"))
    end

    if _has(HyperFactorial, n)
        return _get(HyperFactorial, n)
    end
    start = _max_arg(HyperFactorial) + 1
    for m = start:n
        val = (big(m))^m * _get(HyperFactorial, m - 1)
        _save(HyperFactorial, m, val)
    end
    return _get(HyperFactorial, n)
end

function HyperFactorial()
    _make(HyperFactorial, Integer)
    _save(HyperFactorial, 0, big(1))
    _save(HyperFactorial, 1, big(1))
end
HyperFactorial()

"""
    Binomial(n,k)

returns the binomial coefficient `n`-choose-`k`.
This is the number of `k`-element subsets of an `n`-element set.
See https://oeis.org/A007318
"""
Binomial(n::Integer, k::Integer) = binomial(big(n), big(k))::BigInt

"""
    Multinomial(vec)

returns the multinomial coefficient whose
top index is the sum of `vec` (an array of `Int`s) and whose
bottom indices are given by `vec`.

This may also be called with a common-separated list of arguments,
that is, either of `Multinomial([1,2,3])` or `Multinomial(1,2,3)`.
The result is `60` in both cases as these equal `6!/(1! 2! 3!)`.

**Warning**: This is not the same as `MultiChoose`.
"""
function Multinomial(v...)::BigInt
    nv = length(v)
    for i = 1:nv
        typeof(v[i]) <: Integer || throw(DomainError(v, "arguments must be integers"))
        v[i] >= 0 || throw(DomainError(v, "arguments must be nonngative"))
    end
    vals = [t for t in v]
    return Multinomial(vals)
end

Multinomial() = big(1)::BigInt

function Multinomial(vals::Vector{T})::BigInt where {T<:Integer}
    if any([t < 0 for t in vals])
        throw(DomainError(vals, "arguments must be nonnegative"))
    end

    nv = length(vals)
    n = sum(vals)
    # base cases
    if nv <= 1 || n == 0
        return big(1)
    end
    # reduce
    return Binomial(n, vals[end]) * Multinomial(vals[1:nv-1])
end

"""
    MultiChoose(n,k)

returns the number of `k`-element
*multisets* that can be formed using the elements of an
`n`-element set.

**Warning**: This is not the same as `Multinomial`.
"""
function MultiChoose(n::Integer, k::Integer)::BigInt
    return Binomial(n + k - 1, k)
end

"""
    Catalan(n)

returns the `n`-th Catalan number.
See https://oeis.org/A000108
"""
function Catalan(n::Integer)::BigInt
    n >= 0 || throw(DomainError(n, "argument must be nonnegative"))
    return div(Binomial(2n, n), n + 1)
end

"""
    Derangements(n)

returns the number of permutations of
an `n`-set that have no fixed point.
See https://oeis.org/A000166
"""
function Derangements(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonnegative"))
    end
    if _has(Derangements, n)
        return _get(Derangements, n)
    end

    start = _max_arg(Derangements) + 1
    for m = start:n
        s = (m % 2 == 0) ? 1 : -1
        val = m * _get(Derangements, m - 1) + s
        _save(Derangements, m, val)
    end
    return _get(Derangements, n)
end

function Derangements()
    _make(Derangements, Integer)
    _save(Derangements, 0, big(1))
    _save(Derangements, 1, big(0))
end
Derangements()

"""
    Bell(n)

gives the `n`-th Bell number, that is,
the number of partitions of an `n`-element set.
See https://oeis.org/A000110
"""
function Bell(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonnegative"))
    end

    if _has(Bell, n)
        return _get(Bell, n)
    end

    start = _max_arg(Bell) + 1
    for m = start:n
        val = big(0)
        for k = 0:m-1
            val += Binomial(m - 1, k) * Bell(k)
        end
        _save(Bell, m, val)
    end

    return _get(Bell, n)
end

function Bell()
    _make(Bell, Integer)
    _save(Bell, 0, big(1))
    _save(Bell, 1, big(1))
end
Bell()

"""
    Stirling2(n,k)

gives the Stirling number of the second kind,
that is, the number of paritions of an `n`-set into `k`-parts."
See https://oeis.org/A008277
"""
function Stirling2(n::Integer, k::Integer)::BigInt
    # special cases
    if k < 0 || n < 0
        throw(DomainError((n, k), "arguments must be nonnegative"))
    end

    if k > n
        return big(0)
    end

    if n == 0  # and by logic, k==0
        return big(1)
    end

    if k == 0
        return big(0)
    end

    if n == k
        return big(1)
    end
    # END OF SPECIAL CASES, invoke recursion
    if _has(Stirling2, (n, k))
        return _get(Stirling2, (n, k))
    end

    val = Stirling2(n - 1, k - 1) + Stirling2(n - 1, k) * k
    _save(Stirling2, (n, k), val)
    return val
end
function Stirling2()
    _make(Stirling2, Tuple{Integer,Integer})
end
Stirling2()

"""
    Stirling1(n,k)

gives the (signed) Stirling number
of the first kind, that is, the coefficient of `x^k`
in the poynomial `x(x-1)(x-2)...(x-n+1)`.
See https://oeis.org/A008275
"""
function Stirling1(n::Integer, k::Integer)::BigInt
    # special cases
    if k < 0 || n < 0
        throw(DomainError((n, k), "arguments must be nonnegative"))
    end

    if k > n
        return big(0)
    end

    if n == 0  # and, by logic, k==0
        return big(1)
    end

    if k == 0  # and, by logic, n>0
        return big(0)
    end

    if _has(Stirling1, (n, k))
        return _get(Stirling1, (n, k))
    end

    # end of special cases, invoke recursion

    val = Stirling1(n - 1, k - 1) - (n - 1) * Stirling1(n - 1, k)
    _save(Stirling1, (n, k), val)
    return val
end
function Stirling1()
    _make(Stirling1, Tuple{Integer,Integer})
end
Stirling1()

"""
`IntPartitions(n)` is the number of partitions of the integer `n`.
See https://oeis.org/A000041

`IntPartitions(n,k)` is the number of partitions of the integer
`n` with exactly `k` (nonzero) parts.
"""
function IntPartitions(n::Integer, k::Integer)::BigInt
    if n < 0 || k < 0
        throw(DomainError((n, k), "arguments must be nonnegative"))
    end
    # lots of special cases
    if k > n
        return big(0)
    end
    if n == 0
        return big(1)
    end
    if k == 0
        return big(0)
    end
    if k == n || k == 1
        return big(1)
    end
    if _has(IntPartitions, (n, k))
        return _get(IntPartitions, (n, k))
    end

    val = sum([IntPartitions(n - k, i) for i = 0:k])
    _save(IntPartitions, (n, k), val)
    return val
end

function IntPartitions(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonnegative"))
    end
    if _has(IntPartitions, n)
        return _get(IntPartitions, n)
    end
    val = sum([IntPartitions(n, k) for k = 0:n])
    _save(IntPartitions, n, val)
    return val
end
function IntPartitions()
    _make(IntPartitions, Union{Tuple{Integer,Integer},Integer})
end
IntPartitions()

"""
`IntPartitionsDistinct(n,k)` is the number of partitions of
the integer `n` into exactly `k` *distinct* parts.

`IntPartitionsDistinct(n)` is the number of partitions of `n`
into *distinct* parts.
"""
function IntPartitionsDistinct(n::Integer, k::Integer)::BigInt
    if n < 0 || k < 0
        throw(DomainError((n, k), "arguments must be nonnegative"))
    end
    Ck2 = div(k * (k - 1), 2)
    if n < Ck2
        return big(0)
    end
    return IntPartitions(n - Ck2, k)
end

function IntPartitionsDistinct(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n, "argument must be nonngative"))
    end
    result = big(0)
    for k = 1:n
        s = IntPartitionsDistinct(n, k)
        if s == 0
            break
        end
        result += s
    end
    return result
end

"""
    Euler(n)

returns the `n`-th Euler number. Starting with `n=0`
this is the sequence
1, 0, -1, 0, 5, 0, -61, 0, 1385 and so on.
See https://oeis.org/A122045

Not to be confused with `Eulerian`.
"""
function Euler(n::Integer)::BigInt
    n >= 0 || throw(DomainError(n, "argument must be nonngative"))
    if n % 2 == 1
        return big(0)
    end

    if _has(Euler, n)
        return _get(Euler, n)
    end

    start = _max_arg(Euler) + 2

    for m = start:2:n
        last = div(m, 2) - 1
        val = -sum([Binomial(m, 2k) * Euler(2k) for k = 0:last])
        _save(Euler, m, val)
    end

    return _get(Euler, n)
end

function Euler()
    _make(Euler, Integer)
    _save(Euler, 0, big(1))
end
Euler()

"""
    PowerSum(n,k)

returns the sum of the `k`-th powers of the
integers `1` through `n`, i.e.,
`1^k + 2^k + 3^k + ... + n^k`.
"""
function PowerSum(n::Integer, k::Integer)::BigInt
    (n >= 0 && k >= 0) || throw(DomainError((n, k), "arguments must be nonngative"))
    # Base and special cases
    if n == 0
        return big(0)
    end
    if k == 0
        return big(n)
    end
    if k == 1
        return Binomial(n + 1, 2)
    end
    if _has(PowerSum, (n, k))
        return _get(PowerSum, (n, k))
    end

    val = big(n)^k + PowerSum(n - 1, k)
    _save(PowerSum, (n, k), val)
    return val
end
function PowerSum()
    _make(PowerSum, Tuple{Integer,Integer})
end
PowerSum()

export Eulerian
"""
    Eulerian(n,k)

returns the number of permutations of `{1,2,...,n}`
with `k` ascents. See https://oeis.org/A008292

Not to be confused with `Euler`.
"""
function Eulerian(n::Integer, k::Integer)::BigInt
    @assert (n >= 0 && k >= 0) "$n,$k must both be nonnegative"

    if n == 0
        return big(0)
    end

    if k > n || k == 0
        return big(0)
    end

    if n < 2
        return big(1)
    end

    if k == n      # includes the case n=k=0
        return big(1)
    end

    if k == 1
        return big(1)
    end

    if _has(Eulerian, (n, k))
        return _get(Eulerian, (n, k))
    end

    val = (n - k + 1) * Eulerian(n - 1, k - 1) + k * Eulerian(n - 1, k)
    _save(Eulerian, (n, k), val)
    return val
end
function Eulerian()
    _make(Eulerian, Tuple{Integer,Integer})
end
Eulerian()

"""
    Menage(n)

Number of solutions to the Menage problem with `n` male/female couples.
See https://oeis.org/A059375
"""
function Menage(n::Int)::BigInt
    if n < 0
        throw(DomainError(n, "Argument must be nonnegative"))
    end
    if _has(Menage, n)
        return _get(Menage, n)
    end
    start = _max_arg(Menage) + 1
    for m = start:n
        s = (m % 2 == 0) ? -1 : 1
        a1 = _get(Menage, m - 1) ÷ (2 * Factorial(m - 1))
        a2 = _get(Menage, m - 2) ÷ (2 * Factorial(m - 2))

        am = m * a1 + (m * a2 + 4s) ÷ (m - 2)
        val = 2 * Factorial(m) * am

        _save(Menage, m, val)
    end
    return _get(Menage, n)
end
function Menage()
    _make(Menage, Int)
    _save(Menage, 0, big(1))
    _save(Menage, 1, big(0))
    _save(Menage, 2, big(0))
    _save(Menage, 3, big(12))
end
Menage()

end  #end of module


text:
include("../mathutils/Helpers.jl")
include("../termstructures/YieldCurve.jl")
include("../models/StochasticProcess.jl")

using Roots

# the actual data structure for the model
struct HullWhiteModel{T<:AbstractFloat} <: StochasticProcess
    yieldCurve::YieldCurve
    meanReversion::T
    volatilityTimes::Array{T}
    volatilityValues::Array{T}
    y_::Array{T}
end

# constructor

function HullWhiteModel(yieldCurve, meanReversion, volatilityTimes::Array, volatilityValues::Array)
    #
    t0 = 0.0
    y0 = 0.0
    y_ = zeros(size(volatilityTimes))
    for i = 1:size(y_)[1]
        y_[i] = (GPrime(meanReversion,t0,volatilityTimes[i])^2) * y0 +
            (volatilityValues[i]^2) *
            (1.0 - exp(-2*meanReversion*(volatilityTimes[i]-t0))) /
            (2.0 * meanReversion)
        t0 = volatilityTimes[i]
        y0 = y_[i]
    end
    return HullWhiteModel(yieldCurve,meanReversion,volatilityTimes,volatilityValues,y_)
end

# auxilliary functions/methods

G(meanReversion,t,T) =
    (1.0 - exp(-meanReversion*(T-t))) / meanReversion

GPrime(meanReversion,t,T) = exp(-meanReversion*(T-t))

G(self::HullWhiteModel,t,T) = G(self.meanReversion,t,T)

GPrime(self::HullWhiteModel,t,T) = GPrime(self.meanReversion,t,T)

function y(self::HullWhiteModel,t)
    idx = searchsortedfirst(self.volatilityTimes,t)  # first element of tupel
    t0 = (idx==1) ? 0.0 : self.volatilityTimes[idx-1]
    y0 = (idx==1) ? 0.0 : self.y_[idx-1]
    s1 = self.volatilityValues[min(idx,size(self.volatilityValues)[1])]  # flat extrapolation
    y1 = (GPrime(self,t0,t)^2) * y0 +
            s1^2 * (1.0 - exp(-2*self.meanReversion*(t-t0))) /
            (2.0 * self.meanReversion)
    return y1
end

function riskNeutralExpectationX(self::HullWhiteModel, t, xt, T)
    # E[x] = G'(t,T)x + \int_t^T G'(u,T)y(u)du
    f(u) = GPrime(self,u,T)*y(self,u)
    # use Simpson's rule to approximate integral, this should better be solved analytically
    integral = (T-t) / 6 * (f(t) + 4*f((t+T)/2) + f(T))
    return GPrime(self,t,T)*xt + integral
end

function sigma(self::HullWhiteModel,t)   # Todo test this method
    # find idx s.t. t[idx-1] < t <= t[idx]
    idx = searchsortedfirst(self.volatilityTimes,t)
    return self.volatilityValues[min(idx,size(self.volatilityValues)[1])]
end

# bond volatility is used in hybrid model vol adjuster

function zeroBondVolatility(self::HullWhiteModel, t, T)
    # we wrap scalar bond volatility into array to allow
    # for generalisation to multi-factor models
    return [ sigma(self,t)*G(self,t,T) ]
end

function zeroBondVolatilityPrime(self::HullWhiteModel, t, T)
    # we wrap scalar bond volatility into array to allow
    # for generalisation to multi-factor models
    return [ sigma(self,t)*GPrime(self,t,T) ]
end

# model methods

# conditional expectation in T-forward measure
function expectationX(self::HullWhiteModel, t, xt, T)
    return GPrime(self,t,T)*(xt + G(self,t,T)*y(self,t))
end

function varianceX(self::HullWhiteModel, t, T)
    return y(self,T) - GPrime(self,t,T)^2 * y(self,t)
end

function zeroBondPrice(self::HullWhiteModel, t, T, xt)
    G_ = G(self,t,T)
    return discount(self.yieldCurve,T) / discount(self.yieldCurve,t) *
        exp(-G_*xt - 0.5 * G_^2 * y(self,t) )
end

function zeroBondOption(self::HullWhiteModel, expiryTime, maturityTime, strikePrice, callOrPut)
    nu2 = G(self,expiryTime,maturityTime)^2 * y(self,expiryTime)
    P0  = discount(self.yieldCurve, expiryTime)
    P1  = discount(self.yieldCurve, maturityTime)
    return P0 * Black(strikePrice,P1/P0,sqrt(nu2),1.0,callOrPut)
end

function couponBondOption(self::HullWhiteModel, expiryTime, payTimes, cashFlows, strikePrice, callOrPut)
    function objective(x)
        bond = 0
        for i = 1:size(payTimes)[1]
            bond += cashFlows[i] * zeroBondPrice(self,expiryTime,payTimes[i],x)
        end
        return bond - strikePrice
    end
    xStar = find_zero(objective,(-1.0,1.0),Roots.Brent(),xatol=1.0e-8)
    bondOption = 0.0
    for i = 1:size(payTimes)[1]
        strike = zeroBondPrice(self,expiryTime,payTimes[i],xStar)
        bondOption += cashFlows[i] * zeroBondOption(self,expiryTime,payTimes[i],strike,callOrPut)
    end
    return bondOption
end

# stochastic process interface

function stateSize(self::HullWhiteModel)   # dimension of X(t)
    return 2      # [x, s], we also need for numeraire s = \int_0^t r dt
end

function factors(self::HullWhiteModel)   # dimension of W(t)
    return 1
end

function initialValues(self::HullWhiteModel)
    return [ 0.0, 0.0 ]
end

function zeroBond(self::HullWhiteModel, t, T, X, alias)
    return zeroBondPrice(self, t, T, X[1])
end

function numeraire(self::HullWhiteModel, t, X)
    return exp(X[2]) / discount(self.yieldCurve,t)
end

# evolve X(t0) -> X(t0+dt) using independent Brownian increments dW
# t0, dt are assumed float, X0, X1, dW are np.array
function evolve(self::HullWhiteModel, t0, X0, dt, dW, X1)
    x1 = riskNeutralExpectationX(self,t0,X0[1],t0+dt)
    # x1 = X0[0] + (self.y(t0) - self.meanReversion*X0[0])*dt
    nu = sqrt(varianceX(self,t0,t0+dt))
    x1 = x1 + nu*dW[1]
    # s1 = s0 + \int_t0^t0+dt x dt via Trapezoidal rule
    s1 = X0[2] + (X0[1] + x1) * dt / 2
    # gather results
    X1[1] = x1
    X1[2] = s1
    return nothing
end

# the short rate over an integration time period
# this is required for drift calculation in multi-asset and hybrid models
function shortRateOverPeriod(self::HullWhiteModel, t0, dt, X0, X1)
    B_d = discount(self.yieldCurve,t0) / discount(self.yieldCurve,t0 + dt)  # deterministic drift part for r_d
    x_av = 0.5 * (X0[1] + X1[1])
    return log(B_d) / dt + x_av
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Classify buildings ###
function classify(buildings::Dict{Int,Building})
    bdgs = Dict{Int,Int}()

    for (key, building) in buildings
        if haskey(BUILDING_CLASSES, building.class)
            bdgs[key] = BUILDING_CLASSES[building.class]
        end
    end

    return bdgs
end


text:
@testset "read_surfaces.jl" begin

end


text:
export initsoln, easy_sdp

function initsoln(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix})
    m = length(As)
    X, y, Z = initsoln(Cint(size(C, 1)), Cint(m), C.csdp, fptr(b), fptr(As))
    mywrap(X), mywrap(y, m), mywrap(Z)
end
function initsoln(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix})
    initsoln(C, b, [A.csdp for A in As])
end

function easy_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, verbose=false)
    # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are use in the code
    # so no need to worry, they won't change :)
    Xcsdp = X.csdp
    ycsdp = fptr(y)
    Zcsdp = Z.csdp
    # https://thenewphalls.wordpress.com/2014/03/21/capturing-output-in-julia/
    if !verbose
        oldstdout = STDOUT
        rd, wr = redirect_stdout()
    end
    status, pobj, dobj = easy_sdp(Cint(size(C, 1)), Cint(length(As)), C.csdp, fptr(b), fptr(As), 0.0, ptr(Xcsdp), ptr(ycsdp), ptr(Zcsdp))
    if !verbose
        redirect_stdout(oldstdout)
    end
    # I can even assert that they won't change
    @assert Xcsdp == X.csdp
    @assert ycsdp == fptr(y)
    @assert Zcsdp == Z.csdp
    status, pobj, dobj
end

function write_prob(fname::String, C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix})
    write_prob(fname, Cint(size(C, 1)), Cint(length(As)), C.csdp, fptr(b), fptr(As))
end

function write_sol(fname::String, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix)
    write_sol(fname, Cint(size(X, 1)), Cint(length(y)), X.csdp, fptr(y), Z.csdp)
end


text:
"""
A matrix of colors is interpreted as an image
"""
_default(::Node{Array{RGBA{N0f8}, 2}}, ::Style{:default}, ::Dict{Symbol,Any})

function _default(main::MatTypes{T}, ::Style, data::Dict) where T <: Colorant
    @gen_defaults! data begin
        spatialorder = "yx"
    end
    if !(spatialorder in ("xy", "yx"))
        error("Spatial order only accepts \"xy\" or \"yz\" as a value. Found: $spatialorder")
    end
    ranges = get(data, :ranges) do
        const_lift(main, spatialorder) do m, s
            (0:size(m, s == "xy" ? 1 : 2), 0:size(m, s == "xy" ? 2 : 1))
        end
    end
    delete!(data, :ranges)
    @gen_defaults! data begin
        image = main => (Texture, "image, can be a Texture or Array of colors")
        primitive::GLUVMesh2D = const_lift(ranges) do r
            x, y = minimum(r[1]), minimum(r[2])
            xmax, ymax = maximum(r[1]), maximum(r[2])
            SimpleRectangle{Float32}(x, y, xmax - x, ymax - y)
        end
        preferred_camera      = :orthographic_pixel
        fxaa                  = false
        shader                = GLVisualizeShader(
            "fragment_output.frag", "uv_vert.vert", "texture.frag",
            view = Dict("uv_swizzle" => "o_uv.$(spatialorder)")
        )
    end
end
function _default(main::VectorTypes{T}, ::Style, data::Dict) where T <: Colorant
    @gen_defaults! data begin
        image                 = main => (Texture, "image, can be a Texture or Array of colors")
        primitive::GLUVMesh2D = SimpleRectangle{Float32}(0f0, 0f0, length(to_value(main)), 50f0) => "the 2D mesh the image is mapped to. Can be a 2D Geometry or mesh"
        preferred_camera      = :orthographic_pixel
        fxaa                  = false
        shader                = GLVisualizeShader(
            "fragment_output.frag", "uv_vert.vert", "texture.frag",
            view = Dict("uv_swizzle" => "o_uv.xy")
        )
    end
end

"""
A matrix of Intensities will result in a contourf kind of plot
"""
function gl_heatmap(main::MatTypes{T}, data::Dict) where T <: AbstractFloat
    main_v = to_value(main)
    @gen_defaults! data begin
        ranges = (0:size(main_v, 1), 0:size(main_v, 2))
    end
    prim = const_lift(data[:ranges]) do ranges
        x, y, xw, yh = minimum(ranges[1]), minimum(ranges[2]), maximum(ranges[1]), maximum(ranges[2])
        SimpleRectangle{Float32}(x, y, xw-x, yh-y)
    end
    delete!(data, :ranges)
    @gen_defaults! data begin
        intensity             = main => Texture
        color_map             = default(Vector{RGBA{N0f8}},s) => Texture
        primitive::GLUVMesh2D = prim
        color_norm            = const_lift(extrema2f0, main)
        stroke_width::Float32 = 0.05f0
        levels::Float32       = 5f0
        stroke_color          = RGBA{Float32}(1,1,1,1)
        shader                = GLVisualizeShader("fragment_output.frag", "uv_vert.vert", "intensity.frag")
        fxaa                  = false
    end
end

"""
Float matrix with the style distancefield will be interpreted as a distancefield.
A distancefield is describing a shape, with positive values denoting the inside
of the shape, negative values the outside and 0 the border
"""
function _default(main::MatTypes{T}, s::style"distancefield", data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        distancefield = main => Texture
        shape         = DISTANCEFIELD
        fxaa          = false
    end
    rect = SimpleRectangle{Float32}(0f0,0f0, size(to_value(main))...)
    _default((rect, Point2f0[0]), s, data)
end

export play
"""
    play(img, timedim, t)

Slice a 3D array along axis `timedim` at time `t`.
This can be used to treat a 3D array like a video and create an image stream from it.
"""
function play(array::Array{T, 3}, timedim::Integer, t::Integer) where T
    index = ntuple(dim-> dim == timedim ? t : Colon(), Val(3))
    array[index...]
end

"""
    play(buffer, video_stream, t)

Plays a video stream from VideoIO.jl. You need to supply the image `buffer`,
which will be reused for better performance.
"""
function play(buffer::Array{T, 2}, video_stream, t) where T
    eof(video_stream) && seekstart(video_stream)
    w, h = size(buffer)
    buffer = reinterpret(UInt8, buffer, (3, w,h))
    read!(video_stream, buffer) # looses type and shape
    return reinterpret(T, buffer, (w,h))
end

unwrap(img::AxisArrays.AxisArray) = unwrap(img.data)
unwrap(img::AbstractArray) = img

GLAbstraction.gl_convert(::Type{T}, img::AbstractArray) where {T} = _gl_convert(T, unwrap(img))
_gl_convert(::Type{T}, img::Array) where {T} = gl_convert(T, img)

"""
    play(img)

Turns an Image into a video stream
"""
function play(img::HasAxesArray{T, 3}) where T
    ax = ImageAxes.timeaxis(img)
    if ImageAxes.timeaxis(img) != nothing
        return const_lift(play, unwrap(img), ImageAxes.timedim(img), loop(1:length(ax)))
    end
    error("Image has no time axis: axes(img) = $(axes(img))")
end

"""
Takes a 3D image and decides if it is a volume or an animated Image.
"""
function _default(img::HasAxesArray{T, 3}, s::Style, data::Dict) where T
    # We could do this as a @traitfn, except that those don't
    # currently mix well with non-trait specialization.
    if ImageAxes.timeaxis(img) != nothing
        data[:spatialorder] = "yx"
        td = ImageAxes.timedim(img)
        video_signal = const_lift(play, unwrap(img), td, loop(1:size(img, timedim)))
        return _default(video_signal, s, data)
    else
        ps = ImageAxes.pixelspacing(img)
        spacing = Vec3f0(map(x-> x / maximum(ps), ps))
        pdims   = Vec3f0(map(length, indices(img)))
        dims    = pdims .* spacing
        dims    = dims/maximum(dims)
        data[:dimensions] = dims
        _default(unwrap(img), s, data)
    end
end
function _default(img::TOrSignal{T}, s::Style, data::Dict) where T <: AxisMatrix
    @gen_defaults! data begin
        ranges = const_lift(img) do img
            ps = ImageAxes.pixelspacing(img)
            spacing = Vec2f0(map(x-> x / maximum(ps), ps))
            pdims   = Vec2f0(map(length, indices(img)))
            dims    = pdims .* spacing
            dims    = dims / maximum(dims)
            (0:dims[1], 0:dims[2])
        end
    end
    _default(const_lift(unwrap, img), s, data)
end

"""
Displays 3D array as movie with 3rd dimension as time dimension
"""
function _default(img::AbstractArray{T, 3}, s::Style, data::Dict) where T
    video_signal = const_lift(play, unwrap(img), 3, loop(1:size(img, 3)))
    return _default(video_signal, s, data)
end

"""
Takes a shader as a parametric function. The shader should contain a function stubb
like this:
```GLSL
uniform float arg1; // you can add arbitrary uniforms and supply them via the keyword args
float function(float x) {
 return arg1*sin(1/tan(x));
}
```
"""
_default(func::String, s::Style{:shader}, data::Dict) = @gen_defaults! data begin
    color                 = default(RGBA, s) => Texture
    dimensions            = (120f0, 120f0)
    primitive::GLUVMesh2D = SimpleRectangle{Float32}(0f0,0f0, dimensions...)
    preferred_camera      = :orthographic_pixel
    fxaa                  = false
    shader                = GLVisualizeShader(
        "fragment_output.frag", "parametric.vert", "parametric.frag",
        view = Dict("function" => func)
    )
end

#Volumes
const VolumeElTypes = Union{Gray, AbstractFloat, Intensity}

const default_style = Style{:default}()

function _default(a::VolumeTypes{T}, s::Style{:iso}, data::Dict) where T <: VolumeElTypes
    data = @gen_defaults! data begin
        isovalue  = 0.5f0
        algorithm = IsoValue
    end
     _default(a, default_style, data)
end

function _default(a::VolumeTypes{T}, s::Style{:absorption}, data::Dict) where T<:VolumeElTypes
    data = @gen_defaults! data begin
        absorption = 1f0
        algorithm  = Absorption
    end
    _default(a, default_style, data)
end

function _default(a::VolumeTypes{T}, s::Style{:absorption}, data::Dict) where T<:RGBA
    data = @gen_defaults! data begin
        algorithm  = AbsorptionRGBA
    end
    _default(a, default_style, data)
end

using .GLAbstraction: StandardPrerender

struct VolumePrerender
    sp::StandardPrerender
end
VolumePrerender(a, b) = VolumePrerender(StandardPrerender(a, b))

function (x::VolumePrerender)()
    x.sp()
    glEnable(GL_CULL_FACE)
    glCullFace(GL_FRONT)
end

function _default(main::VolumeTypes{T}, s::Style, data::Dict) where T <: VolumeElTypes
    @gen_defaults! data begin
        volumedata       = main => Texture
        hull::GLUVWMesh  = AABB{Float32}(Vec3f0(0), Vec3f0(1))
        light_position   = Vec3f0(0.25, 1.0, 3.0)
        light_intensity  = Vec3f0(15.0)
        model            = Mat4f0(I)
        modelinv         = const_lift(inv, model)

        color_map        = default(Vector{RGBA}, s) => Texture
        color_norm       = color_map == nothing ? nothing : const_lift(extrema2f0, main)
        color            = color_map == nothing ? default(RGBA, s) : nothing

        algorithm        = MaximumIntensityProjection
        absorption       = 1f0
        isovalue         = 0.5f0
        isorange         = 0.01f0
        shader           = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
        prerender        = VolumePrerender(data[:transparency], data[:overdraw])
        postrender       = () -> begin
            glDisable(GL_CULL_FACE)
        end
    end
end

function _default(main::VolumeTypes{T}, s::Style, data::Dict) where T <: RGBA
    @gen_defaults! data begin
        volumedata       = main => Texture
        hull::GLUVWMesh  = AABB{Float32}(Vec3f0(0), Vec3f0(1))
        model            = Mat4f0(I)
        modelinv         = const_lift(inv, model)

        # These don't do anything but are needed for type specification in the frag shader
        color_map        = nothing
        color_norm       = nothing
        color            = color_map == nothing ? default(RGBA, s) : nothing

        algorithm        = AbsorptionRGBA
        shader           = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
        prerender        = VolumePrerender()
        postrender       = () -> begin
            glDisable(GL_CULL_FACE)
        end
    end
end

function _default(main::IndirectArray{T}, s::Style, data::Dict) where T <: RGBA
    @gen_defaults! data begin
        volumedata       = main.index => Texture
        hull::GLUVWMesh  = AABB{Float32}(Vec3f0(0), Vec3f0(1))
        model            = Mat4f0(I)
        modelinv         = const_lift(inv, model)

        color_map        = main.values => TextureBuffer
        color_norm       = nothing
        color            = nothing

        algorithm        = IndexedAbsorptionRGBA
        shader           = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
        prerender        = VolumePrerender()
        postrender       = () -> begin
            glDisable(GL_CULL_FACE)
        end
    end
end


text:
using SeaPearl

struct Item
    id      :: Int
    value   :: Int
    weight  :: Int
end

mutable struct Solution
    content     :: AbstractArray{Bool}
    value       :: Int
    weight      :: Int
    optimality  :: Bool
end

struct InputData
    items               :: AbstractArray{Union{Item, Nothing}}
    sortedItems         :: AbstractArray{Union{Item, Nothing}}
    numberOfItems       :: Int
    capacity            :: Int
end

include("IOmanager.jl")

function solve_knapsack(filename::String; benchmark=false)

    input = parseFile!(filename)

    permutation = sortperm(input.items; by=(x) -> x.value/x.weight, rev=true)

    n = input.numberOfItems

    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)

    ### Variable declaration ###
    x_s = SeaPearl.IntVar[]
    x_a = SeaPearl.IntVar[]
    for i in 1:n
        push!(x_s, SeaPearl.IntVar(0, input.capacity, "x_s[" * string(i) * "]", trailer))
        push!(x_a, SeaPearl.IntVar(0, 1, "x_a[" * string(i) * "]", trailer))
        SeaPearl.addVariable!(model, last(x_s))
        SeaPearl.addVariable!(model, last(x_a))
    end

    push!(x_s, SeaPearl.IntVar(0, input.capacity, "x_s[" * string(n+1) * "]", trailer))
    SeaPearl.addVariable!(model, last(x_s))

    ### Constraints ###
    # Initial state: x_s[1] = 0
    initial = SeaPearl.EqualConstant(x_s[1], 0, trailer)
    push!(model.constraints, initial)

    # Transition: x_s[i+1] = x_s[i] + w[i]*x_a[i]
    for i in 1:n
        w_x_a_i = SeaPearl.IntVarViewMul(x_a[i], input.items[permutation[i]].weight, "w["*string(i)*"]*x_a["*string(i)*"]")
        minusX_s = SeaPearl.IntVarViewOpposite(x_s[i+1], "-x_s["*string(i+1)*"]")
        SeaPearl.addVariable!(model, w_x_a_i)
        SeaPearl.addVariable!(model, minusX_s)
        vars = SeaPearl.AbstractIntVar[w_x_a_i, minusX_s, x_s[i]]
        transition = SeaPearl.SumToZero(vars, trailer)
        push!(model.constraints, transition)
    end

    ### Objective ### minimize: -sum(v[i]*x_a[i])
    vars = SeaPearl.AbstractIntVar[]
    maxValue = 0
    for i in 1:n
        vx_a_i = SeaPearl.IntVarViewMul(x_a[i], input.items[permutation[i]].value, "v["*string(i)*"]*x_a["*string(i)*"]")
        push!(vars, vx_a_i)
        maxValue += input.items[permutation[i]].value
    end
    y = SeaPearl.IntVar(-maxValue, 0, "y", trailer)
    SeaPearl.addVariable!(model, y)
    push!(vars, y)
    objective = SeaPearl.SumToZero(vars, trailer)
    push!(model.constraints, objective)
    SeaPearl.addObjective!(model, y)

    status = SeaPearl.solve!(model; variableHeuristic=selectVariable)

    if !benchmark
        print(status)
        for oneSolution in model.statistics.solutions
            output = solutionFromSeaPearl(oneSolution, input, permutation)
            printSolution(output)
        end
    end
    return status
end

function solve_knapsack_without_dp(filename::String; benchmark=false)
    input = parseFile!(filename)

    permutation = sortperm(input.items; by=(x) -> x.value/x.weight, rev=true)

    n = input.numberOfItems

    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)

    ### Variable declaration ###
    x = SeaPearl.IntVar[]
    for i in 1:n
        push!(x, SeaPearl.IntVar(0, 1, "x[" * string(i) * "]", trailer))
        SeaPearl.addVariable!(model, last(x))
    end

    ### Constraints ###

    # Creating the totalWeight variable
    varsWeight = SeaPearl.AbstractIntVar[]
    maxWeight = 0
    for i in 1:n
        wx_i = SeaPearl.IntVarViewMul(x[i], input.items[permutation[i]].weight, "w["*string(i)*"]*x["*string(i)*"]")
        push!(varsWeight, wx_i)
        maxWeight += input.items[permutation[i]].weight
    end
    totalWeight = SeaPearl.IntVar(0, maxWeight, "totalWeight", trailer)
    minusTotalWeight = SeaPearl.IntVarViewOpposite(totalWeight, "-totalWeight")
    SeaPearl.addVariable!(model, totalWeight)
    SeaPearl.addVariable!(model, minusTotalWeight)
    push!(varsWeight, minusTotalWeight)
    weightEquality = SeaPearl.SumToZero(varsWeight, trailer)
    push!(model.constraints, weightEquality)

    # Making sure it is below the capacity
    weightConstraint = SeaPearl.LessOrEqualConstant(totalWeight, input.capacity, trailer)
    push!(model.constraints, weightConstraint)

    ### Objective ### minimize: -sum(v[i]*x_a[i])

    # Creating the sum
    varsValue = SeaPearl.AbstractIntVar[]
    maxValue = 0
    for i in 1:n
        vx_i = SeaPearl.IntVarViewMul(x[i], input.items[permutation[i]].value, "v["*string(i)*"]*x["*string(i)*"]")
        push!(varsValue, vx_i)
        maxValue += input.items[permutation[i]].value
    end
    totalValue = SeaPearl.IntVar(-maxValue, 0, "totalValue", trailer)
    SeaPearl.addVariable!(model, totalValue)
    push!(varsValue, totalValue)
    valueEquality = SeaPearl.SumToZero(varsValue, trailer)
    push!(model.constraints, valueEquality)

    # Setting it as the objective
    SeaPearl.addObjective!(model, totalValue)

    status = SeaPearl.solve!(model; variableHeuristic=selectVariableWithoutDP)

    if !benchmark
        print(status)
        for oneSolution in model.statistics.solutions
            output = solutionFromSeaPearlWithoutDP(oneSolution, input, permutation)
            printSolution(output)
        end
    end
    return status
end

function selectVariable(model::SeaPearl.CPModel)
    i = 1
    while SeaPearl.isbound(model.variables["x_a[" * string(i) * "]"])
        i += 1
    end
    return model.variables["x_a[" * string(i) * "]"]
end

function selectVariableWithoutDP(model::SeaPearl.CPModel)
    i = 1
    while SeaPearl.isbound(model.variables["x[" * string(i) * "]"])
        i += 1
    end
    return model.variables["x[" * string(i) * "]"]
end

function solutionFromSeaPearl(SeaPearlSol::SeaPearl.Solution, input::InputData, permutation::Array{Int})
    taken = falses(input.numberOfItems)
    value = 0
    weight = 0
    for i in 1:input.numberOfItems
        if haskey(SeaPearlSol, "x_a[" * string(i) * "]")
            taken[permutation[i]] = convert(Bool, SeaPearlSol["x_a[" * string(i) * "]"])
            if taken[permutation[i]]
                value += input.items[permutation[i]].value
                weight += input.items[permutation[i]].weight
            end
        end
    end
    return Solution(taken, value, weight, false)
end

function solutionFromSeaPearlWithoutDP(SeaPearlSol::SeaPearl.Solution, input::InputData, permutation::Array{Int})
    taken = falses(input.numberOfItems)
    value = 0
    weight = 0
    for i in 1:input.numberOfItems
        if haskey(SeaPearlSol, "x[" * string(i) * "]")
            taken[permutation[i]] = convert(Bool, SeaPearlSol["x[" * string(i) * "]"])
            if taken[permutation[i]]
                value += input.items[permutation[i]].value
                weight += input.items[permutation[i]].weight
            end
        end
    end
    return Solution(taken, value, weight, false)
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
#!/usr/local/bin/julia
# author: zhwang
# date: 2022-03-01
# description: 将模板代码转化为提交文件，模板放置路径为 `仓库/src/reviewed`，导出路径为 `仓库/submit/problems` 和 `仓库/submit/test`

### 检查和准备阶段 ###
cd(Base.source_dir()[1:end-5])
current_files = readdir()
!all(file in current_files for file in ["src", "tools"]) && throw("脚本文件位置有误")

reviewed_path = "src/reviewed/"
reviewed_files = readdir(reviewed_path)
isempty(reviewed_files) && throw("请将文件放在 src/reviewed 目录下")
mkpath.(["submit", "submit/problems", "submit/test"]) # 创建目录

function get_data(content::String)
    ## 获取测试代码
    test_reg = r"## @lc test=start\n([\s\S]*)\n## @lc test=end"
    test = match(test_reg, content)
    isnothing(test) && throw("未识别到测试信息")
    ## 获取正文代码
    body_reg = r"[\s\S]*## @lc code=end"
    body = match(body_reg, content)
    isnothing(body) && throw("未识别到正文内容")
    return body.match, test.captures[1]
end

### 生成文件 ###
problem_path = "submit/problems/"
test_path = "submit/test/"
for file in reviewed_files
    println("正在处理文件 $file")
    content = read(open(reviewed_path * file, "r"), String)
    body, test = get_data(content)
    open(problem_path * file, "w") do io
        write(io, body)
    end
    open(test_path * file, "w") do io
        write(io, test)
    end
end
println("文件处理完毕，一共写入 $(length(reviewed_files)) 个文件")


text:
function get_counts(args...)
    key = ENV["USDA_QUICK_SURVEY_KEY"]
    header = string(usda_url, "/api/get_counts/?key=", key)

    query = ""

    for i in args
        arg = string("&", i)
        query *= arg
    end

    request("GET", string(header, query))
end


text:
module JuliaExampleTest

using ForwardDiff

greet() = print("hello word!")

include("subfunction.jl")

export subfunction, derivative_of_my_f

end


text:
# Helper data for color parsing
include("names_data.jl")

# Color Parsing
# -------------

const col_pat_hex1 = r"(#|0x)([[:xdigit:]])([[:xdigit:]])([[:xdigit:]])"
const col_pat_hex2 = r"(#|0x)([[:xdigit:]]{2})([[:xdigit:]]{2})([[:xdigit:]]{2})"
const col_pat_rgb  = r"rgb\((\d+%?),(\d+%?),(\d+%?)\)"
const col_pat_hsl  = r"hsl\((\d+%?),(\d+%?),(\d+%?)\)"
const col_pat_rgba = r"rgba\((\d+%?),(\d+%?),(\d+%?),(\d+(?:\.\d*)?%?)\)"
const col_pat_hsla = r"hsla\((\d+%?),(\d+%?),(\d+%?),(\d+(?:\.\d*)?%?)\)"

# Parse a number used in the "rgb()" or "hsl()" color.
function parse_rgb(num::AbstractString)
    if num[end] == '%'
        return clamp(parse(Int, num[1:end-1], 10) / 100, 0, 1)
    else
        return clamp(parse(Int, num, 10) / 255, 0, 1)
    end
end

function parse_hsl_hue(num::AbstractString)
    if num[end] == '%'
        error("hue cannot end in %")
    else
        return parse(Int, num, 10)
    end
end

function parse_hsl_sl(num::AbstractString)
    if num[end] != '%'
        error("saturation and lightness must end in %")
    else
        return parse(Int, num[1:end-1], 10) / 100
    end
end

# Parse a number used in the alpha field of "rgba()" and "hsla()".
function parse_alpha_num(num::AbstractString)
    if num[end] == '%'
        return parse(Int, num[1:end-1]) / 100
    else
        return parse(Float32, num)
    end
end

function _parse_colorant(desc::AbstractString)
    desc_ = replace(desc, " ", "")
    mat = match(col_pat_hex2, desc_)
    if mat != nothing
        return RGB{N0f8}(parse(Int, mat.captures[2], 16) / 255,
                       parse(Int, mat.captures[3], 16) / 255,
                       parse(Int, mat.captures[4], 16) / 255)
    end

    mat = match(col_pat_hex1, desc_)
    if mat != nothing
        return RGB{N0f8}((16 * parse(Int, mat.captures[2], 16)) / 255,
                       (16 * parse(Int, mat.captures[3], 16)) / 255,
                       (16 * parse(Int, mat.captures[4], 16)) / 255)
    end

    mat = match(col_pat_rgb, desc_)
    if mat != nothing
        return RGB{N0f8}(parse_rgb(mat.captures[1]),
                       parse_rgb(mat.captures[2]),
                       parse_rgb(mat.captures[3]))
    end

    mat = match(col_pat_hsl, desc_)
    if mat != nothing
        return HSL{ColorTypes.eltype_default(HSL)}(parse_hsl_hue(mat.captures[1]),
                                        parse_hsl_sl(mat.captures[2]),
                                        parse_hsl_sl(mat.captures[3]))
    end

    mat = match(col_pat_rgba, desc_)
    if mat != nothing
        return RGBA{N0f8}(parse_rgb(mat.captures[1]),
                        parse_rgb(mat.captures[2]),
                        parse_rgb(mat.captures[3]),
                        parse_alpha_num(mat.captures[4]))
    end

    mat = match(col_pat_hsla, desc_)
    if mat != nothing
        return HSLA{ColorTypes.eltype_default(HSLA)}(parse_hsl_hue(mat.captures[1]),
                                          parse_hsl_sl(mat.captures[2]),
                                          parse_hsl_sl(mat.captures[3]),
                                          parse_alpha_num(mat.captures[4]))
    end

    desc_ = lowercase(desc_)

    if desc_ == "transparent"
        return RGBA{N0f8}(0,0,0,0)
    end

    if !haskey(color_names, desc_)
        error("Unknown color: ", desc)
    end

    c = color_names[desc_]
    return RGB{N0f8}(c[1] / 255, c[2] / 255, c[3] / 255)
end

# note: these exist to enable proper dispatch, since super(Colorant) == Any
_parse_colorant{C<:Colorant,SUP<:Any}(::Type{C}, ::Type{SUP}, desc::AbstractString) = _parse_colorant(desc)
_parse_colorant{C<:Colorant,SUP<:Colorant}(::Type{C}, ::Type{SUP}, desc::AbstractString) = convert(C, _parse_colorant(desc))::C

"""
    parse(Colorant, desc)

Parse a color description.

This parses subset of HTML/CSS color specifications. In particular, everything
is supported but: "currentColor".

It does support named colors (though it uses X11 named colors, which are
slightly different than W3C named colors in some cases), "rgb()", "hsl()",
"#RGB", and "#RRGGBB' syntax.

Args:
- `Colorant`: literal "Colorant" will parse according to the `desc`
string (usually returning an `RGB`); any more specific choice will
return a color of the specified type.
- `desc`: A color name or description.

Returns:
  An `RGB{N0f8}` color, unless:
    - "hsl(h,s,l)" was used, in which case an `HSL` color;
    - "rgba(r,g,b,a)" was used, in which case an `RGBA` color;
    - "hsla(h,s,l,a)" was used, in which case an `HSLA` color;
    - a specific `Colorant` type was specified in the first argument
"""
Base.parse{C<:Colorant}(::Type{C}, desc::AbstractString) = _parse_colorant(C, supertype(C), desc)
Base.parse{C<:Colorant}(::Type{C}, desc::Symbol) = parse(C, string(desc))
Base.parse{C<:Colorant}(::Type{C}, c::Colorant) = c

macro colorant_str(ex)
    isa(ex, AbstractString) || error("colorant requires literal strings")
    col = parse(Colorant, ex)
    :($col)
end

@noinline function ColorTypes.color(str::AbstractString)
    Base.depwarn("color(\"$str\") is deprecated, use colorant\"$str\" or parse(Colorant, \"$str\")", :color)
    parse(Colorant, str)
end


text: struct NoConfig <: AbstractConfig end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
function process_index(request::HTTP.Request, state::HandlerState)
    get_cache(state).need_recache && rebuild_cache!(state)
    return HTTP.Response(
        200,
        ["Content-Type" => "text/html"],
        body = state.cache.index_string
        )
end


text:
using Gridap
using Gridap.Io
using GridapGmsh

model = GmshDiscreteModel("elasticFlagFine.msh")

writevtk(model,"elasticFlagFine")

fn = "elasticFlagFine.json"
to_json_file(model,fn)


text:
abstract type LossUnsupervised <: Loss end

##############################################################################
# todo: loss should be a function of x and xhat vectors

#########################################
# Frobenius norm loss
#########################################

struct FrobeniusLoss <: LossUnsupervised
end

function eval(L::FrobeniusLoss, C, X, Y)
    return vecnorm(X*Y' - C)^2
end

# TODO: Find some way of more nicely enumerating these things
function deriv(L::FrobeniusLoss, C, X, Y, which="X")
    if which=="X"
        return 2*(X*Y' - C)*Y
    end
    if which == "Y"
        return 2*(Y*X' - C')*X
    end
    throw("which is defined as $(which), not one of X or Y")
end


text:
export track

"""
Use: track(   username::String,
                delay::Number,
                iterations::Int,
                command::String,
                logfile::String)

Tracks CPU stats on local processor using the `top` command.

Returns the path to the run's log file.
"""
function track( username::String,
                delay::Number,
                iterations::Int,
                command::String,
                logfile::String,
                verbose::Bool)

    # Generate commands from input
    top_cmd = `top -b -u $username -d $delay -n $iterations`
    grep_cmd = `grep "julia"`

    verbose && println("CPUmonitor: Begin Tracking...")

    # top -b -u $username -d $delay -n $iterations | grep "julia" > logfile
    run(pipeline(pipeline(top_cmd, grep_cmd), stdout = logfile))

    verbose && println("CPUmonitor: End Tracking...")

    return logfile
end


text:
using Test
using Printf
using Vizinanigans

function run_script(filepath)
    try
        include(filepath)
    catch err
        @warn "Error while testing script: " * sprint(showerror, err)

        # Print the content of the file to the test log, with line numbers, for debugging
        file_content = read(filepath, String)
        delineated_file_content = split(file_content, '\n')
        for (number, line) in enumerate(delineated_file_content)
            @printf("% 3d %s\n", number, line)
        end

        return false
    end
    return true
end

@testset "Vizinanigans" begin
    @testset "2D" begin
        @test run_script(joinpath(@__DIR__, "..", "examples", "visualize_2D.jl"))
    end

    @testset "3D" begin
        @test run_script(joinpath(@__DIR__, "..", "examples", "visualize_3D.jl"))
    end

    @testset "Volume slices" begin
        @test run_script(joinpath(@__DIR__, "..", "examples", "volumeslice.jl"))
    end
end


text:
using HelFEM
using Test
using LinearAlgebra: I
using SparseArrays: SparseVector, SparseMatrixCSC

@testset "HelFEM.jl" begin
    @testset "helfem" begin
        using HelFEM: helfem
        b = helfem.basis(10, 20, 4, 40, 2, 2.0, 0)
        @test helfem.get_nquad(b) == 50
        @test helfem.get_poly_id(b) == 4
        @test helfem.get_poly_order(b) == 10
        m = helfem.radial_integral(b, b, 0, false, false)
        @test size(m, 1) === 179
        @test size(m, 2) === 179
        @test size(m, 3) === 1
        @test size(m) == (179, 179)
        @test_throws ArgumentError size(m, 0)
        @test_throws ArgumentError size(m, -100)

        # Wrapping of arma::vec and arma::mat
        @test collect(helfem.ArmaVector([1,2,3])) == [1,2,3]
        @test collect(helfem.ArmaMatrix([1 2 3; 4 5 6])) == [1 2 3; 4 5 6]
        # Make sure we accept general vectors and matrices
        @test collect(helfem.ArmaVector(SparseVector([1,0,0,2]))) == [1,0,0,2]
        @test collect(helfem.ArmaMatrix(SparseMatrixCSC([1 0 0; 0 5 0]))) == [1 0 0; 0 5 0]
    end

    b1 = HelFEM.RadialBasis(16, 10)
    b2 = HelFEM.RadialBasis(6, 170; rmax = 50.0)

    # Make sure that our basis function evaluation functions give the same results.
    b1(HelFEM.quadraturepoints(b1)) ≈ HelFEM.basisvalues(b1)
    b2(HelFEM.quadraturepoints(b2)) ≈ HelFEM.basisvalues(b2)

    # Each element has the same number of basis functions as the order of the polynomials
    # (nnodes), but then we merge the ones on the element boundaries together and also
    # remove one function from each edge to satisfy boundary conditions
    nbf(nnodes, nelems) = nnodes * nelems - (nelems - 1) - 2
    @test length(b1) == nbf(16, 10)
    @test length(b2) == nbf(6, 170)

    # Introspection for element boundaries
    @test_throws Exception HelFEM.elementrange(b1, 0)
    @test_throws Exception HelFEM.elementrange(b1, 11)
    @test HelFEM.elementrange(b1, 1)[1] == 0.0
    @test HelFEM.elementrange(b2, 1)[1] == 0.0
    @test HelFEM.elementrange(b1, 10)[2] == 40.0
    @test HelFEM.elementrange(b2, 170)[2] == 50.0

    let bs = HelFEM.boundaries(b1)
        @test length(bs) == 11
        @test bs[1] == 0.0
        @test bs[end] == 40.0
    end
    let bs = HelFEM.boundaries(b2)
        @test length(bs) == 171
        @test bs[1] == 0.0
        @test bs[end] == 50.0
    end

    # Check properties of control points. Specifically, we're checking that for LIPs, the
    # basis function values at control points are equal to 1/r and zero otherwise.
    function check_control_point_matrix(b)
        cs = HelFEM.controlpoints(b)
        @test b(cs) .* cs ≈ [
            # The ones are on the first lower diagonal (i.e. below the main diagonal),
            # hence the +1:
            i == j + 1 ? 1 : 0
            for i = 1:length(cs), j = 1:length(b)
        ]
    end
    check_control_point_matrix(b1)
    check_control_point_matrix(b2)

    # Scaling between polynomial coordinates to radial coordinates
    @test_throws Exception HelFEM.scale_to_element(b1, 0, 0)
    @test_throws Exception HelFEM.scale_to_element(b1, 11, 0)
    let rs = HelFEM.scale_to_element(b1, 1, -1)
        @test isa(rs, Number)
        @test rs == 0.0
    end
    let rs = HelFEM.scale_to_element(b1, 10, 1)
        @test isa(rs, Number)
        @test rs == 40.0
    end
    let rs = HelFEM.scale_to_element(b1, 1, [-1, 0, 1])
        @test length(rs) == 3
        @test rs[1] == 0.0
    end
    let rs = HelFEM.scale_to_element(b2, 170, [1, -0.5, 0, 0.5, 1])
        @test length(rs) == 5
        @test rs[1] == 50.0
        @test rs[5] == 50.0
        @test rs[2]+rs[4] ≈ 2*rs[3]
    end

    S1 = HelFEM.overlap(b1)
    @test size(S1) == (length(b1), length(b1))
    S2, S2invh = HelFEM.overlap(b2, invh=true)
    @test size(S2) == (length(b2), length(b2))
    @test size(S2invh) == (length(b2), length(b2))

    @test S2invh' * S2 * S2invh ≈ I
    let S2invh_julia = sqrt(inv(S2))
        @test isreal(S2invh_julia)
        S2invh_julia = real(S2invh_julia)
        @test S2invh_julia' * S2 * S2invh_julia ≈ I
    end

    S1_rint = HelFEM.radial_integral(b1, 0)
    @test S1_rint ≈ S1
    S1_rint = HelFEM.radial_integral(b1, 0, b1)
    @test S1_rint ≈ S1

    S12_rint = HelFEM.radial_integral(b1, 0, b2)
    @test size(S12_rint) == (length(b1), length(b2))

    @test HelFEM.add_boundary!(b1, 1.23) === nothing
    @test length(HelFEM.boundaries(b1)) == 12

    @testset "FEMBasis" begin
        nquad = length(HelFEM.quadraturepoints(b1))
        grid = HelFEM.boundaries(b1)
        pb = HelFEM.PolynomialBasis(:lip, b1.nnodes)
        FB = HelFEM.FEMBasis(pb, grid, nquad=nquad)

        @test HelFEM.boundaries(FB) == grid
        @test HelFEM.elementrange(FB, 1) == (grid[1], grid[2])
        @test HelFEM.elementrange(FB, HelFEM.nelements(FB)) == (grid[end-1], grid[end])
        @test_throws DomainError HelFEM.elementrange(FB, 0)
        @test_throws DomainError HelFEM.elementrange(FB, -1)
        @test_throws DomainError HelFEM.elementrange(FB, HelFEM.nelements(FB)+1)
        @test_throws DomainError HelFEM.elementrange(FB, HelFEM.nelements(FB)+1000)

        @test length(HelFEM.controlpoints(FB)) == (b1.nnodes - 1)*HelFEM.nelements(FB) + 1

        let rs = HelFEM.quadraturepoints(b1)
            @test FB(rs) ./ rs ≈ HelFEM.basisvalues(b1)
        end
        @test HelFEM.radial_integral(FB, r -> 1.0) ≈ HelFEM.overlap(b1)
        @test HelFEM.radial_integral(FB, r -> 1/r) ≈ HelFEM.radial_integral(b1, -1)
        @test HelFEM.radial_integral(FB, r -> r, rderivative=true) ≈ HelFEM.radial_integral(b1, 1; rderivative=true)
        @test HelFEM.radial_integral(FB, r -> r, lderivative=true) ≈ HelFEM.radial_integral(b1, 1; lderivative=true)
        @test HelFEM.radial_integral(FB, r -> r, lderivative=true, rderivative=true) ≈ HelFEM.radial_integral(b1, 1; lderivative=true, rderivative=true)
    end
end


text:
abstract type AbstractReservesFormulation <: AbstractServiceFormulation end
struct RangeUpwardReserve <: AbstractReservesFormulation end
############################### Reserve Variables` #########################################
"""
This function add the variables for reserves to the model
"""
function activeservice_variables!(canonical::Canonical,
                                  service::SR,
                                  devices::Vector{<:PSY.Device}) where SR<:PSY.Reserve
    add_variable(canonical,
                 devices,
                 Symbol("R$(PSY.get_name(service))_$SR"),
                 false;
                 ub_value = d -> d.tech.activepowerlimits.max,
                 lb_value = d -> 0 )
    return
end

################################## Reserve Requirement Constraint ##########################
# This function can be generalized later for any constraint of type Sum(req_var) >= requirement,
# it will only need to be specific to the names and get forecast string.
function service_requirement_constraint!(canonical::Canonical,
                                         service::SR) where {SR<:PSY.Reserve}
    time_steps = model_time_steps(canonical)
    parameters = model_has_parameters(canonical)
    forecast = model_uses_forecasts(canonical)
    initial_time = model_initial_time(canonical)
    reserve_variable = get_variable(canonical, Symbol("R$(PSY.get_name(service))_$SR"))
    constraint_name = Symbol(PSY.get_name(service), "_requirement_$SR")
    constraint = add_cons_container!(canonical, constraint_name, time_steps)
    requirement = PSY.get_requirement(service)
    if forecast
        ts_vector = TS.values(PSY.get_data(PSY.get_forecast(PSY.Deterministic,
                                                            service,
                                                            initial_time,
                                                            "get_requirement")))
    else
        ts_vector = ones(time_steps[end])
    end
    if parameters
        param = include_parameters(canonical, ts_vector,
                                   UpdateRef{SR}("get_requirement"), time_steps)
        for t in time_steps
            constraint[t] = JuMP.@constraint(canonical.JuMPmodel,
                                         sum(reserve_variable[:,t]) >= param[t]*requirement)
        end
    else
        for t in time_steps
            constraint[t] = JuMP.@constraint(canonical.JuMPmodel,
                                    sum(reserve_variable[:,t]) >= ts_vector[t]*requirement)
        end
    end
    return
end

# This function can also be generalized.
function add_to_service_expression!(canonical::Canonical,
                                    model::ServiceModel{SR, RangeUpwardReserve},
                                    service::SR,
                                    expression_list::Vector{Symbol}) where {SR<:PSY.Reserve}
    # Container
    time_steps = model_time_steps(canonical)
    devices = PSY.get_contributingdevices(service)
    var_type = JuMP.variable_type(canonical.JuMPmodel)
    expression = :upward_reserve
    expression ∉ expression_list && push!(expression, expression_list)
    expression_dict = get_expression!(canonical,
                                      expression,
                                      Dict{ServiceExpressionKey, Array{GAE{var_type}}}())
    reserve_variable = get_variable(canonical, Symbol("R$(PSY.get_name(service))_$SR"))
    #fill container
    for d in devices
        T = typeof(d)
        name = PSY.get_name(d)
        expressions = get!(expression_dict,
                           ServiceExpressionKey(name, T),
                           get_variable(canonical, Symbol("P_$(T)"))[name, :].data)
        for t in time_steps
            expressions[t] = JuMP.add_to_expression!(expressions[t], reserve_variable[name, t])
        end
    end
    return
end

function add_device_constraints!(canonical::Canonical,
                                 sys::PSY.System,
                                 expression_list::Vector{Symbol})
    for service_expression in expression_list
        expression_dictionary = get_expression(canonical, service_expression)
        for (k, v) in expression_dictionary
            device = get_component(k.device_type, sys, k.name)
            service_constraints(device, v, )
        end
    end
    return
end