text:
######################################################################
# data structures
######################################################################

"""
Keep track of spell durations using an interval type.
"""
typealias DateInterval ClosedInterval{Date}

"""
Define a constructor for `T`, which accepts DateInterval fields as two
arguments.
"""
macro constructor_DateInterval(t::Symbol)
    T = eval(t)
    @assert isa(T, DataType) "$t needs to be a type."
    outer = []
    inner = []
    for name in fieldnames(T)
        ft = fieldtype(T, name)
        if ft == DateInterval
            arg_start = gensym(name)
            arg_end = gensym(name)
            push!(outer, :($arg_start::Date), :($arg_end::Date))
            push!(inner, :(DateInterval($arg_start, $arg_end)))
        else
            push!(outer, :($name::$ft))
            push!(inner, name)
        end
    end
    @eval $(t)($(outer...))=$(t)($(inner...))
end

"""
Information on (spells of) personal characteristics (nationality and
academic degree may change over time).

- based on HV data
- Smooting: No; Overlap: No
"""
@auto_hash_equals immutable PersonSpell
    birth_year::Int
    gender::Gender
    nationality::Int
    # NOTE ignoring education::String
    death_date::Nullable{Date}
    interval::DateInterval
    # NOTE ignoring timestamp::Date
end

@constructor_DateInterval PersonSpell

"""
Parser for person spells (`hv_pn`). Currently ignore education
(unreliable), infer death date availability from whether death date is
missing, ignore timestamp.
"""
const PersonSpellParser =
    LineParser(PersonSpell,
               1 => Int,        # 1 birth year
               2 => Gender,     # 2 gender
               3 => Int,        # 3 nationality
               # 4 education levels: ignored
               # 5 death date availability: ignored
               6 => MaybeEmpty{AMDB_Date}, # 6 death date
               7 => AMDB_Date,             # 7 spell start
               8 => AMDB_Date)             # 8 spell end

"""
A spell from the AMP datasets, eg `hv_epi_uni_roh`.
"""
@auto_hash_equals immutable AMPSpell
    employer::Int
    status::AMP.Spell
    interval::DateInterval
end

@constructor_DateInterval AMPSpell

"""
Parser for AMP spells (eg `hv_epi_uni_roh`). Ignore timestamp, and
last field (imputation indicator?).
"""
const AMPSpellParser =
    LineParser(AMPSpell,
               1 => Int,            # employed ID
               2 => AMP.spell_dict, # spell status
               3 => AMDB_Date,      # spell start
               4 => AMDB_Date)      # spell end
"""
A spell from the `vmz_v2_int` dataset.
"""
@auto_hash_equals immutable VMZSpell
    status::VMZ.Spell
    discontinuation::VMZ.Einstellgrund
    interval::DateInterval
end

@constructor_DateInterval VMZSpell

"""
Parser for VMZ spells (`vmz_v2_int`). Ignore case id (first field).
"""
const VMZSpellParser =
    LineParser(VMZSpell,
               3 => VMZ.spell_dict,         # spell status
               4 => VMZ.einstellgrund_dict, # einstellgrund
               1 => AMDB_Date,              # spell start
               2 => AMDB_Date)              # spell end

"""
Data matched for an individual, from various data files.
"""
@auto_hash_equals immutable IndividualData
    person_spells::Vector{PersonSpell}
    AMP_spells::Vector{AMPSpell}
    VMZ_spells::Vector{VMZSpell}
end

IndividualData() = IndividualData(Vector{PersonSpell}(),
                                  Vector{AMPSpell}(),
                                  Vector{VMZSpell}())


text:
"""
    TaylorVector{N,T}

A data structure representing a vector of Taylor series with `N` terms.
Each element is an `NTuple{N,T}`.

    TaylorVector(T, N::Integer, n::Integer)

Create a vector of `n` `NTuple{N,T}`s.
"""
struct TaylorVector{N,T} <: AbstractVector{NTuple{N,T}}
    data::LinearAlgebra.Transpose{T,Matrix{T}}
    views::NTuple{N,SubArray{T,1,Matrix{T},Tuple{Int,UnitRange{Int}},true}}
end
TaylorVector{N}(TV::TaylorVector{M,T}) where {N,M,T} =
    TaylorVector{N,T}(TV.data, TV.views[1:N])
function TaylorVector{N}(data::Matrix) where {N}
    views = tuple((view(data, i, 1:size(data, 2)) for i = 1:N)...)
    TaylorVector(LinearAlgebra.transpose(data), views)
end
function TaylorVector{N}(T, n::Integer) where {N}
    TaylorVector{N}(zeros(T, N, n))
end

function Base.show(io::IO, ::MIME"text/plain", X::TaylorVector)
    summary(io, X)
    isempty(X) && return
    println(io, ":")
    Base.print_array(io, map(i -> X[i], 1:length(X)))
end

Base.length(TV::TaylorVector) = size(TV.data, 1)
Base.size(TV::TaylorVector) = (length(TV),)
Base.eltype(::Type{TaylorVector{N,T}}) where {N,T} = NTuple{N,T}
Base.IndexStyle(::TaylorVector) = Base.IndexLinear()

"""
    vectors(TV::TaylorVec{N})

Return the Taylor series as `N` seperate vectors.
"""
vectors(TV::TaylorVector) = TV.views

@generated function Base.getindex(TV::TaylorVector{N}, i::Integer) where {N}
    quote
        Base.@_propagate_inbounds_meta
        x = TV.data
        $(Expr(:tuple, (:(x[i, $k]) for k = 1:N)...))
    end
end
Base.getindex(TV::TaylorVector, i::Integer, j::Integer) = getindex(TV.data, i, j)

function Base.setindex!(TV::TaylorVector{N,T}, x, i::Integer) where {N,T}
    setindex!(TV, convert(NTuple{N,T}, x), i)
end
function Base.setindex!(TV::TaylorVector{1}, x::Number, i::Integer)
    TV.data[i] = x
    x
end
@generated function Base.setindex!(
    TV::TaylorVector{N,T},
    x::NTuple{N,T},
    i::Integer,
) where {N,T}
    quote
        Base.@_propagate_inbounds_meta
        d = TV.data
        $(Expr(:tuple, (:(d[i, $k]) for k = 1:N)...)) = x
        x
    end
end
Base.setindex!(TV::TaylorVector, x, i::Integer, j::Integer) = setindex!(TV.data, x, i, j)

#####################
## CompiledSystem  ##
#####################
const TSYSTEM_TABLE = Dict{UInt,Vector{System}}()

struct CompiledSystem{HI}
    nexpressions::Int
    nvariables::Int
    nparameters::Int
    system::System
end

function CompiledSystem(F::System)
    n = length(F)
    nvars = nvariables(F)
    nparams = nparameters(F)

    # We substitute all variable and parameter names away such that two systems
    # compile to the same code independently of the names
    @var _x_[1:nvars] _p_[1:nparams]
    D = ExpressionMap()
    for (x, y) in zip(variables(F), _x_)
        D[x] = y
    end
    for (x, y) in zip(parameters(F), _p_)
        D[x] = y
    end
    cleard_exprs = subs(expressions(F), D)
    sys = System(cleard_exprs, _x_, _p_)
    h = hash(cleard_exprs)
    k = 0
    if haskey(TSYSTEM_TABLE, h)
        # check that it is identical
        for (i, vi) in enumerate(TSYSTEM_TABLE[h])
            if vi == sys
                k = i
                break
            end
        end
        if k == 0
            push!(TSYSTEM_TABLE[h], sys)
            k = length(TSYSTEM_TABLE[h])
        end
    else
        k = 1
        TSYSTEM_TABLE[h] = [sys]
    end
    return CompiledSystem{(h, k)}(n, nvars, nparams, F)
end

function Base.show(io::IO, TS::CompiledSystem)
    print(io, "Compiled: ")
    show(io, TS.system)
end

interpret(TS::CompiledSystem) = TS.system
interpret(::Type{CompiledSystem{HI}}) where {HI} = TSYSTEM_TABLE[first(HI)][last(HI)]

Base.size(CS::CompiledSystem) = (CS.nexpressions, CS.nvariables)
Base.size(CS::CompiledSystem, i::Integer) = size(CS)[i]
Base.length(CS::CompiledSystem) = CS.nexpressions
nparameters(CS::CompiledSystem) = CS.nparameters
nvariables(CS::CompiledSystem) = CS.nvariables

######################
## CompiledHomotopy ##
######################

const THOMOTOPY_TABLE = Dict{UInt,Vector{Homotopy}}()

struct CompiledHomotopy{HI}
    nexpressions::Int
    nvariables::Int
    nparameters::Int
    homotopy::Homotopy
end

function CompiledHomotopy(H::Homotopy)
    n = length(H)
    nvars = nvariables(H)
    nparams = nparameters(H)

    # We substitute all variable and parameter names away such that two homotopies
    # compile to the same code independently of the names
    @var _x_[1:nvars] _t_ _p_[1:nparams]
    D = ExpressionMap()
    for (x, y) in zip(variables(H), _x_)
        D[x] = y
    end
    for (x, y) in zip(parameters(H), _p_)
        D[x] = y
    end
    D[H.t] = _t_
    cleard_exprs = subs(expressions(H), D)
    homotopy = Homotopy(cleard_exprs, _x_, _t_, _p_)
    h = hash(cleard_exprs)

    k = 0
    if haskey(THOMOTOPY_TABLE, h)
        # check that it is identical
        for (i, vi) in enumerate(THOMOTOPY_TABLE[h])
            if vi == homotopy
                k = 1
                break
            end
        end
        if k == 0
            push!(THOMOTOPY_TABLE[h], homotopy)
            k = 1
        end
    else
        k = 1
        THOMOTOPY_TABLE[h] = [homotopy]
    end
    return CompiledHomotopy{(h, 1)}(n, nvars, nparams, H)
end

Base.size(CH::CompiledHomotopy) = (CH.nexpressions, CH.nvariables)
Base.size(CH::CompiledHomotopy, i::Integer) = size(CH)[i]
Base.length(CH::CompiledHomotopy) = CH.nexpressions
nparameters(CH::CompiledHomotopy) = CH.nparameters
nvariables(CH::CompiledHomotopy) = CH.nvariables

function Base.show(io::IO, TH::CompiledHomotopy)
    print(io, "Compiled: ")
    show(io, TH.homotopy)
end

interpret(CH::CompiledHomotopy) = CH.homotopy
interpret(::Type{CompiledHomotopy{HI}}) where {HI} = THOMOTOPY_TABLE[first(HI)][last(HI)]

"""
    compile(F::System; optimizations = true)

Compile the given system. Returns a `CompiledSystem`. If `optimizations == true` then
the given system is optimized for evaluation efficiency.

    compile(H::Homotopy)

Compile the given homotopy. Returns a `CompiledHomotopy`.
"""
compile(F::System; optimizations::Bool = true) =
    CompiledSystem(optimizations ? optimize(F) : F)
compile(H::Homotopy) = CompiledHomotopy(H)

#############
## CODEGEN ##
#############

boundscheck_var_map(F::System; kwargs...) =
    boundscheck_var_map(F.expressions, F.variables, F.parameters; kwargs...)
boundscheck_var_map(H::Homotopy; kwargs...) =
    boundscheck_var_map(H.expressions, H.variables, H.parameters, H.t; kwargs...)
function boundscheck_var_map(
    exprs,
    vars,
    params,
    t = nothing;
    taylor = false,
    jacobian::Bool = false,
)
    n = length(exprs)
    m = length(vars)
    l = length(params)
    var_map = Dict{Symbol,Union{Symbol,Expr}}()
    for i = 1:m
        if taylor
            var_map[Symbol(vars[i])] = :(x[$i, 1])
        else
            var_map[Symbol(vars[i])] = :(x[$i])
        end
    end
    for i = 1:l
        if taylor
            var_map[Symbol(params[i])] = :(p[$i, 1])
        else
            var_map[Symbol(params[i])] = :(p[$i])
        end
    end
    if t !== nothing
        var_map[Symbol(t)] = :(t)
    end

    checks = Expr[]
    if jacobian
        push!(checks, :(@boundscheck checkbounds(U, 1:$n, 1:$m)))
    else
        push!(checks, :(@boundscheck checkbounds(u, 1:$n)))
    end
    push!(checks, :(@boundscheck checkbounds(x, 1:$m)))
    push!(checks, :(@boundscheck p === nothing || checkbounds(p, 1:$l)))

    Expr(:block, checks...), var_map
end

function add_assignement!(D::Dict{Symbol,Vector{Expr}}, id::Symbol, e::Expr)
    if haskey(D, id)
        push!(D[id], e)
    else
        D[id] = [e]
    end
    D
end

function _evaluate!_impl(::Type{T}) where {T<:Union{CompiledSystem,CompiledHomotopy}}
    I = interpret(T)
    checks, var_map = boundscheck_var_map(I)
    slp = let
        list, ids = instruction_list(I.expressions)
        assignements = Dict{Symbol,Vector{Expr}}()
        for (i, id) in enumerate(ids)
            add_assignement!(assignements, id, :(u[$i] = $id))
        end
        to_expr(list, var_map, assignements)
    end

    quote
        $checks
        @inbounds $slp
        u
    end
end

function _jacobian!_impl(::Type{T}) where {T<:Union{CompiledSystem,CompiledHomotopy}}
    I = interpret(T)
    checks, var_map = boundscheck_var_map(I; jacobian = true)

    slp = let
        list, ids = instruction_list(I.expressions)
        vars = Symbol.(I.variables)
        params = Symbol.(I.parameters)
        dlist, J = diff(list, vars, ids)

        assignements = Dict{Symbol,Vector{Expr}}()

        U_constants = Expr[]
        for j = 1:size(J, 2), i = 1:size(J, 1)
            if J[i, j] isa Symbol
                if J[i, j] ∉ vars && J[i, j] ∉ params
                    add_assignement!(assignements, J[i, j], :(U[$i, $j] = $(J[i, j])))
                else
                    push!(U_constants, :(U[$i, $j] = $(var_map[J[i, j]])))
                end
            elseif J[i, j] !== nothing
                push!(U_constants, :(U[$i, $j] = $(J[i, j])))
            end
        end
        expr = to_expr(dlist, var_map, assignements)
        append!(expr.args, U_constants)
        expr
    end
    quote
        $checks
        U .= zero(eltype(x))
        @inbounds $slp
        U
    end
end

function _evaluate_and_jacobian!_impl(
    ::Type{T},
) where {T<:Union{CompiledSystem,CompiledHomotopy}}
    I = interpret(T)
    checks, var_map = boundscheck_var_map(I; jacobian = true)

    slp = let
        list, ids = instruction_list(I.expressions)
        vars = Symbol.(I.variables)
        params = Symbol.(I.parameters)
        dlist, J = diff(list, vars, ids)

        assignements = Dict{Symbol,Vector{Expr}}()
        for (i, id) in enumerate(ids)
            add_assignement!(assignements, id, :(u[$i] = $id))
        end

        U_constants = Expr[]
        for j = 1:size(J, 2), i = 1:size(J, 1)
            if J[i, j] isa Symbol
                if J[i, j] ∉ vars && J[i, j] ∉ params
                    add_assignement!(assignements, J[i, j], :(U[$i, $j] = $(J[i, j])))
                else
                    push!(U_constants, :(U[$i, $j] = $(var_map[J[i, j]])))
                end
            elseif J[i, j] !== nothing
                push!(U_constants, :(U[$i, $j] = $(J[i, j])))
            end
        end
        expr = to_expr(dlist, var_map, assignements)
        append!(expr.args, U_constants)
        expr
    end
    quote
        $checks
        U .= zero(eltype(x))
        @inbounds $slp
        nothing
    end
end

Base.@propagate_inbounds function set_row!(u::AbstractMatrix, t::Tuple{A}, i) where {A}
    u[i, 1] = first(t)
end
Base.@propagate_inbounds function set_row!(u::AbstractMatrix, t::Tuple{A,B}, i) where {A,B}
    a, b = t
    u[i, 1] = a
    u[i, 2] = b
end
Base.@propagate_inbounds function set_row!(
    u::AbstractMatrix,
    t::Tuple{A,B,C},
    i,
) where {A,B,C}
    a, b, c = t
    u[i, 1] = a
    u[i, 2] = b
    u[i, 3] = c
end
Base.@propagate_inbounds function set_row!(
    u::AbstractMatrix,
    t::Tuple{A,B,C,D},
    i,
) where {A,B,C,D}
    a, b, c, d = t
    u[i, 1] = a
    u[i, 2] = b
    u[i, 3] = c
    u[i, 4] = d
end
Base.@propagate_inbounds function set_row!(
    u::AbstractMatrix,
    t::Tuple{A,B,C,D,E},
    i,
) where {A,B,C,D,E}
    a, b, c, d, e = t
    u[i, 1] = a
    u[i, 2] = b
    u[i, 3] = c
    u[i, 4] = d
    u[i, 5] = e
end

function _functions_taylor!_impl(
    ::Type{T},
    K::Int;
    highest_order_only::Bool,
) where {T<:Union{CompiledSystem,CompiledHomotopy}}
    I = interpret(T)
    checks, var_map = boundscheck_var_map(I)
    list, ids = instruction_list(I.expressions)
    assignements = Dict{Symbol,Vector{Expr}}()
    for (i, id) in enumerate(ids)
        if highest_order_only
            add_assignement!(assignements, id, :(u[$i] = last($id)))
        else
            add_assignement!(assignements, id, :(set_row!(u, $id, $i)))
        end
    end

    block = Expr(:block)
    exprs = block.args
    for (id, (op, arg1, arg2)) in list.instructions
        a = get(var_map, arg1, arg1)
        if op == :^
            r::Int = arg2
            if r == 2
                push!(exprs, :($id = taylor(Val{:sqr}, Val{$K}, $a)))
            else
                push!(exprs, :($id = taylor(Val{:^}, Val{$K}, $a, $r)))
            end
        elseif arg2 !== nothing
            b = get(var_map, arg2, arg2)
            push!(exprs, :($id = taylor(Val{$(QuoteNode(op))}, Val{$K}, $a, $b)))
        else
            push!(exprs, :($id = taylor(Val{$(QuoteNode(op))}, Val{$K}, $a)))
        end
        if haskey(assignements, id)
            append!(exprs, assignements[id])
        end
    end
    # TODO: let block only if homotopy....
    if I isa Homotopy
        quote
            $checks
            let t = (t, one(t))
                @inbounds $block
            end
            u
        end
    else
        quote
            $checks
            @inbounds $block
            u
        end
    end
end

function _inline_taylor!_impl(
    T::Type{<:Union{CompiledHomotopy,CompiledSystem}},
    K,
    dx,
    dp;
    highest_order_only::Bool,
)
    H = interpret(T)
    # @show H
    checks, var_map = boundscheck_var_map(H; taylor = true)

    list, ids = instruction_list(H.expressions)

    vars = Symbol.(H.variables)
    params = Symbol.(H.parameters)
    # @show vars, dx, dp
    # @show params
    diff_map = DiffMap()
    for (i, v) in enumerate(vars)
        for k = 1:dx
            diff_map[v, k] = :(x[$i, $(k + 1)])
        end
    end

    for (i, v) in enumerate(params)
        for k = 1:dp
            diff_map[v, k] = :(p[$i, $(k + 1)])
        end
    end

    if H isa Homotopy
        diff_map[Symbol(H.t), 1] = 1
    end
    # @show diff_map
    dlist = univariate_diff!(list, K, diff_map)

    assignements = Dict{Symbol,Vector{Expr}}()
    u_constants = Expr[]

    if highest_order_only
        for (i, id) in enumerate(ids)
            k = K
            k_id = diff_map[id, k]
            if k_id isa Symbol
                if k_id ∉ vars && k_id ∉ params
                    add_assignement!(assignements, k_id, :(u[$i] = $k_id))
                else
                    push!(u_constants, :(u[$i] = $(var_map[k_id])))
                end
            elseif k_id isa Nothing
                push!(u_constants, :(u[$i] = zero(eltype(u))))
            else
                push!(u_constants, :(u[$i] = $k_id))
            end
        end
    else
        for (i, id) in enumerate(ids)
            add_assignement!(assignements, id, :(u[$i, 1] = $id))
            for k = 1:K
                k_id = diff_map[id, k]
                if k_id isa Symbol
                    if k_id ∉ vars && k_id ∉ params
                        add_assignement!(assignements, k_id, :(u[$i, $(k + 1)] = $k_id))
                    else
                        push!(u_constants, :(u[$i, $(k + 1)] = $(var_map[k_id])))
                    end
                elseif k_id isa Nothing
                    push!(u_constants, :(u[$i, $(k + 1)] = zero(u[$i, $(k + 1)])))
                else
                    push!(u_constants, :(u[$i, $(k + 1)] = $k_id))
                end
            end
        end
    end
    slp = to_expr(dlist, var_map, assignements)
    append!(slp.args, u_constants)

    quote
        $checks
        @inbounds $slp
        u
    end
end

function _taylor!_impl(T, K, D, DP; kwargs...)
    # Experiments show that for K < 2 (i.e. u and u̇) the fully inlined version has better
    # compilation times, for K ≥ 2 the function version starts to become
    # significantly faster
    if K ≤ 1
        _inline_taylor!_impl(T, K, D, DP; kwargs...)
    else
        _functions_taylor!_impl(T, K; kwargs...)
    end
end

################
## EVALUATION ##
################

# inplace (generated)
"""
    evaluate!(u, T::CompiledSystem, x, p = nothing)

Evaluate `T` for variables `x` and parameters `p` and store result in `u`.
"""
@generated function evaluate!(u, T::CompiledSystem, x, p = nothing)
    _evaluate!_impl(T)
end

"""
    evaluate!(u, T::CompiledHomotopy, x, t, p = nothing)

Evaluate `T` for variables `x`, `t` and parameters `p` and store result in `u`.
"""
@generated function evaluate!(u, T::CompiledHomotopy, x, t, p = nothing)
    _evaluate!_impl(T)
end

"""
    jacobian!(U, T::CompiledHomotopy, x, p = nothing)

Evaluate the Jacobian of `T` for variables `x`, `t` and parameters `p`
and store result in `u`.
"""
@generated function jacobian!(U, T::CompiledSystem, x, p = nothing)
    _jacobian!_impl(T)
end

"""
    jacobian!(U, T::CompiledHomotopy, x, t, p = nothing)

Evaluate the Jacobian of `T` for variables `x`, `t` and parameters `p` and
store result in `u`.
"""
@generated function jacobian!(U, T::CompiledHomotopy, x, t, p = nothing)
    _jacobian!_impl(T)
end

"""
    evaluate_and_jacobian!(u, U, T::CompiledHomotopy, x, p = nothing)

Evaluate `T` and its Jacobian for variables `x` and parameters `p` and
store result in `u`.
"""
@generated function evaluate_and_jacobian!(u, U, T::CompiledSystem, x, p = nothing)
    _evaluate_and_jacobian!_impl(T)
end

"""
    evaluate_and_jacobian!(u, U, T::CompiledHomotopy, x, t, p = nothing)

Evaluate `T` and its Jacobian for variables `x`, `t` and parameters `p` and
store result in `u`.
"""
@generated function evaluate_and_jacobian!(u, U, T::CompiledHomotopy, x, t, p = nothing)
    _evaluate_and_jacobian!_impl(T)
end

"""
    taylor!(
        u::AbstractMatrix,
        v::Val{M}
        F::CompiledSystem,
        x::TaylorVector{D},
        p::Union{Nothing,Vector,TaylorVector} = nothing,
    )
    taylor!(
        u::TaylorVector{M},
        v::Val{M}
        F::CompiledSystem,
        x::TaylorVector{D},
        p::Union{Nothing,Vector,TaylorVector} = nothing,
    )

Compute the Taylor series of ``u = F(x,p)`` given the Taylor series `x` and `p`.
"""
@generated function taylor!(
    u::AbstractMatrix,
    v::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::Nothing = nothing,
) where {M,D}
    _taylor!_impl(T, M, D - 1, -1; highest_order_only = false)
end
@generated function taylor!(
    u::AbstractMatrix,
    v::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::AbstractVector,
) where {M,D}
    _taylor!_impl(T, M, D - 1, 0; highest_order_only = false)
end
@generated function taylor!(
    u::AbstractMatrix,
    v::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::TaylorVector{DP},
) where {M,D,DP}
    _taylor!_impl(T, M, D - 1, DP - 1; highest_order_only = false)
end

function taylor!(
    u::TaylorVector{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p = nothing,
) where {M,D}
    taylor!(u.data, Val(M - 1), T, x, p)
    u
end

@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::Nothing = nothing,
) where {M,D}
    _taylor!_impl(T, M, D - 1, -1; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::AbstractVector,
) where {M,D}
    _taylor!_impl(T, M, D - 1, 0; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::TaylorVector{D},
    p::TaylorVector{DP},
) where {M,D,DP}
    _taylor!_impl(T, M, D - 1, DP - 1; highest_order_only = true)
end

@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::AbstractVector,
    p::Nothing = nothing,
) where {M}
    _taylor!_impl(T, M, 0, -1; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::AbstractVector,
    p::AbstractVector,
) where {M}
    _taylor!_impl(T, M, 0, 0; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledSystem,
    x::AbstractVector,
    p::TaylorVector{DP},
) where {M,DP}
    _taylor!_impl(T, M, 0, DP - 1; highest_order_only = true)
end

"""
    taylor!(
        u::AbstractVector,
        ::Val{M}
        H::Homotopy,
        x::TaylorVector{D},
        t::Number,
        p::Union{Nothing,Vector,TaylorVector} = nothing,
    )

Compute the `M`-th derivative of ``u=H(x,t,p)`` given the taylor series `x` and `p`.
"""
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledHomotopy,
    x::TaylorVector{D},
    t::Number,
    p::Nothing = nothing,
) where {M,D}
    _taylor!_impl(T, M, D - 1, -1; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledHomotopy,
    x::TaylorVector{D},
    t::Number,
    p::AbstractVector,
) where {M,D}
    _taylor!_impl(T, M, D - 1, 0; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledHomotopy,
    x::TaylorVector{D},
    t::Number,
    p::TaylorVector{DP},
) where {M,D,DP}
    _taylor!_impl(T, M, D - 1, DP - 1; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledHomotopy,
    x::AbstractVector,
    t::Number,
    p::Nothing = nothing,
) where {M}
    _taylor!_impl(T, M, 0, -1; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{M},
    T::CompiledHomotopy,
    x::AbstractVector,
    t::Number,
    p::AbstractVector,
) where {M}
    _taylor!_impl(T, M, 0, 0; highest_order_only = true)
end
@generated function taylor!(
    u::AbstractVector,
    ::Val{1},
    T::CompiledHomotopy,
    x::AbstractVector,
    t::Number,
    p::TaylorVector{DP},
) where {M,DP}
    _taylor!_impl(T, M, 0, DP - 1; highest_order_only = true)
end

# non-inplace

"""
    to_smallest_eltype(A::AbstractArray)

Convert an array to the smallest eltype such that all elements still fit.

## Example
```julia
typeof(to_smallest_elype(Any[2,3])) == Vector{Int}
```
"""
function to_smallest_eltype(A::AbstractArray)
    T = typeof(first(A))
    for a in A
        T = promote_type(T, typeof(a))
    end
    convert.(T, A)
end

function evaluate(T::CompiledSystem, x, p = nothing)
    to_smallest_eltype(evaluate!(Vector{Any}(undef, size(T, 1)), T, x, p))
end
function evaluate(T::CompiledHomotopy, x, t, p = nothing)
    to_smallest_eltype(evaluate!(Vector{Any}(undef, size(T, 1)), T, x, t, p))
end

(T::CompiledSystem)(x, p = nothing) = evaluate(x, p)
(T::CompiledHomotopy)(x, t, p = nothing) = evaluate(x, t, p)

function jacobian(T::CompiledSystem, x, p = nothing)
    n, m = size(T)
    U = Matrix{Any}(undef, n, m)
    to_smallest_eltype(jacobian!(U, T, x, p))
end
function jacobian(T::CompiledHomotopy, x, t, p = nothing)
    n, m = size(T)
    U = Matrix{Any}(undef, n, m)
    to_smallest_eltype(jacobian!(U, T, x, t, p))
end

function evaluate_and_jacobian(T::CompiledSystem, x, p = nothing)
    n, m = size(T)
    u = Vector{Any}(undef, n)
    U = Matrix{Any}(undef, n, m)
    evaluate_and_jacobian!(u, U, T, x, p)
    to_smallest_eltype(u), to_smallest_eltype(U)
end
function evaluate_and_jacobian(T::CompiledHomotopy, x, t, p = nothing)
    n, m = size(T)
    u = Vector{Any}(undef, n)
    U = Matrix{Any}(undef, n, m)
    evaluate_and_jacobian!(u, U, T, x, t, p)
    to_smallest_eltype(u), to_smallest_eltype(U)
end

ext: jl


text:
module EPOCHInput

export parse_input

using Unitful
using PhysicalConstants.CODATA2018: c_0, ε_0, μ_0, m_e, e
using BangBang

include("parse.jl")
include("constant.jl")
include("control.jl")
include("no_globals.jl")
include("no_parsing.jl")
include("output.jl")

end


text:
"""
Package: Forecast

Store results from the function stl

# Arguments
    `decomposition::DataFrame`    A time array with three time series from a fitted STL model
    `call::String`                method called to generate ta
"""
mutable struct STL{T<:DataFrame}
    decomposition::T  # A DataFrame with three time series from a fitted STL model
    call::String      # method called
end

function Base.show(io::IO, stlx::STL)
    printstyled("STL Object:",bold=true,color=:underline)
    print(" ", stlx.call)
end


text:
using Dimensionless, Unitful
using Test

@testset "DimBasis construction" begin
    # Check invalid basis construction exceptions
    @test_throws ErrorException DimBasis(9.81u"m/s^2", 1u"mm", 1u"s")
    @test_throws ErrorException DimBasis(9.81u"m", 9.81u"m")
    @test_throws ErrorException DimBasis(9.81u"m/s^2")

    # Check exception in dim_matrix() in case of unvalid dimension
    basis = DimBasis(9.81u"m/s^2", 6371u"km", 1420788u"kg", 1u"mol/g")
    @test_throws ErrorException dim_matrix(basis.basis_dims, 1u"A")

    @test basis.dim_mat == Rational.([1 1 0 0; -2 0 0 0; 0 0 1 -1; 0 0 0 1])
    for dim in [:Mass, :Length, :Time, :Amount]
        @test Unitful.Dimension{dim} in basis.basis_dims
    end
    for dim in [:Temperature, :Current, :Luminosity]
        @test !(Unitful.Dimension{dim} in basis.basis_dims)
    end

    # Check all constructor combinations for DimBasis
    basis_vector_names_template = ["α", "β", "γ"]
    basis_vectors_template = [1u"s^2", 2u"mm^-1", 3u"bar^3"]
    for (broadcast_fcn, isnamed) in Iterators.product((identity, unit, dimension), (false, true))
        # Construct bases
        input_vectors = broadcast(broadcast_fcn, basis_vectors_template)
        basis = isnamed ?
        DimBasis([Pair(basis_vector_names_template[i], input_vectors[i]) for i in eachindex(basis_vectors_template)]...) :
        DimBasis(input_vectors...)

        # Check validity
        isnamed && @test basis.basis_vector_names == basis_vector_names_template
        @test basis.dim_mat == Rational.([2 0 -6; 0 -1 -3; 0 0 3])
        @test basis.basis_dims == [Unitful.Dimension{:Time}, Unitful.Dimension{:Length}, Unitful.Dimension{:Mass}]
        @test basis.basis_vectors == input_vectors
    end
end

@testset "Change of basis" begin
    # Test change of basis and broadcast
    dimless_dimful_basis = DimBasis(1u"mA", 2u"K", 3u"Pa*s", 4u"h", 5u"kg")
    vars = [0u"A/m^2", 1u"nm", 2u"K"]
    dim_less = dimensionless.(vars, dimless_dimful_basis)
    @test dimensionful.(dim_less, [u"A/m^2", u"nm", u"K"], dimless_dimful_basis) == vars

    # Test counting of dimensions and dimensionless variables for named basis
    quantities_named = ["a"=>100.0u"m", "b"=>100u"kg", "c"=>0u"s"]
    units_named = [Pair(var.first, unit(var.second)) for var in quantities_named]
    dims_named = [Pair(var.first, dimension(var.second)) for var in quantities_named]
    @test number_of_dimensions(quantities_named...) == 3
    @test number_of_dimensions(units_named...) == 3
    @test number_of_dimensions(dims_named...) == 3
    @test number_of_dimensionless(quantities_named...) == 0
    @test number_of_dimensionless(units_named...) == 0
    @test number_of_dimensionless(dims_named...) == 0

    # Test counting of dimensions and dimensionless variables for unnamed basis
    quantities_unnamed = [var.second for var in quantities_named]
    units_unnamed = [var.second for var in units_named]
    dims_unnamed = [var.second for var in dims_named]
    @test number_of_dimensions(quantities_unnamed...) == 3
    @test number_of_dimensions(units_unnamed...) == 3
    @test number_of_dimensions(dims_unnamed...) == 3
    @test number_of_dimensionless(quantities_unnamed...) == 0
    @test number_of_dimensionless(units_unnamed...) == 0
    @test number_of_dimensionless(dims_unnamed...) == 0

    # change_basis(quantity/unit, old_b, new_b) for named and unnamed quantity bases
    basis_named = DimBasis("x"=>1u"m", "y"=>1u"kg", "z"=>1u"s")
    new_basis_named = DimBasis("x_new"=>2u"m", "y_new"=>2u"kg", "z_new"=>2u"s")
    @test Dimensionless.current_to_new_fac(Unitful.𝐋*Unitful.𝐓*Unitful.𝐌, basis_named, new_basis_named) == 2^3
    @test change_basis(100u"cm*kg*s", basis_named, new_basis_named) == 8u"m*kg*s"
    @test change_basis(u"m*kg*s", basis_named, new_basis_named) == 8

    basis_unnamed = DimBasis(1u"m", 1u"kg", 1u"s")
    new_basis_unnamed = DimBasis(2u"m", 2u"kg", 2u"s")
    @test Dimensionless.current_to_new_fac(Unitful.𝐋*Unitful.𝐓*Unitful.𝐌, basis_unnamed, new_basis_unnamed) == 2^3
    @test change_basis(100u"cm*kg*s", basis_unnamed, new_basis_unnamed) == 8u"m*kg*s"
    @test change_basis(u"m*kg*s", basis_unnamed, new_basis_unnamed) == 8
end

@testset "Utils" begin
    var = "L"=>0.2u"m"
    quantities_named = ["ρ"=>1400.0u"kg/m^3", "v"=>0.5u"m/s", "η"=>(10^4)u"Pa*s"]
    units_named = [Pair(var.first, unit(var.second)) for var in quantities_named]
    dims_named = [Pair(var.first, dimension(var.second)) for var in quantities_named]
    res_desired = "L ρ v η^-1"
    res_buf = IOBuffer()

    # Test dimensionless string for named basis
    print_dimensionless(res_buf, var, DimBasis(quantities_named...))
    @test String(take!(res_buf)) == res_desired
    print_dimensionless(res_buf, var, DimBasis(units_named...))
    @test String(take!(res_buf)) == res_desired
    print_dimensionless(res_buf, var, DimBasis(dims_named...))
    @test String(take!(res_buf)) == res_desired
end


text:
"""
    check_feasibility(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point) where T<:AbstractNodeLabel

Evaluate the constraints defined by `elemid` of type `elem` at `bus` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if there is not violated constraints.

"""
function check_feasibility(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:AbstractNodeLabel
    return
end

"""
    check_feasibility(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point, epsilon::Float64) where T<:AbstractLinkLabel

Evaluate the constraints defined by `elemid` of type `elem` at `link` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if there is not violated constraints.

"""

function check_feasibility(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point, epsilon::Float64) where T<:AbstractLinkLabel
    return
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
@testset "FerriteProperties" begin
test_spldata = PlanerTransformer.SplInput((
  25e3,(100,14),(300,270),
  100e3,(62,20,),(240,750),
  200e3,(38,8),(190,1100),
  400e3,(20,17),(140,1700),
  700e3,(12,20),(80,1750)
  ))
test_fp = FerriteProperties("3f3",0.2e6,0.5e6,25,100,
                            BHloop(15,.44,.15),
                            BHloop(10,.345,.120),
                            SpecificPowerLossData(test_spldata)) # 3f3
@test specific_power_loss(test_fp.spl_hot,0.2,15e3)≈72926.68854848042
@test specific_power_loss(test_fp.spl_hot,0.1,50e3)≈24158.470844381052
@test specific_power_loss(test_fp.spl_hot,0.07,150e3)≈37883.89793616807
@test specific_power_loss(test_fp.spl_hot,0.05,300e3)≈52475.72366004552
@test specific_power_loss(test_fp.spl_hot,0.02,600e3)≈42277.685735775805
@test specific_power_loss(test_fp.spl_hot,0.01,1e6)≈51375.309340858454
@test specific_power_loss(test_fp.spl_hot,0.04,400e3)≈87672.48079453735
@test specific_power_loss(test_fp,0.04,400e3)≈87672.48079453735
@test length(ferrite_dict)>0
# need better inverse function at some point
@test flux_density(test_fp.spl_hot,42277.685735775805,600e3)≈0.019989634564937644
@test flux_density(test_fp,42277.685735775805,600e3)≈0.019989634564937644

@test_throws ArgumentError SpecificPowerLossData((),())
@test_throws ArgumentError SpecificPowerLossData((1),(2,3))
@test_throws ArgumentError SpecificPowerLossData((1,2,3),((1,10),(2,20)))

@test_throws ArgumentError SpecificPowerLossData(PlanerTransformer.SplInput((1,2,3,4)))
end # testset


text:
module MsgPack

using Compat

export pack, unpack, Ext
import Base: ==

const INT_FP   = 0x00 # - 0xf7
const MAP_F    = 0x80 # - 0x8f
const ARR_F    = 0x90 # - 0x9f
const STR_F    = 0xa0 # - 0xbf
const EXT_F    = 0xd4 # - 0xd8

const NIL      = 0xc0
const UNUSED   = 0xc1
const FALSE    = 0xc2
const TRUE     = 0xc3
const BIN_8    = 0xc4
const BIN_16   = 0xc5
const BIN_32   = 0xc6
const EXT_8    = 0xc7
const EXT_16   = 0xc8
const EXT_32   = 0xc9
const FLOAT_32 = 0xca
const FLOAT_64 = 0xcb
const UINT_8   = 0xcc
const UINT_16  = 0xcd
const UINT_32  = 0xce
const UINT_64  = 0xcf
const INT_8    = 0xd0
const INT_16   = 0xd1
const INT_32   = 0xd2
const INT_64   = 0xd3
const STR_8    = 0xd9
const STR_16   = 0xda
const STR_32   = 0xdb
const ARR_16   = 0xdc
const ARR_32   = 0xdd
const MAP_16   = 0xde
const MAP_32   = 0xdf

const INT_FN   = 0xe0 # - 0xff

immutable Ext
    typecode::Int8
    data::Vector{UInt8}

    function Ext(t::Integer, d::Vector{UInt8}; impltype=false)
        # -128 to -1 reserved for implementation
        if -128 <= t <= -1
            impltype || error("MsgPack Ext typecode -128 through -1 reserved by implementation")
        elseif !(0 <= t <= 127)
            error("MsgPack Ext typecode must be in the range [-128, 127]")
        end

        new(t, d)
    end
end
==(a::Ext, b::Ext) = a.typecode == b.typecode && a.data == b.data

# return Ext where Ext.data is a serialized object
function extserialize(t::Integer, d)
    i = IOBuffer()
    serialize(i, d)
    return Ext(t, takebuf_array(i))
end

# return (typecode, object) from an Ext where Ext.data is a serialized object
extdeserialize(e::Ext) = (e.typecode, deserialize(IOBuffer(e.data)))

readn(s, t) = ntoh(read(s, t))
readi(s, t) = @compat Int64(readn(s, t))

readu64(s, t) = begin
    v = @compat UInt64(readn(s, t))
    if v > 2^63-1
        v
    else
        @compat Int64(v)
    end
end

const DISPATCH =
    @compat Dict( NIL      => s -> nothing
                 ,UNUSED   => s -> error("unused")
                 ,FALSE    => s -> false
                 ,TRUE     => s -> true
                 ,BIN_8    => s -> unpack_bin(s, readn(s, UInt8))
                 ,BIN_16   => s -> unpack_bin(s, readn(s, UInt16))
                 ,BIN_32   => s -> unpack_bin(s, readn(s, UInt32))
                 ,EXT_8    => s -> unpack_ext(s, readn(s, UInt8))
                 ,EXT_16   => s -> unpack_ext(s, readn(s, UInt16))
                 ,EXT_32   => s -> unpack_ext(s, readn(s, UInt32))
                 ,FLOAT_32 => s -> readn(s, Float32)
                 ,FLOAT_64 => s -> readn(s, Float64)
                 ,UINT_8   => s -> readi(s, UInt8)
                 ,UINT_16  => s -> readi(s, UInt16)
                 ,UINT_32  => s -> readi(s, UInt32)
                 ,UINT_64  => s -> readu64(s, UInt64)
                 ,INT_8    => s -> readi(s, Int8)
                 ,INT_16   => s -> readi(s, Int16)
                 ,INT_32   => s -> readi(s, Int32)
                 ,INT_64   => s -> readi(s, Int64)
                 ,STR_8    => s -> unpack_str(s, readn(s, UInt8))
                 ,STR_16   => s -> unpack_str(s, readn(s, UInt16))
                 ,STR_32   => s -> unpack_str(s, readn(s, UInt32))
                 ,ARR_16   => s -> unpack_arr(s, readn(s, UInt16))
                 ,ARR_32   => s -> unpack_arr(s, readn(s, UInt32))
                 ,MAP_16   => s -> unpack_map(s, readn(s, UInt16))
                 ,MAP_32   => s -> unpack_map(s, readn(s, UInt32))
                )

unpack(s) = unpack(IOBuffer(s))
unpack(s::IO) = begin
    b = read(s, UInt8)

    if b <= 0x7f
        # positive fixint
        @compat Int64(b)

    elseif b <= 0x8f
        # fixmap
        unpack_map(s, b $ MAP_F)

    elseif b <= 0x9f
        # fixarray
        unpack_arr(s, b $ ARR_F)

    elseif b <= 0xbf
        # fixstr
        unpack_str(s, b $ STR_F)

    elseif 0xd4 <= b <= 0xd8
        # fixext
        unpack_ext(s, 2^(b - EXT_F))

    elseif b <= 0xdf
        DISPATCH[b](s)

    else
        # negative fixint
        @compat Int64(reinterpret(Int8, b))
    end
end

unpack_map(s, n) = begin
    out = Dict()
    for i in 1:n
        k = unpack(s)
        v = unpack(s)
        out[k] = v
    end
    out
end

unpack_arr(s, n) = begin
    out = Array(Any, n)
    for i in 1:n
        out[i] = unpack(s)
    end
    out
end

unpack_str(s, n) = utf8(readbytes(s, n))
unpack_ext(s, n) = Ext(read(s, Int8), readbytes(s, n), impltype=true)
unpack_bin(s, n) = readbytes(s, n)

wh(io, head, v) = begin
    write(io, head)
    write(io, hton(v))
end

pack(v) = begin
    s = IOBuffer()
    pack(s, v)
    takebuf_array(s)
end

@compat pack(s, ::Void) = write(s, NIL)
pack(s, v::Bool)   = if v write(s, TRUE) else write(s, FALSE) end

pack(s, v::Integer) = begin
    if v < 0
        if v >= -32
            write(s, @compat Int8(v))
        elseif v >= -2^7
            wh(s, INT_8, @compat Int8(v))
        elseif v >= -2^15
            wh(s, INT_16, @compat Int16(v))
        elseif v >= -2^31
            wh(s, INT_32, @compat Int32(v))
        elseif v >= -2^63
            wh(s, INT_64, @compat Int64(v))
        else
            error("MsgPack signed int overflow")
        end
    else
        if v <= 127
            write(s, @compat UInt8(v))
        elseif v <= 2^8-1
            wh(s, UINT_8, @compat UInt8(v))
        elseif v <= 2^16-1
            wh(s, UINT_16, @compat UInt16(v))
        elseif v <= 2^32-1
            wh(s, UINT_32, @compat UInt32(v))
        elseif v <= @compat UInt64(2)^64-1
            wh(s, UINT_64, @compat UInt64(v))
        else
            error("MsgPack unsigned int overflow")
        end
    end
end

pack(s, v::Float32) = wh(s, 0xca, v)
pack(s, v::Float64) = wh(s, 0xcb, v)

# str format
pack(s, v::AbstractString) = begin
    n = sizeof(v)
    if n < 2^5
        write(s, STR_F | @compat UInt8(n))
    ## Note: with this section commented out, we do not have
    ##       the most compact format for a string.  However,
    ##       the string is still in spec, and some other
    ##       msgpack libaries (*ahem* Python) can't decode
    ##       strings created with this rule.
    #elseif n < 2^8
    #    wh(s, 0xd9, @compat UInt8(n))
    elseif n < 2^16
        wh(s, 0xda, @compat UInt16(n))
    elseif n < 2^32
        wh(s, 0xdb, @compat UInt32(n))
    else
        error("MsgPack str overflow: ", n)
    end
    write(s, v)
end

# ext format
pack(s, v::Ext) = begin
    n = sizeof(v.data)
    if n == 1
        write(s, 0xd4)
    elseif n == 2
        write(s, 0xd5)
    elseif n == 4
        write(s, 0xd6)
    elseif n == 8
        write(s, 0xd7)
    elseif n == 16
        write(s, 0xd8)
    elseif n < 2^8
        wh(s, 0xc7, @compat UInt8(n))
    elseif n < 2^16
        wh(s, 0xc8, @compat UInt16(n))
    elseif n < 2^32
        wh(s, 0xc9, @compat UInt32(n))
    else
        error("MsgPack ext overflow: ", n)
    end
    write(s, v.typecode)
    write(s, v.data)
end

# bin format
pack(s, v::Vector{UInt8}) = begin
    n = length(v)
    if n < 2^8
        wh(s, 0xc4, @compat UInt8(n))
    elseif n < 2^16
        wh(s, 0xc5, @compat UInt16(n))
    elseif n < 2^32
        wh(s, 0xc6, @compat UInt32(n))
    else
        error("MsgPack bin overflow: ", n)
    end
    write(s, v)
end

# Simple arrays
pack(s, @compat v::Union{Vector, Tuple}) = begin
    n = length(v)
    if n < 2^4
        write(s, ARR_F | @compat UInt8(n))
    elseif n < 2^16
        wh(s, 0xdc, @compat UInt16(n))
    elseif n < 2^32
        wh(s, 0xdd, @compat UInt32(n))
    else
        error("MsgPack array overflow: ", n)
    end

    for x in v
        pack(s, x)
    end
end

# Maps
pack(s, v::Dict) = begin
    n = length(v)
    if n < 2^4
        write(s, MAP_F | @compat UInt8(n))
    elseif n < 2^16
        wh(s, 0xde, @compat UInt16(n))
    elseif n < 2^32
        wh(s, 0xdf, @compat UInt32(n))
    else
        error("MsgPack map overflow: ", n)
    end

    for (k, x) in v
        pack(s, k)
        pack(s, x)
    end
end

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -6.148474362391325 <= q <= 10.677081718106185)

                     add_NL_constraint(m, :(gelu(-0.2518902526786948 + 0.9847866884384731*gelu(0.2752793536861313 + -0.1568397657479923*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.7649756572019379*gelu(0.9918790053549298 + 0.85981905523253*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.10087195663512549*gelu(-0.8202416682813327 + -0.17766510212211006*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + gelu(0.638435553115452 + 0.09453389081965424*gelu(0.2752793536861313 + -0.1568397657479923*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.07927014075141203*gelu(0.9918790053549298 + 0.85981905523253*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.2624391124261729*gelu(-0.8202416682813327 + -0.17766510212211006*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + gelu(0.5273095272255199 + 0.4025366346817978*gelu(0.2752793536861313 + -0.1568397657479923*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.585290488444234*gelu(0.9918790053549298 + 0.85981905523253*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.1647823489958915*gelu(-0.8202416682813327 + -0.17766510212211006*gelu(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*gelu(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*gelu(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    QuadRefinement()

Refinement of polygonal meshes into quadrangles.
A n-gon is subdivided into n quadrangles.
"""
struct QuadRefinement <: RefinementMethod end

function refine(mesh, ::QuadRefinement)
  # retrieve geometry and topology
  points = vertices(mesh)
  connec = topology(mesh)

  # convert to half-edge structure
  t = convert(HalfEdgeTopology, connec)

  # add centroids of elements
  ∂₂₀ = Boundary{2,0}(t)
  epts = map(1:nelements(t)) do elem
    ps = view(points, ∂₂₀(elem))
    cₒ = sum(coordinates, ps) / length(ps)
    Point(cₒ)
  end

  # add midpoints of edges
  ∂₁₀ = Boundary{1,0}(t)
  fpts = map(1:nfacets(t)) do edge
    ps = view(points, ∂₁₀(edge))
    cₒ = sum(coordinates, ps) / length(ps)
    Point(cₒ)
  end

  # original vertices
  vpts = points

  # new points in refined mesh
  newpoints = [vpts; epts; fpts]

  offset₁ = length(vpts)
  offset₂ = offset₁ + length(epts)

  # connect vertices into new quadrangles
  ∂₂₁ = Boundary{2,1}(t)
  newconnec = Connectivity{Quadrangle,4}[]
  for elem in 1:nelements(t)
    verts = CircularVector(∂₂₀(elem))
    edges = CircularVector(∂₂₁(elem))
    for i in 1:length(edges)
      u = elem       + offset₁
      v = edges[i]   + offset₂
      w = verts[i+1]
      z = edges[i+1] + offset₂
      quad = connect((u, v, w, z))
      push!(newconnec, quad)
    end
  end

  SimpleMesh(newpoints, newconnec)
end


text:
module LPSjl

include("BayesianNetwork/BayesianNetwork.jl")

include("ParameterLearning/ParameterLearning.jl")

export
    BayesianNetworks,
    ParameterLearning

end


text:
using POMDPs

POMDPs.add("GenerativeModels")
POMDPs.add("POMDPToolbox")
POMDPs.add("MCTS")
POMDPs.add("POMCP")

try
    Pkg.clone("https://github.com/tawheeler/Vec.jl.git")
catch
    println("already installed.")
end

try
    Pkg.clone("https://github.com/tawheeler/AutomotiveDrivingModels.jl.git")
catch
    println("already installed.")
end

try
    Pkg.clone("https://github.com/tawheeler/AutoViz.jl.git")
catch
    println("already installed.")
end

try
    Pkg.clone("https://github.com/slundberg/PmapProgressMeter.jl.git")
catch
    println("already installed.")
end

try
    Pkg.clone("https://github.com/JuliaPlots/StatPlots.jl.git")
catch
    println("already installed.")
end


text:
using Statistics
using StatsBase: sample
using Revise
using EvoTrees
using BenchmarkTools

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

#############################
# CPU - linear
#############################
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=0.5, nbins=64)

# asus laptopt: for 1.25e6 no eval: 12.805 s (196419 allocations: 5.92 GiB)
@time model = fit_evotree(params1, X_train, Y_train);
@btime model = fit_evotree($params1, $X_train, $Y_train);
@time pred_train = predict(model, X_train);
@btime pred_train = predict(model, X_train);
gain = importance(model, 1:100)

@time model, cache = EvoTrees.init_evotree(params1, X_train, Y_train);
@time EvoTrees.grow_evotree!(model, cache);

using BSON: @save, @load
@save "model.bson" model
@load "model.bson" model
pred_train = predict(model, X_train)

#############################
# CPU - Logistic
#############################
params1 = EvoTreeRegressor(T=Float64,
    loss=:logistic, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=4, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64)
@time model = fit_evotree(params1, X_train, Y_train);
@time pred = predict(model, X_train);

#############################
# CPU - Gaussian
#############################
params1 = EvoTreeGaussian(T=Float64,
    loss=:gaussian, metric=:none,
    nrounds=100,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

# Asus laptop: 25.726237 seconds (205.93 k allocations: 8.028 GiB, 3.04% gc time)
@time model = fit_evotree(params1, X_train, Y_train);
# Asus laptop: 1.667185 seconds (2.41 k allocations: 1.658 GiB)
@time model, cache = EvoTrees.init_evotree(params1, X_train, Y_train);

################################
# GPU - Linear
################################
# train model
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=0.5, nbins=64,
    device="gpu")

# Asus laptop:  10.103512 seconds (20.83 M allocations: 7.921 GiB, 12.15% gc time)
@time model = EvoTrees.fit_evotree(params1, X_train, Y_train);
@btime model = EvoTrees.fit_evotree(params1, X_train, Y_train);
# @time model, cache = EvoTrees.init_evotree_gpu(params1, X_train, Y_train);
# @time EvoTrees.grow_evotree!(model, cache);

# X_train_32 = Float32.(X_train)
@time pred_train = EvoTrees.predict(model, X_train);
@btime pred_train = EvoTrees.predict(model, X_train);
mean(pred_train)

################################
# GPU - Logistic
################################
# train model
params1 = EvoTreeRegressor(T=Float32,
    loss=:logistic, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64,
    device="gpu")
@time model = fit_evotree(params1, X_train, Y_train);
@time pred_train = predict(model, X_train)

################################
# GPU - Gaussian
################################
params1 = EvoTreeGaussian(T=Float64,
    loss=:gaussian, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32,
    device="gpu")
# Asus laptop: 8.214754 seconds (14.79 M allocations: 6.101 GiB, 12.42% gc time)
@time model = EvoTrees.fit_evotree(params1, X_train, Y_train);
# Auss laptop:  1.888472 seconds (8.40 k allocations: 1.613 GiB, 14.86% gc time)
@time model, cache = EvoTrees.init_evotree(params1, X_train, Y_train);

############################
# xgboost
############################
using XGBoost
num_round = 100
param = ["max_depth" => 5,
         "eta" => 0.05,
         "objective" => "reg:squarederror",
         "print_every_n" => 5,
         "subsample" => 0.5,
         "colsample_bytree" => 0.5,
         "tree_method" => "hist",
         "nthread" => 16,
         "max_bin" => 32]
metrics = ["rmse"]
@time xgboost(X_train, num_round, label=Y_train, param=param, metrics=metrics, silent=1);
@time dtrain = DMatrix(X_train, label=Y_train)
@time model_xgb = xgboost(dtrain, num_round, param=param, silent=1);
@btime model_xgb = xgboost(dtrain, num_round, param=param, silent=1);
@time pred_train = XGBoost.predict(model_xgb, X_train)

@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=9999, early_stopping_rounds=9999);
@btime model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=9999, early_stopping_rounds=9999);

@time model = fit_evotree(params1, X_train, Y_train, early_stopping_rounds=10);
@time model = fit_evotree(params1, X_train, Y_train, print_every_n=2);

# @time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 5);
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval);
@time pred_train = predict(model, X_train)

@code_warntype predict(model, X_train)
@time pred = zeros(SVector{1,Float64}, size(X_train, 1))
@time EvoTrees.predict!(pred, model.trees[2], X_train)

@time predict(model, X_train)
@btime pred_train = predict($model, $X_train)
mean(abs.(pred_train .- Y_train))

# logistic
params1 = EvoTreeRegressor(
    loss=:logistic, metric=:logloss,
    nrounds=100,
    λ=0.0f0, γ=0.0f0, η=0.1f0,
    max_depth=6, min_weight=1.0f0,
    rowsample=0.5f0, colsample=0.5f0, α=0.5f0, nbins=32)
@time model = fit_evotree(params1, X_train, Y_train);
@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=10)
@time pred_train = predict(model, X_train)

# Quantile
params1 = EvoTreeRegressor(
    loss=:quantile, metric=:quantile, α=0.80f0,
    nrounds=100,
    λ=0.1f0, γ=0.0f0, η=0.1f0,
    max_depth=6, min_weight=1.0f0,
    rowsample=0.5f0, colsample=0.5f0, nbins=32)
@time model = fit_evotree(params1, X_train, Y_train);
@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=10)
@time pred_train = predict(model, X_train)

# gaussian
params1 = EvoTreeGaussian(
    loss=:gaussian, metric=:gaussian,
    nrounds=100, α=0.5f0,
    λ=0.0f0, γ=0.0f0, η=0.1f0,
    max_depth=6, min_weight=10.0f0,
    rowsample=0.5f0, colsample=0.5f0, nbins=32)
@time model = fit_evotree(params1, X_train, Y_train);
@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=10)
@time pred_train = predict(model, X_train)

# softmax
params1 = EvoTreeClassifier(
    loss=:softmax, metric=:mlogloss,
    nrounds=100, α=0.5f0,
    λ=0.0f0, γ=0.0f0, η=0.1f0,
    max_depth=6, min_weight=10.0f0,
    rowsample=0.5f0, colsample=0.5f0, nbins=32)

Y_train_int = UInt32.(round.(Y_train * 2) .+ 1)
Y_eval_int = UInt32.(round.(Y_eval * 2) .+ 1)
Y_train_int = Int.(Y_train_int)
@time model = fit_evotree(params1, X_train, Y_train_int, print_every_n=10);
@time model = fit_evotree(params1, X_train, Y_train_int, X_eval=X_eval, Y_eval=Y_eval_int, print_every_n=10)
@time pred_train = predict(model, X_train)


text:
__precompile__()

module DiffEqOptimalControl

# package code goes here

end # module


text:
using LinearAlgebra: Symmetric, dot

@testset "Kunisch-Rendl" begin
    using NLS_Solver: BoundConstraint_Enum, INACTIVE_BC, ACTIVE_LB, ACTIVE_UB, restrict_to_inactive!

    n=10
    A=[Rational{BigInt}(1,i+j-1) for i in 1:n, j in 1:n]

    Z=BoundConstraint_Enum[INACTIVE_BC,
                           ACTIVE_UB,
                           INACTIVE_BC,
                           INACTIVE_BC,
                           ACTIVE_UB,
                           ACTIVE_UB,
                           ACTIVE_LB,
                           ACTIVE_LB,
                           ACTIVE_LB,
                           INACTIVE_BC]

    lb=zeros(Rational{BigInt},n)
    ub=ones(Rational{BigInt},n)

    Q_result = Rational{Int64}[1//1 0//1 1//3 1//4 0//1 0//1 0//1 0//1 0//1 1//10; 0//1 1//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1; 1//3 0//1 1//5 1//6 0//1 0//1 0//1 0//1 0//1 1//12; 1//4 0//1 1//6 1//7 0//1 0//1 0//1 0//1 0//1 1//13; 0//1 0//1 0//1 0//1 1//1 0//1 0//1 0//1 0//1 0//1; 0//1 0//1 0//1 0//1 0//1 1//1 0//1 0//1 0//1 0//1; 0//1 0//1 0//1 0//1 0//1 0//1 1//1 0//1 0//1 0//1; 0//1 0//1 0//1 0//1 0//1 0//1 0//1 1//1 0//1 0//1; 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 1//1 0//1; 1//10 0//1 1//12 1//13 0//1 0//1 0//1 0//1 0//1 1//19]
    q_result = Rational{Int64}[28//15, -1//1, 197//56, 1597//360, -1//1, -1//1, 0//1, 0//1, 0//1, 23629//2310]

    # Store upper part
    #
    Q=Symmetric(A,:U)
    q=Rational{BigInt}[i for i in 1:n]

    restrict_to_inactive!(Q, q, Z, lb, ub)

    @test Q==Q_result
    @test q==q_result

    # Store lower part
    #
    Q=Symmetric(A,:L)
    q=Rational{BigInt}[i for i in 1:n]

    restrict_to_inactive!(Q, q, Z, lb, ub)

    @test Q==Q_result
    @test q==q_result

end

@testset "problem 1 uplow = $UL" for UL ∈ (:U,:L)
    using NLS_Solver: check_first_order, Kunisch_Rendl_Conf

    n=5
    A=[Rational{BigInt}(1,i+j-1) for i in 1:n, j in 1:n]

    Q=Symmetric(A,UL)
    q=Rational{BigInt}[-1 for i in 1:n]
    bc=BoundConstraints(Int,n)
    x_init=zeros(Int,n)

    conf = Kunisch_Rendl_Conf()
    result = solve(Q,q,x_init,bc,conf)

    x_sol = solution(result)
    @test converged(result)
    @test iteration_count(result)==7
    @test objective_value(result) ≈ dot(x_sol,Q*x_sol)/2+dot(q,x_sol)
    @test 1+check_first_order(Q,q,x_sol,bc) ≈ 1+0
end

@testset "problem 2" begin
     using NLS_Solver: check_first_order, Kunisch_Rendl_Conf, multiplier_τ

    Q=Symmetric(Float64[[30 20 15]
                        [20 15 12]
                        [15 12 10]],:L)

    q=-Float64[1:3;]
    bc=BoundConstraints(zeros(3),Float64[1:3;])
    x_init = zeros(3)

    conf = Kunisch_Rendl_Conf()
    result = solve(Q,q,x_init,bc,conf)

    x_sol = solution(result)
    τ = multiplier_τ(result)

    @test converged(result)
    @test iteration_count(result) == 4
    @test τ[1] ≈ -3.5
    @test τ[2] ≈ -1.6
    @test τ[3] == 0
    @test objective_value(result) ≈ dot(x_sol,Q*x_sol)/2+dot(q,x_sol)
    @test 1+check_first_order(Q,q,x_sol,bc) ≈ 1+0

end


text:
using Documenter
using DocumenterCitations
using LatSpec

bib = CitationBibliography(joinpath(@__DIR__,"LatSpec.bib"))
const PAGES = [
               "Home"       => "index.md",
               "API"        => "api.md",
               "References" => "references.md"
              ]

makedocs(bib;
    authors = "Bernd Riederer and Fabian Zierler",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets=String[],
    ),
    pages = PAGES,
    repo="https://github.com/bernd1995/LatSpec.jl/blob/{commit}{path}#{line}",
    sitename="LatSpec.jl"
)

deploydocs(;
    repo = "github.com/bernd1995/LatSpec.jl",
    target = "build",
    push_preview = true
)


text:
#==========================================================================================
                                Constructors
==========================================================================================#
function SymSemiseparableMatrix(U::AbstractArray, V::AbstractArray)
    if size(U,1) == size(V,1) && size(U,2) == size(V,2)
        return SymSemiseparableMatrix(size(U,2),size(U,1),U,V)
    else
        error("Dimension mismatch between generators U and V")
    end
end
function SymSemiseparableMatrix(L::SymSemiseparableCholesky)
    return SymSemiseparableMatrix(L.n, L.p, L.Ut, ss_create_v(L.Ut, L.Wt))
end
#==========================================================================================
                    Defining Matrix Properties / Overloading Base
==========================================================================================#
Base.Matrix(K::SymSemiseparableMatrix)   = tril(K.Ut'*K.Vt) + triu(K.Vt'*K.Ut,1)
Base.size(K::SymSemiseparableMatrix)     = (K.n,K.n)
function Base.getindex(K::SymSemiseparableMatrix{T}, i::Int, j::Int) where T
	i > j && return dot(K.Ut[:,i],K.Vt[:,j])
	return dot(K.Ut[:,j],K.Vt[:,i])
end
function Base.propertynames(F::SymSemiseparableMatrix, private::Bool=false)
   return (private ? fieldnames(typeof(F)) : ())
end
#==========================================================================================
                    Overloading LinearAlgebra routines
==========================================================================================#
LinearAlgebra.cholesky(K::SymSemiseparableMatrix) = SymSemiseparableCholesky(K)
LinearAlgebra.logdet(K::SymSemiseparableMatrix)   = 2.0*logdet(cholesky(K))
LinearAlgebra.det(K::SymSemiseparableMatrix)      = det(cholesky(K))^2
function LinearAlgebra.cholesky(K::SymSemiseparableMatrix, sigma)
	Wt, dbar = dss_create_wdbar(K.Ut,K.Vt,ones(K.n)*sigma)
	return DiaSymSemiseparableCholesky(K.n,K.p,K.Ut,Wt,dbar)
end
#==========================================================================================
                        Defining multiplication and inverse
==========================================================================================#
function LinearAlgebra.mul!(y::AbstractArray, L::SymSemiseparableMatrix, x::AbstractArray)
    ss_mul_mat!(y, L.Ut, L.Vt, x)
    return y
end
function LinearAlgebra.mul!(y::AbstractArray, L::Adjoint{<:Any,<:SymSemiseparableMatrix}, x::AbstractArray)
    ss_mul_mat!(y, L.parent.Ut, L.parent.Vt, x)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray,K::SymSemiseparableMatrix, x::AbstractArray)
	L  = cholesky(K)
    y .= L'\(L\x)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray,K::Adjoint{<:Any,<:SymSemiseparableMatrix}, x::AbstractArray)
	L  = cholesky(K.parent)
	y .= L'\(L\x)
    return y
end
function (Base.:\)(K::SymSemiseparableMatrix, x::AbstractVecOrMat)
    y  = similar(x)
	L  = cholesky(K)
    y .= L'\(L\x)
    return y
end
function (Base.:\)(K::Adjoint{<:Any,<:SymSemiseparableMatrix}, x::AbstractVecOrMat)
    y  = similar(x)
	L  = cholesky(K.parent)
	y .= L'\(L\x)
    return y
end
#==========================================================================================
                        Defining the Linear Algebra
==========================================================================================#
"""
    ss_mul_mat!(Y, U, V, X)

Computes the matrix-matrix product `Y = (tril(U'*V) + triu(V'*U,1))*X` in linear complexity.
"""
function ss_mul_mat!(Y, U, V, X)
    p     = size(U,1)
    n_rhs = size(X,2)
    Vbar  = zeros(p,n_rhs)
    Ubar  = U*X
    for (u,v,y,x) in zip(eachcol(U),eachcol(V),eachrow(Y),eachrow(X))
        add_inner_minus!(Ubar,u,x)    # Ubar -= u .* x
        add_inner_plus!(Vbar,v,x)     # Vbar += v .* x
        add_Y!(y,Ubar,Vbar,u,v)       # Y[i,:] = Vbar'*u + Ubar'*v
    end
end

"""
    ss_create_v(U, W)

Using `L = tril(U*W')`, compute `V` such that `LL = tril(UV') + triu(V'*U,1)`.
"""
function ss_create_v(U, W)
    m,n = size(U)
    V   = zeros(n,m)
    P   = zeros(m,m)
    @inbounds for i = 1:n
        tmpW = @view W[:,i]
        tmpU = @view U[:,i]
        tmpV = tmpW*(tmpU'*tmpW)
        V[i,:] = tmpV + P*tmpU
        add_outer_product!(P,tmpW)
    end
    return V
end


text:
struct ContinuousFunction{D,T<:Function}
    f::T
end

ContinuousFunction{D}(f::T) where {D,T} = ContinuousFunction{D,T}(f)

(func::ContinuousFunction)(args...; kwargs...) = func.f(args...; kwargs...)

### OneDimFunction

OneDimFunction(f::T) where {T} = ContinuousFunction{1}(f::T)

@recipe function f(func::ContinuousFunction{1,F}, start, stop; length=100) where {F}
    label --> :none

    x = range(start, stop; length=length)
    y = map(xi -> func(xi), x)
    x, y
end

### TwoDimFunction

TwoDimFunction(f::T) where {T} = ContinuousFunction{2}(f::T)

@recipe function f(func::ContinuousFunction{2}, xstart, xstop, ystart, ystop; length=100)
    label --> :none

    xbins = range(xstart, xstop; length=length)
    ybins = range(ystart, ystop; length=length)
    xygrid = [[xy...] for xy in Iterators.product(xbins, ybins)]
    zgrid = func(hcat(xygrid[:]...))
    zgrid = reshape(zgrid, size(xygrid)...)
    xbins, ybins, zgrid
end

### Distributions

using Distributions: ContinuousUnivariateDistribution, ContinuousMultivariateDistribution, pdf

function default_range(dist::ContinuousUnivariateDistribution, n::Int=3)
    μ, σ = mean(dist), std(dist)
    return (μ - n * σ, μ + n * σ)
end

@recipe function f(dist::T, args...) where {T<:ContinuousUnivariateDistribution}
    start, stop = isempty(args) ? default_range(dist) : args
    ContinuousFunction{1}(x -> pdf(dist, x)), start, stop
end

function default_range(dist::ContinuousMultivariateDistribution, n::Int=3)
    μ, σ = mean(dist), sqrt.(var(dist))
    start, stop = μ - n * σ, μ + n * σ
    return (start[1], stop[1], start[2], stop[2])
end

@recipe function f(dist::T, args...) where {T<:ContinuousMultivariateDistribution}
    xstart, xstop, ystart, ystop = isempty(args) ? default_range(dist) : args
    ContinuousFunction{2}(x -> pdf(dist, x)), xstart, xstop, ystart, ystop
end


text:
using MultipleScattering

function run_time_response_single_particle(;
        ωs = LinRange(0.0,50.0,2000),
        particle_x = 100.0
    )

    # Define the particle
    radius = 1.0
    particle_medium = Acoustic(2; ρ = 10.0, c = 0.5)

    # Define positions and radii of particle
    particles = [Particle(particle_medium, Sphere([particle_x,0.0],radius))]

    # Define source wave
    host_medium = Acoustic(1.0, 1.0, 2)
    source =  plane_source(host_medium; position = [0.0,0.0], direction = [1.0,0.0])

    # Simulate a single particle in frequency space
    freq_result = run(particles, source, [0.0,0.0], ωs; min_basis_order = 1, basis_order = 14)

    # Convert the frequency simulation into a time simulation
    time_result = frequency_to_time(freq_result; impulse = GaussianImpulse(1.0; σ = 10/maximum(ωs)^2))

    return freq_result, time_result
end

function plot_time_response_single_particle()
    freq_result, time_result = run_time_response_single_particle()

    plot(
        plot(freq_result),
        plot(time_result)
    )
end


text:
#####################
# Compilers to ProbabilisticCircuits data structures starting from already parsed line objects
#####################

# reuse some internal infrastructure of LogicCircuits' LoadSave module
using LogicCircuits.LoadSave: CircuitFormatLines, CircuitFormatLine, lnconstant,
VtreeFormatLines, CircuitHeaderLine, UnweightedLiteralLine, WeightedLiteralLine,
DecisionLine, LCElement, BiasLine, WeightedNamedConstantLine, PSDDElement,
CircuitCommentLine, ID, compile_smooth_struct_logical_m, compile_smooth_logical_m

"""
Compile lines into a probabilistic circuit
"""
function compile_prob(lines::CircuitFormatLines)::ProbCircuit
    # first compile a logic circuit
    logic_circuit, id2lognode = compile_smooth_logical_m(lines)
    decorate_prob(lines, logic_circuit, id2lognode)
end

"""
Compile lines into a logistic circuit.
"""
function compile_logistic(lines::CircuitFormatLines, classes::Int)::LogisticCircuit
    # first compile a logic circuit
    logic_circuit, id2lognode = compile_smooth_logical_m(lines)
    decorate_logistic(lines, logic_circuit, classes, id2lognode)
end

"""
Compile circuit and vtree lines into a structured probabilistic circuit (one whose logic circuit origin is structured).
"""
function compile_struct_prob(circuit_lines::CircuitFormatLines, vtree_lines::VtreeFormatLines)
    logic_circuit, vtree, id2lognode, id2vtree = compile_smooth_struct_logical_m(circuit_lines, vtree_lines)
    prob_circuit = decorate_prob(circuit_lines, logic_circuit, id2lognode)
    return prob_circuit, vtree
end

function decorate_prob(lines::CircuitFormatLines, logic_circuit::LogicCircuit, id2lognode::Dict{ID,<:LogicCircuit})::ProbCircuit
    # set up cache mapping logic circuit nodes to their probabilistic decorator

    prob_circuit = ProbCircuit(logic_circuit)
    lognode2probnode = Dict{LogicCircuit, ProbCircuit}()

    prob_lin = linearize(prob_circuit) # TODO better implementation
    logic_lin = linearize(logic_circuit)

    foreach(i -> lognode2probnode[logic_lin[i]] = prob_lin[i], 1 : num_nodes(logic_circuit))

    # map from line node ids to probabilistic circuit nodes
    id2probnode(id) = lognode2probnode[id2lognode[id]]

    root = nothing

    # go through lines again and update the probabilistic circuit node parameters

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln into probabilistic circuit is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine,UnweightedLiteralLine})
        # do nothing
    end
    function compile(ln::WeightedNamedConstantLine)
        @assert lnconstant(ln) == true
        root = id2probnode(ln.node_id)
        root.log_probs .= [ln.weight, log1p(-exp(ln.weight))]
    end
    function compile(ln::DecisionLine{<:PSDDElement})
        root = id2probnode(ln.node_id)
        root.log_probs .= [x.weight for x in ln.elements]
    end

    foreach(compile, lines)

    root
end

function decorate_logistic(lines::CircuitFormatLines, logic_circuit::LogicCircuit,
                            classes::Int, id2lognode::Dict{ID,<:LogicCircuit})::LogisticCircuit

    # set up cache mapping logic circuit nodes to their logistic decorator
    logistic_circuit = LogisticCircuit(logic_circuit, classes)
    log2logistic = Dict{LogicCircuit, LogisticCircuit}()
    logistic_lin = linearize(logistic_circuit)
    logic_lin = linearize(logic_circuit)

    foreach(i -> log2logistic[logic_lin[i]] = logistic_lin[i], 1 : length(logic_lin))
    id2logisticnode(id) = log2logistic[id2lognode[id]]

    root = nothing
    # go through lines again and update the probabilistic circuit node parameters

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln into logistic circuit is not supported")
    end

    function compile(::Union{CircuitHeaderLine,CircuitCommentLine,UnweightedLiteralLine})
        # do nothing
    end

    function compile(ln::WeightedLiteralLine)
        root = id2logisticnode(ln.node_id)::Logistic⋁Node
        root.thetas[1, :] .= ln.weights
    end

    function compile(ln::DecisionLine{<:LCElement})
        root = id2logisticnode(ln.node_id)::Logistic⋁Node
        for (ind, elem) in enumerate(ln.elements)
            root.thetas[ind, :] .= elem.weights
        end
    end

    function compile(ln::BiasLine)
        root = id2logisticnode(ln.node_id)::Logistic⋁Node
        # @assert length(root.thetas) == 1
        root.thetas[1,:] .= ln.weights
    end

    foreach(compile, lines)

    root
end


text:
# Source: https://www.johndcook.com/blog/2020/02/08/arenstorf-orbit/

using Revise
using LinearAlgebra
using Plots

includet("src/InternalArray.jl")
includet("src/FixedStepMethods.jl")
includet("src/MixedStepMethods.jl")
IArray = InternalArray

mmoon = 0.012277471
mearth = 1-mmoon
p0 = [0.994, 0]
v0 = [0, -2.001585106]
target_err = 1e-4

function derivative(t, x)
    if ndims(x) == 2
        r = zeros(size(x))
        for i in 1:size(x, 2)
            r[:,i] = derivative(t, x[:,i])
        end
        return r
    end
    D1 = ((x[1]+mmoon)*(x[1]+mmoon) + x[2]*x[2])
    D1 *= sqrt(D1)
    D2 = ((x[1]-mearth)*(x[1]-mearth) + x[2]*x[2])
    D2 *= sqrt(D2)
    return [
        x[3],
        x[4],
        x[1] + 2*x[4] - mearth*(x[1]+mmoon)/D1 - mmoon*(x[1]-mearth)/D2,
        x[2] - 2*x[3] - mearth*x[2]/D1 - mmoon*x[2]/D2
    ]
end

t, x = mixed_step_ode_solver(
    derivative,
    0,
    [p0; v0],
    (t, x) ->
        last(t) > 5 && norm(p0 .- last(x)[1:2]) < target_err*last(t) || last(t) > 50,
    min_step = 1e-7,
    max_step = 0.1,
    target_err = target_err,
    method = mod_euler_step,
);

plot(
    IArray(x, 1),
    IArray(x, 2),
    label = "satellite",
    title = "Arenstorf orbit",
    aspect_ratio = :equal,
    framestyle = :box,
    xlims = (-1.4, 1.2),
    legend = :bottomright
)
plot!(cos, sin, 0, 2π, label = "moon", color = :black, style = :dot)
scatter!((0, 0), label = "", color = :black, markersize = 1.5)
# savefig("media/arenstorf_orbit.png")

using IJulia

px = [x[i][1] for i in 1:length(x)]
py = [x[i][2] for i in 1:length(x)]
anim = Animation()
i = 1
for j = 0:0.1:last(t)
    while i < length(t) && t[i] < j
        i += 1
    end
    plot(
        px[1:i],
        py[1:i],
        label = "satellite",
        title = "Arenstorf orbit",
        aspect_ratio = :equal,
        framestyle = :box,
        xlims = (-1.4, 1.2),
        ylims = (-1.3, 1.4),
        legend = :bottomright,
    )
    inds = range(0, 2*π*j/last(t), length = 150)
    plot!(cos.(inds), sin.(inds), label = "moon", color = :black, style = :dot)
    scatter!((0, 0), label = "", color = :black, markersize = 1.5)
    frame(anim)
end
gif(anim, "media/arenstorf_orbit.gif", fps = 15)

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
type DiscretizedVanillaOption <: DiscretizedAsset
  stoppingTimes::Vector{Float64}
  args::VanillaOptionArgs
  common::DiscretizedAssetCommon
end

function DiscretizedVanillaOption(vanillaOption::VanillaOption, process::StochasticProcess, grid::TimeGrid)
  stoppingTimes = zeros(length(vanillaOption.exercise.dates))
  for i in eachindex(stoppingTimes)
    stoppingTimes[i] = get_time(process, vanillaOption.exercise.dates[i])

    if ~is_empty(grid)
      # stuff
      stoppingTimes[i] = closest_time(grid, stoppingTimes[i])
    end
  end

  return DiscretizedVanillaOption(stoppingTimes, VanillaOptionArgs(vanillaOption.payoff, vanillaOption.exercise), DiscretizedAssetCommon())
end

function reset!(dvo::DiscretizedVanillaOption, sz::Int)
  dvo.common.values = zeros(sz)
  adjust_values!(dvo)

  return dvo
end

post_adjust_values_impl!(dvo::DiscretizedVanillaOption) = post_adjust_values_impl!(dvo, dvo.args.exercise, dvo.common.time)

function post_adjust_values_impl!(dvo::DiscretizedVanillaOption, ::AmericanExercise, time_now::Float64)
  if time_now <= dvo.stoppingTimes[2] && time_now >= dvo.stoppingTimes[1]
    apply_specific_condition!(dvo)
  end
  return dvo
end

function post_adjust_values_impl!(dvo::DiscretizedVanillaOption, ::EuropeanExercise, time_now::Float64)
  if is_on_time(dvo, dvo.stoppingTimes[1])
    apply_specific_condition!(dvo)
  end
  return dvo
end

function post_adjust_values_impl!(dvo::DiscretizedVanillaOption, ::BermudanExercise, time_now::Float64)
  for i in eachindex(dvo.stoppingTimes)
    if is_on_time(dvo, dvo.stoppingTimes[i])
      apply_specific_condition!(dvo)
    end
  end
  return dvo
end

function apply_specific_condition!(dvo::DiscretizedVanillaOption)
  grid = get_grid(dvo.common.method, dvo.common.time)
  for j in eachindex(dvo.common.values)
    dvo.common.values[j] = max(dvo.common.values[j], dvo.args.payoff(grid[j]))
  end

  return dvo
end


text:
#Controlled HF noise plot (test data reduction)
using InspectDR
using Colors

#==Input constants
===============================================================================#
dfltline = line(color=RGB24(0, 0, 1), width=1, style=:solid)

#==Input data
===============================================================================#
function gennoise(xv, yv; step=.01, npts=20, lo=0, hi=1)
	xmax = xv[end]
	xnew = xmax+collect(step:step:(npts*step))
	ynew = similar(xnew)
	v = Float64[lo, hi]
	for i in 1:length(ynew)
		ynew[i] = v[(i&0x1)+1]
	end
	append!(xv, xnew)
	append!(yv, ynew)
end

function addpts(ynew, xv, yv; step=1.0)
	xmax = xv[end]
	npts = length(ynew)
	xnew = xmax+collect(step:step:(npts*step))
	append!(xv, xnew)
	append!(yv, ynew)
end

#Generate dataset:
	x = Float64[0]
	y = Float64[0]
	addpts(Float64[1.1,-.1,], x,y) #Force maximum y-extents
	addpts(Float64[1,1,1,0,0], x,y)
	gennoise(x,y)
	addpts(Float64[0,1,1,1,0,0], x,y)
	gennoise(x,y)
	addpts(Float64[1,1,1,0,0], x,y)
	gennoise(x,y)
	addpts(Float64[1,1,1,0,0], x,y)

#==Generate plot
===============================================================================#
plot = InspectDR.Plot2D()
wfrm = add(plot, x, y+1)
wfrm.line = dfltline
plot.xres=1000 #Force resolution

a = plot.annotation
a.title = "Sample Plot (λ)"
a.xlabel = "Time (s)"
a.ylabels = ["Signal Voltage (V)"]

gplot = display(InspectDR.GtkDisplay(), plot)

:Test_Complete


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMDynamics.jl/blob/master/LICENSE

using FEMDynamics
using FEMBase.Test

@testset "LinearImplicit" begin
    include("test_dropping_element.jl")
end


text:
using Stheno: FiniteGP, block_diagonal
using Distributions: MvNormal

_rng() = MersenneTwister(123456)

function generate_noise_matrix(rng::AbstractRNG, N::Int)
    A = randn(rng, N, N)
    return Symmetric(A * A' + I)
end

@testset "finite_gp" begin
    @testset "statistics" begin
        rng, N, N′ = MersenneTwister(123456), 1, 9
        x, x′, Σy, Σy′ = randn(rng, N), randn(rng, N′), zeros(N, N), zeros(N′, N′)
        f = wrap(GP(sin, SEKernel()), GPC())
        fx, fx′ = FiniteGP(f, x, Σy), FiniteGP(f, x′, Σy′)

        @test mean(fx) == mean(f, x)
        @test cov(fx) == cov(f, x)
        @test cov(fx, fx′) == cov(f, x, x′)
        @test mean.(marginals(fx)) == mean(f(x))
        @test var.(marginals(fx)) == var(f, x)
        @test std.(marginals(fx)) == sqrt.(var(f, x))
    end
    @testset "rand (deterministic)" begin
        rng, N, D = MersenneTwister(123456), 10, 2
        X, x, Σy = ColVecs(randn(rng, D, N)), randn(rng, N), zeros(N, N)
        Σy = generate_noise_matrix(rng, N)
        fX = FiniteGP(wrap(GP(1, SEKernel()), GPC()), X, Σy)
        fx = FiniteGP(wrap(GP(1, SEKernel()), GPC()), x, Σy)

        # Check that single-GP samples have the correct dimensions.
        @test length(rand(rng, fX)) == length(X)
        @test size(rand(rng, fX, 10)) == (length(X), 10)

        @test length(rand(rng, fx)) == length(x)
        @test size(rand(rng, fx, 10)) == (length(x), 10)
    end
    @testset "rand (statistical)" begin
        rng, N, D, μ0, S = MersenneTwister(123456), 10, 2, 1, 100_000
        X, Σy = ColVecs(randn(rng, D, N)), 1e-12
        f = FiniteGP(wrap(GP(1, SEKernel()), GPC()), X, Σy)

        # Check mean + covariance estimates approximately converge for single-GP sampling.
        f̂ = rand(rng, f, S)
        @test maximum(abs.(mean(f̂; dims=2) - mean(f))) < 1e-2

        Σ′ = (f̂ .- mean(f)) * (f̂ .- mean(f))' ./ S
        @test mean(abs.(Σ′ - cov(f))) < 1e-2
    end
    @testset "rand (gradients)" begin
        rng, N, S = MersenneTwister(123456), 10, 3
        x = collect(range(-3.0, stop=3.0, length=N))
        Σy = 1e-12

        # Check that the gradient w.r.t. the samples is correct (single-sample).
        adjoint_test(
            x->rand(
                MersenneTwister(123456),
                FiniteGP(wrap(GP(sin, SEKernel()), GPC()), x, Σy),
            ),
            randn(rng, N),
            x;
            atol=1e-9, rtol=1e-9,
        )

        # Check that the gradient w.r.t. the samples is correct (multisample).
        adjoint_test(
            x->rand(
                MersenneTwister(123456),
                FiniteGP(wrap(GP(sin, SEKernel()), GPC()), x, Σy),
                S,
            ),
            randn(rng, N, S),
            x;
            atol=1e-9, rtol=1e-9,
        )
    end
    @testset "logpdf / elbo / dtc" begin
        rng, N, S, σ, gpc = MersenneTwister(123456), 10, 11, 1e-1, GPC()
        x = collect(range(-3.0, stop=3.0, length=N))
        f = wrap(GP(1, SEKernel()), gpc)
        fx, y = FiniteGP(f, x, 0), FiniteGP(f, x, σ^2)
        ŷ = rand(rng, y)

        # Check that logpdf returns the correct type and roughly agrees with Distributions.
        @test logpdf(y, ŷ) isa Real
        @test logpdf(y, ŷ) ≈ logpdf(MvNormal(Vector(mean(y)), cov(y)), ŷ)

        # Check that multi-sample logpdf returns the correct type and is consistent with
        # single-sample logpdf
        Ŷ = rand(rng, y, S)
        @test logpdf(y, Ŷ) isa Vector{Float64}
        @test logpdf(y, Ŷ) ≈ [logpdf(y, Ŷ[:, n]) for n in 1:S]

        # Check gradient of logpdf at mean is zero for `f`.
        adjoint_test(ŷ->logpdf(fx, ŷ), 1, ones(size(ŷ)))
        lp, back = Zygote.pullback(ŷ->logpdf(fx, ŷ), ones(size(ŷ)))
        @test back(randn(rng))[1] == zeros(size(ŷ))

        # Check that gradient of logpdf at mean is zero for `y`.
        adjoint_test(ŷ->logpdf(y, ŷ), 1, ones(size(ŷ)))
        lp, back = Zygote.pullback(ŷ->logpdf(y, ŷ), ones(size(ŷ)))
        @test back(randn(rng))[1] == zeros(size(ŷ))

        # Check that gradient w.r.t. inputs is approximately correct for `f`.
        x, l̄ = randn(rng, N), randn(rng)
        adjoint_test(
            x->logpdf(FiniteGP(f, x, 1e-3), ones(size(x))),
            l̄, collect(x);
            atol=1e-8, rtol=1e-8,
        )
        adjoint_test(
            x->sum(logpdf(FiniteGP(f, x, 1e-3), ones(size(Ŷ)))),
            l̄, collect(x);
            atol=1e-8, rtol=1e-8,
        )

        # Check that the gradient w.r.t. the noise is approximately correct for `f`.
        σ_ = randn(rng)
        adjoint_test((σ_, ŷ)->logpdf(FiniteGP(f, x, exp(σ_)), ŷ), l̄, σ_, ŷ)
        adjoint_test((σ_, Ŷ)->sum(logpdf(FiniteGP(f, x, exp(σ_)), Ŷ)), l̄, σ_, Ŷ)

        # Check that the gradient w.r.t. a scaling of the GP works.
        adjoint_test(
            α->logpdf(FiniteGP(α * f, x, 1e-1), ŷ), l̄, randn(rng);
            atol=1e-8, rtol=1e-8,
        )
        adjoint_test(
            α->sum(logpdf(FiniteGP(α * f, x, 1e-1), Ŷ)), l̄, randn(rng);
            atol=1e-8, rtol=1e-8,
        )

        # Ensure that the elbo is close to the logpdf when appropriate.
        @test elbo(y, ŷ, fx) isa Real
        @test elbo(y, ŷ, fx) ≈ logpdf(y, ŷ)
        @test elbo(y, ŷ, y) < logpdf(y, ŷ)
        @test elbo(y, ŷ, FiniteGP(f, x, 2 * σ^2)) < elbo(y, ŷ, y)

        # Check adjoint w.r.t. elbo is correct.
        adjoint_test(
            (x, ŷ, σ)->elbo(FiniteGP(f, x, σ^2), ŷ, FiniteGP(f, x, 0)),
            randn(rng), x, ŷ, σ;
            atol=1e-6, rtol=1e-6,
        )

        # Ensure that the dtc is close to the logpdf when appropriate.
        @test dtc(y, ŷ, fx) isa Real
        @test dtc(y, ŷ, fx) ≈ logpdf(y, ŷ)

        # Check adjoint w.r.t. dtc is correct.
        adjoint_test(
            (x, ŷ, σ)->dtc(FiniteGP(f, x, σ^2), ŷ, FiniteGP(f, x, 0)),
            randn(rng), x, ŷ, σ;
            atol=1e-6, rtol=1e-6,
        )
    end
    @testset "Type Stability - $T" for T in [Float64, Float32]
        rng = MersenneTwister(123456)
        x = randn(rng, T, 123)
        z = randn(rng, T, 13)
        f = wrap(GP(T(0), SEKernel()), GPC())

        fx = f(x, T(0.1))
        u = f(z, T(1e-4))

        y = rand(rng, fx)
        @test y isa Vector{T}
        @test logpdf(fx, y) isa T
        @test elbo(fx, y, u) isa T
    end
end