text:
"""
Driver to run TRON on GPU. This should be called from a kernel.
"""
@inline function tron_kernel(n::Int, shift::Int, max_feval::Int, max_minor::Int, gtol::Float64, scale::Float64, use_polar::Bool,
                     x::CuDeviceArray{Float64,1}, xl::CuDeviceArray{Float64,1},
                     xu::CuDeviceArray{Float64,1},
                     param::CuDeviceArray{Float64,2},
                     YffR::Float64, YffI::Float64,
                     YftR::Float64, YftI::Float64,
                     YttR::Float64, YttI::Float64,
                     YtfR::Float64, YtfI::Float64)
    tx = threadIdx().x

    g = @cuDynamicSharedMem(Float64, n, (3*n)*sizeof(Float64))
    xc = @cuDynamicSharedMem(Float64, n, (4*n)*sizeof(Float64))
    s = @cuDynamicSharedMem(Float64, n, (5*n)*sizeof(Float64))
    wa = @cuDynamicSharedMem(Float64, n, (6*n)*sizeof(Float64))
    wa1 = @cuDynamicSharedMem(Float64, n, (7*n)*sizeof(Float64))
    wa2 = @cuDynamicSharedMem(Float64, n, (8*n)*sizeof(Float64))
    wa3 = @cuDynamicSharedMem(Float64, n, (9*n)*sizeof(Float64))
    wa4 = @cuDynamicSharedMem(Float64, n, (10*n)*sizeof(Float64))
    wa5 = @cuDynamicSharedMem(Float64, n, (11*n)*sizeof(Float64))
    gfree = @cuDynamicSharedMem(Float64, n, (12*n)*sizeof(Float64))
    dsave = @cuDynamicSharedMem(Float64, n, (13*n)*sizeof(Float64))
    indfree = @cuDynamicSharedMem(Int, n, (14*n)*sizeof(Float64))
    iwa = @cuDynamicSharedMem(Int, 2*n, n*sizeof(Int) + (14*n)*sizeof(Float64))
    isave = @cuDynamicSharedMem(Int, n, (3*n)*sizeof(Int) + (14*n)*sizeof(Float64))

    A = @cuDynamicSharedMem(Float64, (n,n), (14*n)*sizeof(Float64)+(4*n)*sizeof(Int))
    B = @cuDynamicSharedMem(Float64, (n,n), (14*n+n^2)*sizeof(Float64)+(4*n)*sizeof(Int))
    L = @cuDynamicSharedMem(Float64, (n,n), (14*n+2*n^2)*sizeof(Float64)+(4*n)*sizeof(Int))

    if tx <= n
        @inbounds for j=1:n
            A[tx,j] = 0.0
            B[tx,j] = 0.0
            L[tx,j] = 0.0
        end
    end
    CUDA.sync_threads()

    task = 0
    status = 0

    delta = 0.0
    fatol = 0.0
    frtol = 1e-12
    fmin = -1e32
    cgtol = 0.1
    cg_itermax = n

    f = 0.0
    nfev = 0
    ngev = 0
    nhev = 0
    minor_iter = 0
    search = true

    while search

        # [0|1]: Evaluate function.

        if task == 0 || task == 1
            if use_polar
                f = eval_f_polar_kernel(n, shift, scale, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            else
                f = eval_f_kernel(n, scale, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            end
            nfev += 1
            if nfev >= max_feval
                search = false
            end

        end

        # [2] G or H: Evaluate gradient and Hessian.

        if task == 0 || task == 2
            if use_polar
                eval_grad_f_polar_kernel(n, shift, scale, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
                eval_h_polar_kernel(n, shift, scale, x, A, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            else
                eval_grad_f_kernel(n, scale, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
                eval_h_kernel(n, scale, x, A, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            end
            ngev += 1
            nhev += 1
            minor_iter += 1
        end

        # Initialize the trust region bound.

        if task == 0
            gnorm0 = dnrm2(n, g, 1)
            delta = gnorm0
        end

        # Call Tron.

        if search
            delta, task = dtron(n, x, xl, xu, f, g, A, frtol, fatol, fmin, cgtol,
                                cg_itermax, delta, task, B, L, xc, s, indfree, gfree,
                                isave, dsave, wa, iwa, wa1, wa2, wa3, wa4, wa5)
        end

        # [3] NEWX: a new point was computed.

        if task == 3
            gnorm_inf = dgpnorm(n, x, xl, xu, g)

            if gnorm_inf <= gtol
                task = 4
            end

            if minor_iter >= max_minor
                status = 1
                search = false
            end
        end

        # [4] CONV: convergence was achieved.
        # [10] : warning fval is less than fmin

        if task == 4 || task == 10
            search = false
        end
    end

    CUDA.sync_threads()

    return status, minor_iter
end

@inline function tron_qp_kernel(n::Int, max_feval::Int, max_minor::Int, gtol::Float64, scale::Float64,
    x::CuDeviceArray{Float64,1}, xl::CuDeviceArray{Float64,1}, xu::CuDeviceArray{Float64,1},
    A::CuDeviceArray{Float64,2}, c::CuDeviceArray{Float64,1})

    tx = threadIdx().x
    I = blockIdx().x

    g = @cuDynamicSharedMem(Float64, n, (3*n)*sizeof(Float64))
    xc = @cuDynamicSharedMem(Float64, n, (4*n)*sizeof(Float64))
    s = @cuDynamicSharedMem(Float64, n, (5*n)*sizeof(Float64))
    wa = @cuDynamicSharedMem(Float64, n, (6*n)*sizeof(Float64))
    wa1 = @cuDynamicSharedMem(Float64, n, (7*n)*sizeof(Float64))
    wa2 = @cuDynamicSharedMem(Float64, n, (8*n)*sizeof(Float64))
    wa3 = @cuDynamicSharedMem(Float64, n, (9*n)*sizeof(Float64))
    wa4 = @cuDynamicSharedMem(Float64, n, (10*n)*sizeof(Float64))
    wa5 = @cuDynamicSharedMem(Float64, n, (11*n)*sizeof(Float64))
    gfree = @cuDynamicSharedMem(Float64, n, (12*n)*sizeof(Float64))
    dsave = @cuDynamicSharedMem(Float64, 3, (13*n)*sizeof(Float64))
    indfree = @cuDynamicSharedMem(Int, n, (13*n+3)*sizeof(Float64))
    iwa = @cuDynamicSharedMem(Int, 2*n, n*sizeof(Int) + (13*n+3)*sizeof(Float64))
    isave = @cuDynamicSharedMem(Int, 3, (3*n)*sizeof(Int) + (13*n+3)*sizeof(Float64))

    B = @cuDynamicSharedMem(Float64, (n,n), (13*n+3+n^2)*sizeof(Float64)+(3*n+3)*sizeof(Int))
    L = @cuDynamicSharedMem(Float64, (n,n), (13*n+3+2*n^2)*sizeof(Float64)+(3*n+3)*sizeof(Int))

    if tx <= n
        @inbounds begin
            for j=1:n
                B[tx,j] = 0.0
                L[tx,j] = 0.0
            end
        end
    end
    CUDA.sync_threads()

    task = 0
    status = 0

    delta = 0.0
    fatol = 0.0
    frtol = 1e-12
    fmin = -1e32
    cgtol = 0.1
    cg_itermax = n

    f = 0.0
    nfev = 0
    ngev = 0
    nhev = 0
    minor_iter = 0
    search = true

    while search

        # [0|1]: Evaluate function.

        if task == 0 || task == 1
            f = eval_qp_f_kernel(n, x, A, c)
            nfev += 1
            if nfev >= max_feval
                search = false
            end
        end

        # [2] G or H: Evaluate gradient and Hessian.

        if task == 0 || task == 2
            eval_qp_grad_f_kernel(n, x, g, A, c)
            # We do not have to evaluate Hessian since A does not change.
            ngev += 1
            nhev += 1
            minor_iter += 1
        end

        # Initialize the trust region bound.

        if task == 0
            gnorm0 = dnrm2(n, g, 1)
            delta = gnorm0
        end

        # Call Tron.

        if search
            delta, task = dtron(n, x, xl, xu, f, g, A, frtol, fatol, fmin, cgtol,
                                cg_itermax, delta, task, B, L, xc, s, indfree, gfree,
                                isave, dsave, wa, iwa, wa1, wa2, wa3, wa4, wa5)
        end

        # [3] NEWX: a new point was computed.

        if task == 3
            gnorm_inf = dgpnorm(n, x, xl, xu, g)

            if gnorm_inf <= gtol
                task = 4
            end

            if minor_iter >= max_minor
                status = 1
                search = false
            end
        end

        # [4] CONV: convergence was achieved.
        # [10] : warning fval is less than fmin

        if task == 4 || task == 10
            search = false
        end
    end

    CUDA.sync_threads()

    return status, minor_iter
end


text:
############ Symbolic intervals that store intermediate bounds #################

# must have fields sym, lbs, ubs
abstract type AbstractSymbolicIntervalBounds{F<:AbstractPolytope} end

struct SymbolicIntervalBounds{F} <: AbstractSymbolicIntervalBounds{F}
    sym::SymbolicInterval{F}
    lbs::Vector{Vector{Float64}}
    ubs::Vector{Vector{Float64}}
end

domain(s::A) where A<:AbstractSymbolicIntervalBounds = domain(s.sym)
LazySets.radius(s::A) where A<:AbstractSymbolicIntervalBounds = LazySets.radius(domain(s))

function init_symbolic_interval_bounds(network::NetworkNegPosIdx, input::AbstractHyperrectangle)
    sym = init_deep_poly_symbolic_interval(input)
    lbs = [fill(-Inf, n_nodes(l)) for l in network.layers]
    ubs = [fill( Inf, n_nodes(l)) for l in network.layers]
    return SymbolicIntervalBounds(sym, lbs, ubs)
end

function init_symbolic_interval_bounds(s::SymbolicIntervalBounds, input::AbstractHyperrectangle)
    sym = init_deep_poly_symbolic_interval(input)
    lbs = [copy(lb) for lb in s.lbs]
    ubs = [copy(ub) for ub in s.ubs]
    return SymbolicIntervalBounds(sym, lbs, ubs)
end

function bounds(input::SymbolicInterval, lbs::Vector{Float64}, ubs::Vector{Float64})
    l̂, û = bounds(input)
    l = max.(l̂, lbs)
    u = min.(û, ubs)
    return l, u
end

# TODO: maybe just calculate upper bound instead of both lower and upper bound and then just returning upper bound
function LazySets.ρ(d::AbstractArray{T,1} where T, sym::A where A<:SymbolicIntervalBounds)
    d_Low, d_Up = interval_map(reshape(d, 1, length(d)), sym.sym.Low, sym.sym.Up)
    sym_prime = SymbolicInterval(d_Low, d_Up, domain(sym))
    lo, up = bounds(sym_prime)
    # the result should only be one-dimensional
    return up[1]
end

"""
Calculates an input within the hyperrectangle [lbs, ubs] that maximizes the linear
symbolic equation sym_eq.
"""
function maximizer(sym_eq, lbs, ubs)
    if size(sym_eq, 1) == 1
        W⁺ = (sym_eq[1:end-1] .> 0)
        W⁻ = (sym_eq[1:end-1] .< 0)
        maximizer = W⁺ .* ubs + W⁻ .* lbs
    else
        W⁺ = (sym_eq[:, 1:end-1] .> 0)
        W⁻ = (sym_eq[:, 1:end-1] .< 0)
        maximizer =( W⁺' .* ubs + W⁻' .* lbs)'
    end

    return maximizer
end

function minimizer(sym_eq, lbs, ubs)
    # minimizer of eq is just maximizer of -eq
    return maximizer(-sym_eq, lbs, ubs)
end

function maximizer(s::S) where S <: AbstractSymbolicIntervalBounds
    maximizer(s.sym.Up, low(domain(s)), high(domain(s)))
end

function minimizer(s::S) where S <: AbstractSymbolicIntervalBounds
    minimizer(s.sym.Low, low(domain(s)), high(domain(s)))
end

##### Splitting

function Base.split(H::AbstractHyperrectangle, index::Int)
    lbs, ubs = low(H), high(H)
    split_val = 0.5 * (lbs[index] + ubs[index])

    high1 = copy(ubs)
    high1[index] = split_val
    low2 = copy(lbs)
    low2[index] = split_val

    domain1 = Hyperrectangle(low=lbs, high=high1)
    domain2 = Hyperrectangle(low=low2, high=ubs)

    return domain1, domain2
end

function split_symbolic_interval_bounds(s::SymbolicIntervalBounds{<:Hyperrectangle}, index::Int)
    domain1, domain2 = split(domain(s), index)

    s1 = init_symbolic_interval_bounds(s, domain1)
    s2 = init_symbolic_interval_bounds(s, domain2)

    return [s1, s2]
end

function split_largest_interval(s::SymbolicIntervalBounds{<:Hyperrectangle})
    largest_dimension = argmax(high(domain(s)) - low(domain(s)))
    return split_symbolic_interval_bounds(s, largest_dimension)
end

"""
function split_multiple_times(cell::SymbolicIntervalBounds, n; split=split_largest_interval)
    q = Queue{SymbolicIntervalBounds}()
    enqueue!(q, cell)
    for i in 1:n
        new_cells = split(dequeue!(q))
        enqueue!(q, new_cells[1])
        enqueue!(q, new_cells[2])
    end
    return q
end
"""

ext: jl


text:
"""Filtration of an abstract complex

We call this sequence of complexes the **filtration** of `f` and
think of it as a construction by adding chunks of simplices at a time `t::FI`.

∅ = K₀ ⊆ K₁ ⊆ . . . ⊆ Kₙ = K.
"""
mutable struct Filtration{C<:AbstractComplex, FI<:AbstractFloat}
    # underlying abstract cell complex
    complex::C
    # total order of simplices as array of (dim, simplex id, filtration value)
    total::Vector{Tuple{Int,UInt,FI}}
    divisions::Number
end

order(flt::Filtration) = flt.total
complex(flt::Filtration) = flt.complex
show(io::IO, flt::Filtration{C,FI}) where {C <: AbstractComplex, FI} = print(io, "Filtration($(complex(flt)), $FI)")
eltype(flt::Filtration{C,FI}) where {C <: AbstractComplex, FI} = C
valtype(flt::Filtration{C,FI}) where {C <: AbstractComplex, FI} = FI
length(flt::Filtration) = isinf(flt.divisions) ? length(unique(e->e[3], order(flt))) : flt.divisions

#
# Constructors
#
Filtration(::Type{C}, ::Type{FI}) where {C <: AbstractComplex, FI} =
    Filtration(C(), Vector{Tuple{Int,UInt,FI}}(), Inf)
Filtration(::Type{C}) where {C <: AbstractComplex} = Filtration(C, Float64)
similar(flt::Filtration{C,FI}) where {C <: AbstractComplex, FI} = Filtration(C, FI)

"""Construct filtration from a cell complex using the order of their appearence in the complex"""
function filtration(cplx::AbstractComplex)
    idx = Vector{Tuple{Int,UInt,Float64}}()
    i = 1
    for d in 0:dim(cplx)
        for c in cells(cplx, d)
            push!(idx, (d, hash(c), i))
            i += 1
        end
    end
    return Filtration(cplx, idx, Inf)
end

"""Construct filtration from a cell complex and a complex weight function"""
function filtration(cplx::C, w::Dict{Int,Vector{FI}}; divisions::Number=Inf) where {C<:AbstractComplex, FI}
    idx = Vector{Tuple{Int,UInt,FI}}()
    for d in 0:dim(cplx)
        for c in cells(cplx, d)
            ci = position(cplx, c)
            push!(idx, (d, hash(c), w[d][ci]))
        end
    end
    sort!(idx, by=x->(x[3], x[1])) # sort by dimension & filtration value
    return Filtration(cplx, idx, divisions)
end

function push!(flt::Filtration{C,FI}, cl::AbstractCell, v::FI; recursive=false) where {C<:AbstractComplex, FI}
    cplx = complex(flt)
    @assert isa(cl, eltype(cplx)) "Complex $(cplx) does not accept $(typeof(cl))"
    cls = push!(cplx, cl, recursive=recursive)
    ord = order(flt)
    idx = length(ord) == 0 ? 1 : findlast(e->e[3]<=v, ord)
    for c in sort!(cls, by=s->dim(s))
        if idx == length(ord)
            push!(ord, (dim(c), hash(c), v))
        else
            insert!(ord, idx, (dim(c), hash(c), v))
        end
        idx += 1
    end
    return flt
end

#
# Auxiliary function
#
minimum(flt::Filtration) = order(flt)[1][3]
maximum(flt::Filtration) = order(flt)[end][3]

"""
    complex(flt, val)

Return a complex from the filtration `flt` at the filtration value `val`
"""
function complex(flt::Filtration{C,FI}, val::FI) where {C <: AbstractComplex, FI}
    cplx = complex(flt)
    res = eltype(flt)()
    for (d,si,fv) in order(flt)
        fv > val && break
        push!(res, cplx[si, d])
    end
    return res
end

"""
    boundary(flt::Filtration) -> Vector{<:AbsrtactSet}

Generate a boundary matrix from the filtration `flt` for the persistent homology calculations.
"""
function boundary(flt::Filtration; reduced=false)
    ridx = reduced ? 1 : 0
    # initialize boundary matrix
    cplx = complex(flt)
    sz = sum(size(cplx))+ridx

    # find suitable type for BM storage
    STs = [Int8, Int16, Int32, Int64]
    STi = findfirst(i->i>sz, map(typemax, STs))
    BT = STi == 1 ? BitSet : Set{STs[STi]}

    # create empty BM
    bm = map(i->BT(), 1:sz)

    # fill boundary matrix
    ord = order(flt)
    revidx = Dict((ci, d) => i for (i, (d, ci, fv)) in enumerate(ord))
    for (i, (d, ci, fv)) in enumerate(ord)
        if d > 0
            for face in faces(cplx, ci)
                push!(bm[i+ridx], revidx[(face, d-1)]+ridx)
            end
        elseif reduced
            push!(bm[i+ridx], 1)
        end
    end
    return bm
end

function sparse(∂::Vector{<:AbstractSet})
    m = length(∂)
    ret = spzeros(Int, m, m)
    for i in 1:m
        bm = ∂[i]
        for (l, j) in enumerate(bm)
            ret[j,i] = j # (-1)^((l-1)%2) # coefficients require exact order of faces in provides simplex
        end
    end
    return ret
end

"""
    similarity(flt::Filtration)

Returns a similarity matrix created from 1-subcomplex of simplicial complex and distance weights of the filtration `flt`.
"""
function similarity(flt::Filtration)
    cplx = complex(flt)
    ord = order(flt)
    C0 = map(hash, cells(cplx, 0))
    N = length(C0)
    @assert N > 0 "Complex should not be empty"
    adj = spzeros(valtype(flt),N,N)
    for c in cells(cplx, 1)
        i, j = map(h->findfirst(isequal(h), C0), vertices(c))
        idx = findfirst(v->v[1]==1 && v[2] == hash(c), ord)
        adj[i, j] = adj[j, i] = ord[idx][3]
    end
    return adj
end

#
# I/O
#
function write(io::IO, flt::Filtration)
    cplx = complex(flt)
    for (d, ci, fv) in order(flt)
        for k in indices(cplx[ci,d])
            write(io, "$k,")
        end
        write(io, "$fv\n")
    end
end

function read(io::IO, ::Type{Filtration{C,FI}}, ::Type{CL}) where {C <: AbstractComplex, FI, CL<:AbstractCell}
    flt = Filtration(C,FI)
    while !eof(io)
        l = readline(io)
        vals = split(l, ',')
        sval = parse(CL, join(vals[1:end-1], ' '))
        fval = parse(FI, vals[end])
        push!(flt, sval, fval)
    end
    return flt
end

"Write a combined boundary matrix to a text file"
function writeboundarymatrix(io::IO, bm::Vector, zeroindex = true)
    for smplx in bm
        if length(smplx) == 0
            write(io, "0")
        else
            write(io, "$(length(smplx)-1)")
        end
        for i in smplx
            write(io, " $(zeroindex ? i-1 : i)")
        end
        write(io, 0x0A)
    end
end

function writebin(io::IO, flt::Filtration)
    cplx = complex(flt)
    for (d, ci, fv) in order(flt)
        write(io, UInt8(d))
        for k in indices(cplx[ci,d])
            write(io, UInt32(k))
        end
        write(io, fv)
    end
end

function readbin(io::IO, ::Type{Filtration{C,FI}}, ::Type{CL}, ::Type{ST}; maxoutdim=1) where {C <: AbstractComplex, FI, CL<:AbstractCell, ST}
    flt = Filtration(C,FI)
    while !eof(io)
        d = read(io, UInt8)
        vals = ST[read(io, UInt32) for i in 1:(d+1)]
        fval = read(io, FI)
        if d <= maxoutdim
            push!(flt, CL(vals), fval)
        end
    end
    return flt
end


text:
"""
    struct CorrelationFunction

Type representing the two-time first-order correlation function of two operators.
"""
struct CorrelationFunction{OP1,OP2,OP0,DE0,DE,S}
    op1::OP1
    op2::OP2
    op2_0::OP0
    de0::DE0
    de::DE
    steady_state::S
end

"""
    CorrelationFunction(op1,op2,de0;steady_state=false,add_subscript=0,mix_choice=maximum)

The first-order two-time correlation function of two operators.

The first-order two-time correlation function of `op1` and `op2` evolving under
the system `de0`. The keyword `steady_state` determines whether the original
system `de0` was evolved up to steady state. The arguments `add_subscript`
defines the subscript added to the name of `op2` representing the constant time.

Note that the correlation function is stored in the first index of the underlying
system of equations.
"""
function CorrelationFunction(op1,op2,de0::AbstractMeanfieldEquations;
                            steady_state=false, add_subscript=0,
                            filter_func=nothing, mix_choice=maximum,
                            iv=SymbolicUtils.Sym{Real}(:τ),
                            order=nothing,
                            simplify=true, kwargs...)
    h1 = hilbert(op1)
    h2 = _new_hilbert(hilbert(op2), acts_on(op2))
    h = h1⊗h2

    H0 = de0.hamiltonian
    J0 = de0.jumps
    Jd0 = de0.jumps_dagger

    op1_ = _new_operator(op1, h)
    op2_ = _new_operator(op2, h, length(h.spaces); add_subscript=add_subscript)
    op2_0 = _new_operator(op2, h)
    H = _new_operator(H0, h)
    J = [_new_operator(j, h) for j in J0]
    Jd = [_new_operator(j, h) for j in Jd0]
    lhs_new = [_new_operator(l, h) for l in de0.states]

    order_ = if order===nothing
        if de0.order===nothing
            de0.order
            order_lhs = maximum(get_order(l) for l in de0.states)
            order_corr = get_order(op1_*op2_)
            max(order_lhs, order_corr)
        else
            de0.order
        end
    else
        order
    end
    op_ = op1_*op2_
    @assert get_order(op_) <= order_

    de = meanfield(op_,H,J;Jdagger=Jd,rates=de0.rates,iv=iv,order=order_)
    _complete_corr!(de, length(h.spaces), lhs_new, order_, steady_state;
                            filter_func=filter_func,
                            mix_choice=mix_choice,
                            simplify=simplify,
                            kwargs...)

    varmap = make_varmap(lhs_new, de0.iv)
    de0_ = begin
        eqs = Symbolics.Equation[]
        eqs_op = Symbolics.Equation[]
        ops = map(undo_average, lhs_new)
        for i=1:length(de0.equations)
            rhs = _new_operator(de0.equations[i].rhs, h)
            rhs_op = _new_operator(de0.operator_equations[i].rhs, h)
            push!(eqs, Symbolics.Equation(lhs_new[i], rhs))
            push!(eqs_op, Symbolics.Equation(ops[i], rhs_op))
        end
        MeanfieldEquations(eqs,eqs_op,lhs_new,ops,H,J,Jd,de0.rates,de0.iv,varmap,order_)
    end

    return CorrelationFunction(op1_, op2_, op2_0, de0_, de, steady_state)
end

"""
    correlation_u0(c::CorrelationFunction, u_end)

Find the vector containing the correct initial values when numerical solving
the time evolution for the correlation function.

See also: [`CorrelationFunction`](@ref) [`correlation_p0`](@ref)
"""
function correlation_u0(c::CorrelationFunction, u_end)
    a0 = c.op2_0
    a1 = c.op2
    subs = Dict(a1=>a0)
    ops = c.de.operators
    lhs = [average(substitute(op, subs)) for op in ops]
    u0 = complex(eltype(u_end))[]
    lhs0 = c.de0.states
    τ = MTK.get_iv(c.de)
    keys = []
    for j=1:length(lhs)
        l=lhs[j]
        l_adj = _adjoint(l)
        if l ∈ Set(lhs0)
            i = findfirst(isequal(l), lhs0)
            push!(u0, u_end[i])
            push!(keys, make_var(c.de.equations[j].lhs, τ))
        elseif l_adj ∈ Set(lhs0)
            i = findfirst(isequal(l_adj), lhs0)
            push!(u0, conj(u_end[i]))
            push!(keys, make_var(c.de.equations[j].lhs, τ))
        elseif l isa Number
            push!(u0, l)
            push!(keys, make_var(c.de.equations[j].lhs, τ))
        else
            check = false
            for i=1:length(lhs0)
                l_ = substitute(l, Dict(lhs0[i] => u_end[i]))
                check = !isequal(l_, l)
                check && (push!(u0, l_); push!(keys, make_var(c.de.equations[i].lhs, τ)); break)
            end
            check || error("Could not find initial value for $l !")
        end
    end
    return keys .=> u0
end

"""
    correlation_p0(c::CorrelationFunction, u_end, ps=Pair[])

Find all occurring steady-state values and add them to a list of parameters to
pass this to the `ODEProblem`.

See also: [`CorrelationFunction`](@ref) [`correlation_u0`](@ref)
"""
function correlation_p0(c::CorrelationFunction, u_end, ps=Pair{Any,Any}[])
    if c.steady_state
        steady_vals = c.de0.states
        steady_params = map(_make_parameter, steady_vals)
        ps′ = (ps..., (steady_params .=> u_end)...)
    else
        # Check if <a0> is contained
        avg = average(c.op2_0)
        conj_avg = _conj(avg)
        if avg ∈ Set(c.de0.states)
            p = _make_parameter(avg)
            idx = findfirst(isequal(avg), c.de0.states)
            ps′ = (ps..., p=>u_end[idx])
        elseif conj_avg ∈ Set(c.de0.states)
            p = _make_parameter(conj_avg)
            idx = findfirst(isequal(conj_avg), c.de0.states)
            ps′ = (ps..., p=>u_end[idx])
        else
            ps′ = ps
        end
    end
    return ps′
end

"""
    struct Spectrum

Type representing the spectrum, i.e. the Fourier transform of a
[`CorrelationFunction`](@ref) in steady state.

To actually compute the spectrum at a frequency `ω`, construct the type on top
of a correlation function and call it with `Spectrum(c)(ω,usteady,p0)`.
"""
struct Spectrum
    corr
    Afunc
    bfunc
    cfunc
    A
    b
    c
end

"""
    Spectrum(c::CorrelationFunction, ps=[]; kwargs...)

Create an instance of [`Spectrum`](@ref) corresponding to the Fourier transform
of the [`CorrelationFunction`](@ref) `c`.

Examples
========
```
julia> c = CorrelationFunction(a',a,de;steady_state=true)
⟨a′*a_0⟩

julia> S = Spectrum(c)
ℱ(⟨a′*a_0⟩)(ω)
```
"""
function Spectrum(c::CorrelationFunction, ps=[]; w=Parameter(:ω), kwargs...)
    c.steady_state || error("Cannot use Laplace transform when not in steady state! Use `CorrelationFunction(op1,op2,de0;steady_state=true)` or try computing the Fourier transform of the time evolution of the correlation function directly.")
    de = c.de
    de0 = c.de0
    lhs = getfield.(de.equations, :lhs)
    rhs = getfield.(de.equations, :rhs)
    lhs0 = getfield.(de0.equations, :lhs)
    A,b,c_,Afunc,bfunc,cfunc = _build_spec_func(w, lhs, rhs, c.op2_0, c.op2, lhs0, ps; kwargs...)
    return Spectrum(c, Afunc, bfunc, cfunc, A, b, c_)
end

"""
    (s::Spectrum)(ω::Real,usteady,ps=[];wtol=0)

From an instance of [`Spectrum`](@ref) `s`, actually compute the spectral power
density at the frequency `ω`. Numerically solves the equation `x=inv(A)*b` where
`x` is the vector containing the Fourier transformed correlation function, i.e.
the spectrum is given by `real(x[1])`.
`A` and `b` are a matrix and a vector, respectively, describing the linear system
of equations that needs to be solved to obtain the spectrum.
The tolerance `wtol=0` specifies in which range the frequency should be treated
as zero, i.e. whenever `abs(ω) <= wtol` the term proportional to `1/(im*ω)` is
neglected to avoid divergences.
"""
function (s::Spectrum)(ω::Real,usteady,ps=[];wtol=0)
    A = s.Afunc[1](ω,usteady,ps)
    b = s.bfunc[1](usteady,ps)
    if abs(ω) <= wtol
        b_ = b
    else
        c = s.cfunc[1](ω,usteady,ps)
        b_ = b .+ c
    end
    return 2*real(getindex(A \ b_, 1))
end

"""
    (s::Spectrum)(ω_ls,usteady,ps=[];wtol=0)

From an instance of [`Spectrum`](@ref) `s`, actually compute the spectral power
density at all frequencies in `ω_ls`.
"""
function (s::Spectrum)(ω_ls,usteady,ps=[];wtol=0)
    s_ = Vector{real(eltype(usteady))}(undef, length(ω_ls))
    A = s.Afunc[1](ω_ls[1],usteady,ps)
    b0 = s.bfunc[1](usteady,ps)
    b = copy(b0)
    c = s.cfunc[1](ω_ls[1],usteady,ps)

    if abs(ω_ls[1]) <= wtol
        s_[1] = 2*real(getindex(A \ b, 1))
    else
        s_[1] = 2*real(getindex(A \ (b .+ c), 1))
    end

    Afunc! = (A,ω) -> s.Afunc[2](A,ω,usteady,ps)
    cfunc! = (c,ω) -> s.cfunc[2](c,ω,usteady,ps)
    @inbounds for i=2:length(ω_ls)
        Afunc!(A,ω_ls[i])
        if abs(ω_ls[i]) <= wtol
            s_[i] = 2*real(getindex(A \ b0, 1))
        else
            cfunc!(c,ω_ls[i])
            @. b = b0 + c
            s_[i] = 2*real(getindex(A \ b, 1))
        end
    end

    return s_
end

# Convert to ODESystem
function MTK.ODESystem(c::CorrelationFunction; kwargs...)
    τ = MTK.get_iv(c.de)

    ps = []
    for eq∈c.de.equations
        MTK.collect_vars!([],ps,eq.rhs,τ)
    end
    unique!(ps)

    if c.steady_state
        steady_vals = c.de0.states
        steady_hashes = map(hash, steady_vals)
        avg = average(c.op2_0)
        h = hash(avg)
        avg_adj = _adjoint(avg)
        h′ = hash(avg_adj)
        idx = findfirst(isequal(h), steady_hashes)
        if idx === nothing
            idx_ = findfirst(isequal(h′), steady_hashes)
            if idx_ === nothing
                de = c.de
            else
                subs = Dict(average(c.op2) => _adjoint(steady_vals[idx_]))
                de = substitute(c.de, subs)
            end
        else
            subs = Dict(average(c.op2) => steady_vals[idx])
            de = substitute(c.de, subs)
        end
        ps_ = [ps..., steady_vals...]
    else
        avg = average(c.op2_0)
        if avg ∈ Set(c.de0.states)
            avg2 = average(c.op2)
            ps_ = [ps..., avg]
            de = substitute(c.de, Dict(avg2 => avg))
        elseif _conj(avg) ∈ Set(c.de0.states)
            avg2 = average(c.op2)
            ps_ = [ps..., _conj(avg)]
            de = substitute(c.de, Dict(avg2 => avg))
        else
            ps_ = [ps...]
            de = c.de
        end
    end

    ps_avg = filter(x->x isa Average, ps_)
    ps_adj = map(_conj, ps_avg)
    filter!(x->!(x ∈ Set(ps_avg)), ps_adj)
    ps_adj_hash = hash.(ps_adj)

    de_ = deepcopy(de)
    for i=1:length(de.equations)
        lhs = de_.equations[i].lhs
        rhs = substitute_conj(de_.equations[i].rhs, ps_adj, ps_adj_hash)
        de_.equations[i] = Symbolics.Equation(lhs, rhs)
    end

    avg0 = average(c.op2_0)
    if c.steady_state
        steady_params = map(_make_parameter, steady_vals)
        subs_params = Dict(steady_vals .=> steady_params)
        de_ = substitute(de_, subs_params)
    elseif avg0 ∈ Set(c.de0.states)
        avg0_par = _make_parameter(avg0)
        de_ = substitute(de_, Dict(avg0 => avg0_par))
    elseif _conj(avg0) ∈ Set(c.de0.states)
        avg0_par = _make_parameter(_conj(avg0))
        de_ = substitute(de_, Dict(_conj(avg0) => avg0_par))
    end

    eqs = MTK.equations(de_)
    return MTK.ODESystem(eqs, τ; kwargs...)
end

substitute(c::CorrelationFunction, args...; kwargs...) =
    CorrelationFunction(c.op1, c.op2, substitute(c.de0, args...; kwargs...), substitute(c.de, args...; kwargs...))

function _make_parameter(s::Average)
    name = Symbol(string(s))
    sym = Parameter(name)
    return SymbolicUtils.setmetadata(sym, Symbolics.VariableSource, (:_make_parameter, name))
end
_make_parameter(s::SymbolicUtils.Symbolic{<:Parameter}) = s

### Auxiliary functions for CorrelationFunction
function _new_hilbert(h::ProductSpace, aon)
    if length(aon)==1
        return _new_hilbert(h.spaces[aon[1]], 0)
    else
        spaces = [_new_hilbert(h_, 0) for h_ in h.spaces[aon...]]
        return ProductSpace(spaces)
    end
end
_new_hilbert(h::FockSpace, aon) = FockSpace(Symbol(h.name, 0))
_new_hilbert(h::NLevelSpace, aon) = NLevelSpace(Symbol(h.name, 0), h.levels, h.GS)

function _new_operator(op::Destroy, h, aon=op.aon; add_subscript=nothing)
    if isnothing(add_subscript)
        Destroy(h, op.name, aon; op.metadata)
    else
        Destroy(h, Symbol(op.name, :_, add_subscript), aon; op.metadata)
    end
end
function _new_operator(op::Create, h, aon=op.aon; add_subscript=nothing)
    if isnothing(add_subscript)
        Create(h, op.name, aon; op.metadata)
    else
        Create(h, Symbol(op.name, :_, add_subscript), aon; op.metadata)
    end
end
function _new_operator(t::Transition, h, aon=t.aon; add_subscript=nothing)
    if isnothing(add_subscript)
        Transition(h, t.name, t.i, t.j, aon; t.metadata)
    else
        Transition(h, Symbol(t.name, :_, add_subscript), t.i, t.j, aon; t.metadata)
    end
end
_new_operator(x::Number, h, aon=nothing; kwargs...) = x
function _new_operator(t, h, aon=nothing; kwargs...)
    if SymbolicUtils.istree(t)
        args = []
        if isnothing(aon)
            for arg in SymbolicUtils.arguments(t)
                push!(args, _new_operator(arg, h; kwargs...))
            end
        else
            for arg in SymbolicUtils.arguments(t)
                push!(args, _new_operator(arg,h,aon; kwargs...))
            end
        end
        f = SymbolicUtils.operation(t)
        return f(args...)
    else
        return t
    end
end
function _new_operator(avg::SymbolicUtils.Term{<:AvgSym}, h, aon=nothing; kwargs...)
    op = SymbolicUtils.arguments(avg)[1]
    if isnothing(aon)
        _average(_new_operator(op, h; kwargs...))
    else
        _average(_new_operator(op, h, aon; kwargs...))
    end
end

function _complete_corr!(de,aon0,lhs_new,order,steady_state;
                                mix_choice=maximum,
                                simplify=true,
                                filter_func=nothing,
                                kwargs...)
    vs = de.states
    H = de.hamiltonian
    J = de.jumps
    Jd = de.jumps_dagger
    rates = de.rates

    vhash = map(hash, vs)
    vs′ = map(_conj, vs)
    vs′hash = map(hash, vs′)
    filter!(!in(vhash), vs′hash)
    missed = find_missing(de.equations, vhash, vs′hash; get_adjoints=false)

    vhash_new = map(hash, lhs_new)
    vhash_new′ = map(hash, _adjoint.(lhs_new))
    filter!(!in(vhash_new), vhash_new′)

    function _filter_aon(x) # Filter values that act only on Hilbert space representing system at time t0
        aon = acts_on(x)
        if aon0 in aon
            length(aon)==1 && return false
            return true
        end
        if steady_state # Include terms without t0-dependence only if the system is not in steady state
            h = hash(x)
            return !(h∈vhash_new || h∈vhash_new′)
        else
            return true
        end
    end
    filter!(_filter_aon, missed)
    isnothing(filter_func) || filter!(filter_func, missed) # User-defined filter

    filter_reds = de isa ScaledMeanfieldEquations
    filter_reds && filter_redundants!(missed, de.scale_aons, de.names)

    while !isempty(missed)
        ops_ = [SymbolicUtils.arguments(m)[1] for m in missed]
        me = meanfield(ops_,H,J;
                                Jdagger=Jd,
                                rates=rates,
                                simplify=simplify,
                                order=order,
                                iv=de.iv,
                                kwargs...)

        _append!(de, me)

        vhash_ = hash.(me.states)
        vs′hash_ = hash.(_conj.(me.states))
        append!(vhash, vhash_)
        for i=1:length(vhash_)
            vs′hash_[i] ∈ vhash_ || push!(vs′hash, vs′hash_[i])
        end

        missed = find_missing(me.equations, vhash, vs′hash; get_adjoints=false)
        filter!(_filter_aon, missed)
        isnothing(filter_func) || filter!(filter_func, missed) # User-defined filter
        filter_reds && filter_redundants!(missed, de.scale_aons, de.names)
    end

    if !isnothing(filter_func)
        # Find missing values that are filtered by the custom filter function,
        # but still occur on the RHS; set those to 0
        missed = find_missing(de.equations, vhash, vs′hash; get_adjoints=false)
        filter!(!filter_func, missed)
        filter_reds && filter_redundants!(missed, de.scale_aons, de.names)
        missed_adj = map(_adjoint, missed)
        subs = Dict(vcat(missed, missed_adj) .=> 0)
        for i=1:length(de.equations)
            de.equations[i] = substitute(de.equations[i], subs)
            de.states[i] = de.equations[i].lhs
        end
    end

    return de
end

### Auxiliary functions for Spectrum

function _build_spec_func(ω, lhs, rhs, a1, a0, steady_vals, ps=[])
    s = Dict(a0=>a1)
    ops = [SymbolicUtils.arguments(l)[1] for l in lhs]

    b = [average(substitute(op, s)) for op in ops] # Initial values
    c = [SymbolicUtils.simplify(c_ / (1.0im*ω)) for c_ in _find_independent(rhs, a0)]
    aon0 = acts_on(a0)
    @assert length(aon0)==1
    rhs_ = _find_dependent(rhs, aon0[1])
    Ax = [im*ω*lhs[i] - rhs_[i] for i=1:length(lhs)] # Element-wise form of A*x

    # Substitute <a0> by steady-state average <a>
    s_avg = Dict(average(a0) => average(a1))
    Ax = [substitute(A, s_avg) for A∈Ax]
    c = [substitute(c_, s_avg) for c_∈c]

    # Compute symbolic A column-wise by substituting unit vectors into element-wise form of A*x
    A = Matrix{Any}(undef, length(Ax), length(Ax))
    for i=1:length(Ax)
        subs_vals = zeros(length(Ax))
        subs_vals[i] = 1
        subs = Dict(lhs .=> subs_vals)
        A_i = [SymbolicUtils.simplify(substitute(Ax[j],subs)) for j=1:length(Ax)]
        A[:,i] = A_i
    end

    # Substitute conjugates
    vs_adj = map(_conj, steady_vals)
    filter!(x->!(x ∈ Set(steady_vals)), vs_adj)
    vs′hash = map(hash, vs_adj)
    A = [substitute_conj(A_,vs_adj,vs′hash) for A_∈A]
    b = [substitute_conj(b_,vs_adj,vs′hash) for b_∈b]
    c = [substitute_conj(c_,vs_adj,vs′hash) for c_∈c]

    # Keep Symbolics.unflatten_long_ops from stepping into symbolic average
    # functions by substituting. This can be removed once averages store
    # operators as metadata
    A1 = [_substitute_vars(A_) for A_∈A]
    b1 = [_substitute_vars(b_) for b_∈b]
    c1 = [_substitute_vars(c_) for c_∈c]

    # Build functions
    Afunc = Symbolics.build_function(A1, ω, steady_vals, ps; expression=false)
    bfunc = Symbolics.build_function(b1, steady_vals, ps; expression=false)
    cfunc = Symbolics.build_function(c1, ω, steady_vals, ps; expression=false)

    return A, b, c, Afunc, bfunc, cfunc
end

function _substitute_vars(t::SymbolicUtils.Symbolic)
    if SymbolicUtils.istree(t)
        f = SymbolicUtils.operation(t)
        if f === sym_average
            sym = Symbol(string(t))
            return SymbolicUtils.setmetadata(SymbolicUtils.Sym{Complex}(sym),
                Symbolics.VariableSource, (:_substitute_vars, sym))
        else
            args = [_substitute_vars(arg) for arg∈SymbolicUtils.arguments(t)]
            return SymbolicUtils.similarterm(t, f, args)
        end
    else
        return t
    end
end
_substitute_vars(x::Number) = x

_find_independent(rhs::Vector, a0) = [_find_independent(r, a0) for r in rhs]
function _find_independent(r, a0)
    if SymbolicUtils.is_operation(+)(r)
        args_ind = []
        aon0 = acts_on(a0)
        for arg in SymbolicUtils.arguments(r)
            aon = acts_on(arg)
            (aon0 in acts_on(arg) && length(aon)>1) || push!(args_ind, arg)
        end
        isempty(args_ind) && return 0
        return +(args_ind...)
    else
        aon_r = acts_on(r)
        aon0 = acts_on(a0)
        same = (length(aon_r)==length(aon0)) && all(a ∈ aon0 for a∈aon_r) &&
                all(a ∈ aon_r for a ∈ aon0)
        same && return 0
        return r
    end
end

_find_dependent(rhs::Vector, a0) = [_find_dependent(r, a0) for r in rhs]
function _find_dependent(r, a0)
    if SymbolicUtils.is_operation(+)(r)
        args = []
        for arg in SymbolicUtils.arguments(r)
            aon = acts_on(arg)
            (a0 in aon) && length(aon)>1 && push!(args, arg)
        end
        isempty(args) && return 0
        return +(args...)
    else
        aon_r = acts_on(r)
        aon0 = acts_on(a0)
        same = (length(aon_r)==length(aon0)) && all(a ∈ aon0 for a∈aon_r) &&
                all(a ∈ aon_r for a ∈ aon0)
        same || return 0
        return r
    end
end


text:
using Ipopt
using Test
using PowerModels

include("../src/PWF.jl")

include("test_functions.jl")

ipopt = optimizer_with_attributes(Ipopt.Optimizer, "tol" => 0.0001, "print_level" => 0)

@testset "PWF" begin
    include("test_pwf.jl")
    include("test_pwf_to_powermodels.jl")
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#     Functions to manage the Jetson GPIOs.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export init_gpio, gpio_get_mode, gpio_set_mode, gpio_read, gpio_clear, gpio_set,
       gpio_value

################################################################################
#                                Initialization
################################################################################

"""
    init_gpio()

Initialize the GPIO.

"""
function init_gpio()
    # Open and map /dev/gpiomem.
    try
        gpiomem_io  = open("/dev/gpiomem", "w+")
        gpiomem_map = Mmap.mmap(gpiomem_io, Vector{UInt32}, (GPIOMEM_SIZE,), 0,
                                grow = false)

        objects.gpiomem_io  = gpiomem_io
        objects.gpiomem_map = gpiomem_map
        objects.gpio_init   = true
    catch
        error("Could not map `/dev/gpiomem`. Make sure you have the required permissions.")
    end

    return nothing
end

################################################################################
#                                  GPIO State
################################################################################

"""
    gpio_get_mode(gpio::Int)

Return the state of the GPIO `gpio`.

"""
function gpio_get_mode(gpio::Int)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (0 ≤ gpio ≤ 27) "GPIO out of range."

    map = objects.gpiomem_map

    @inbounds begin
        ind   = div(gpio, 10) + 1
        g     = 3*(gpio % 10)
        mask  = 7 << g
        state = (map[ind] & mask) >> g

        if haskey(GPIO_MODE_GET, state)
            return GPIO_MODE_GET[state]
        else
            return :undef
        end
    end
end

"""
    gpio_set_mode(gpio, mode::Symbol)

Set the mode of the GPIO `gpio` to `mode`. `gpio` can be an integer of an
`AbstractVector` with the list of GPIOs that will be set.

`mode` can be:

* `:out`: GPIO direction will be set to out.
* `:in`: GPIO direction will be set to in.
* `:alt0`: GPIO will be set to alternate function 0.
* `:alt1`: GPIO will be set to alternate function 1.
* `:alt2`: GPIO will be set to alternate function 2.
* `:alt3`: GPIO will be set to alternate function 3.
* `:alt4`: GPIO will be set to alternate function 4.
* `:alt5`: GPIO will be set to alternate function 5.

"""
function gpio_set_mode(gpio::Int, mode::Symbol)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (0 ≤ gpio ≤ 27) "GPIO out of range."
    @assert haskey(GPIO_MODE_SET, mode) "Unknown GPIO mode."

    map = objects.gpiomem_map

    @inbounds begin
        ind  = div(gpio, 10) + 1
        g    = 3*(gpio % 10)

        # We must always clear the GPIO mode before setting the new mode. This
        # is done by temporarily setting the GPIO to `:in`.
        map[ind] &= ~(7 << g)

        if (mode != :in)
            code = GPIO_MODE_SET[mode]
            map[ind] |= (UInt32(code) << g)
        end
    end

    return nothing
end

function gpio_set_mode(gpio::AbstractVector{Int}, mode::Symbol)
    for g ∈ gpio
        gpio_set_mode(g, mode)
    end

    return nothing
end

################################################################################
#                                     Read
################################################################################

"""
    gpio_read(gpio::Int)

Read the GPIO `gpio`. The returned value is boolean.

"""
function gpio_read(gpio::Int)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (0 ≤ gpio ≤ 27) "GPIO out of range."

    @inbounds begin
        if (objects.gpiomem_map[13+1] & (1 << gpio)) > 0
            return true
        else
            return false
        end
    end
end

################################################################################
#                                Set and clear
################################################################################

"""
    gpio_clear(gpio::Int)

Clear GPIO `gpio`.

"""
@inline function gpio_clear(gpio::Int)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (0 ≤ gpio ≤ 27) "GPIO out of range."

    @inbounds begin
        objects.gpiomem_map[10+1] = 1 << gpio
    end

    return nothing
end

@inline function gpio_clear(gpio::AbstractVector{Int})
    for g in gpio
        gpio_clear(g)
    end

    return nothing
end

"""
    gpio_set(gpio::Int)

Set GPIO `gpio`.

"""
@inline function gpio_set(gpio::Int)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (0 ≤ gpio ≤ 27) "GPIO out of range."

    @inbounds begin
        objects.gpiomem_map[7+1] = 1 << gpio
    end

    return nothing
end

@inline function gpio_set(gpio::AbstractVector{Int})
    for g in gpio
        gpio_set(g)
    end

    return nothing
end

"""
    gpio_value(v::Integer)

Apply the value `v` to all GPIOs. The bits that are `0` will be cleared, and the
bits that are `1` will be set.

    gpio_value(v::BitVector)

Convert the bit array `v` to `Int` and call `gpio_value(::Int)`.

"""
function gpio_value(v::Integer)
    @assert objects.gpio_init "GPIO not initialized. Run init_gpio()."
    @assert (v ≥ 0) "`v` cannot be negative."

    @inbounds begin
        # Set the GPIO.
        objects.gpiomem_map[7+1]  = v & 0xFFFFFFF

        # Clear the GPIOs.
        objects.gpiomem_map[10+1] = (~UInt32(v)) & 0xFFFFFFF
    end

    return nothing
end

gpio_value(v::BitVector) = gpio_value(v.chunks[1])


text:
module printing

export set_printing

import Base: show, summary

using ..bases, ..states
using ..operators, ..operators_dense, ..operators_sparse
using ..operators_lazytensor, ..operators_lazysum, ..operators_lazyproduct
using ..spin, ..fock, ..nlevel, ..particle, ..subspace, ..manybody, ..sparsematrix
using SparseArrays

"""
    QuantumOptics.set_printing(; standard_order, rounding_tol)

Set options for REPL output.

# Arguments
* `standard_order=false`: For performance reasons, the order of the tensor
    product is inverted, i.e. `tensor(a, b)=kron(b, a)`. When changing this
    to `true`, the output shown in the REPL will exhibit the correct order.
* `rounding_tol=1e-17`: Tolerance for floating point errors shown in the output.
"""
function set_printing(; standard_order::Bool=_std_order, rounding_tol::Real=_round_tol)
    global _std_order = standard_order
    global _round_tol = rounding_tol
    global machineprecorder = Int32(round(-log10(_round_tol), digits=0))
    nothing
end
set_printing(standard_order=false, rounding_tol=1e-17)

function show(stream::IO, x::GenericBasis)
    if length(x.shape) == 1
        write(stream, "Basis(dim=$(x.shape[1]))")
    else
        s = replace(string(x.shape), " " => "")
        write(stream, "Basis(shape=$s)")
    end
end

function show(stream::IO, x::CompositeBasis)
    write(stream, "[")
    for i in 1:length(x.bases)
        show(stream, x.bases[i])
        if i != length(x.bases)
            write(stream, " ⊗ ")
        end
    end
    write(stream, "]")
end

function show(stream::IO, x::SpinBasis)
    d = denominator(x.spinnumber)
    n = numerator(x.spinnumber)
    if d == 1
        write(stream, "Spin($n)")
    else
        write(stream, "Spin($n/$d)")
    end
end

function show(stream::IO, x::FockBasis)
    write(stream, "Fock(cutoff=$(x.N))")
end

function show(stream::IO, x::NLevelBasis)
    write(stream, "NLevel(N=$(x.N))")
end

function show(stream::IO, x::PositionBasis)
    write(stream, "Position(xmin=$(x.xmin), xmax=$(x.xmax), N=$(x.N))")
end

function show(stream::IO, x::MomentumBasis)
    write(stream, "Momentum(pmin=$(x.pmin), pmax=$(x.pmax), N=$(x.N))")
end

function show(stream::IO, x::SubspaceBasis)
    write(stream, "Subspace(superbasis=$(x.superbasis), states:$(length(x.basisstates)))")
end

function show(stream::IO, x::ManyBodyBasis)
    write(stream, "ManyBody(onebodybasis=$(x.onebodybasis), states:$(length(x.occupations)))")
end

summary(io::IO, x::Ket) = print(io, "Ket(dim=$(length(x.basis)))\n  basis: $(x.basis)\n")
function show(stream::IO, x::Ket)
    summary(stream, x)
    if !_std_order
        Base.print_array(stream, round.(x.data; digits=machineprecorder))
    else
        showarray_stdord(stream, round.(x.data; digits=machineprecorder), x.basis.shape, false, header=false)
    end
end

summary(io::IO, x::Bra) = print(io, "Bra(dim=$(length(x.basis)))\n  basis: $(x.basis)\n")
function show(stream::IO, x::Bra)
    summary(stream, x)
    if !_std_order
        Base.print_array(stream, round.(x.data; digits=machineprecorder))
    else
        showarray_stdord(stream, round.(x.data; digits=machineprecorder), x.basis.shape, false, header=false)
    end
end

function summary(stream::IO, x::AbstractOperator)
    print(stream, "$(typeof(x).name.name)(dim=$(length(x.basis_l))x$(length(x.basis_r)))\n")
    if bases.samebases(x)
        print(stream, "  basis: ")
        show(stream, basis(x))
    else
        print(stream, "  basis left:  ")
        show(stream, x.basis_l)
        print(stream, "\n  basis right: ")
        show(stream, x.basis_r)
    end
end

show(stream::IO, x::AbstractOperator) = summary(stream, x)

function show(stream::IO, x::DenseOperator)
    summary(stream, x)
    print(stream, "\n")
    if !_std_order
        if !haskey(stream, :compact)
            stream = IOContext(stream, :compact => true)
        end
        Base.print_array(stream, round.(x.data; digits=machineprecorder))
    else
        showarray_stdord(stream, round.(x.data; digits=machineprecorder), x.basis_l.shape, x.basis_r.shape, false, header=false)
    end
end

function show(stream::IO, x::SparseOperator)
    summary(stream, x)
    if nnz(x.data) == 0
        print(stream, "\n    []")
    else
        if !_std_order
            if !haskey(stream, :compact)
                stream = IOContext(stream, :compact => true)
            end
            show(stream, round.(x.data; digits=machineprecorder))
        else
            showsparsearray_stdord(stream, round.(x.data; digits=machineprecorder), x.basis_l.shape, x.basis_r.shape)
        end
    end
end

function show(stream::IO, x::LazyTensor)
    summary(stream, x)
    print(stream, "\n  operators: $(length(x.operators))")
    s = replace(string(x.indices), " " => "")
    print(stream, "\n  indices: $s")
end

function show(stream::IO, x::Union{LazySum, LazyProduct})
    summary(stream, x)
    print(stream, "\n  operators: $(length(x.operators))")
end

"""
    ind2Nary(m::Int, dims::Vector{Int})

The inverse of `Nary2ind`.

# Example
```
julia> dims = [2,2,3];

julia> for i in 1:prod(dims)
           println(i,": ", ind2Nary(i,dims))
       end
1: [0, 0, 0]
2: [0, 0, 1]
3: [0, 0, 2]
4: [0, 1, 0]
5: [0, 1, 1]
6: [0, 1, 2]
7: [1, 0, 0]
8: [1, 0, 1]
9: [1, 0, 2]
10: [1, 1, 0]
11: [1, 1, 1]
12: [1, 1, 2]
```
"""
function ind2Nary(m::Int, dims::Vector{Int})
    m = m - 1
    nq = length(dims)
    ar = zeros(Int, nq)
    product = prod(dims[2:end])
    for ith in 1:nq-1
        d = div(m, product)
        m = m - d * product
        product = div(product, dims[ith+1])
        ar[ith] = d
    end
    ar[end] = m
    return ar
end

"""
    Nary2ind(x, dims) -> index

Convert composite N-arys to index.

# Example
```
julia> dims = [2,2,3];

julia> Nary2ind([1,1,0], dims)
10

julia> for i in 1:prod(dims)
           println(i,": ", Nary2ind(ind2Nary(i,dims),dims), ": ", ind2Nary(i,dims))
       end
1: 1: [0, 0, 0]
2: 2: [0, 0, 1]
3: 3: [0, 0, 2]
4: 4: [0, 1, 0]
5: 5: [0, 1, 1]
6: 6: [0, 1, 2]
7: 7: [1, 0, 0]
8: 8: [1, 0, 1]
9: 9: [1, 0, 2]
10: 10: [1, 1, 0]
11: 11: [1, 1, 1]
12: 12: [1, 1, 2]
```
"""
function Nary2ind(x::Vector{Int}, dims::Vector{Int})
    tmp = 0
    @assert length(x) == length(dims)
    nterms = length(x)
    tp = prod(dims[2:end])
    for i in 1:nterms-1
        tmp += x[i] * tp
        tp = div(tp, dims[i+1])
    end
    tmp += x[end] + 1
end

"""
    mirror_world_index(idx, dims)

Convert index of standard order to that of inversed order.
This function is named after the book, '鏡の中の物理学' (Physics in the mirror),
written by Tomonaga Shin'ichirō (Japanese physicist).
"""
function mirror_world_index(idx::Int, dims::Vector{Int})
    return Nary2ind( reverse(ind2Nary(idx, dims)), reverse(dims)  )
end

# Following program is modified:
# julia/base/show.jl
# https://github.com/JuliaLang/julia/blob/5cd144ffa328fdd4cd9e983b616c7205f9bb4f51/base/show.jl
# and
# julia/base/sparse/sparsematrix.jl
# https://github.com/JuliaLang/julia/blob/78831902cc412e298c5fcc94dae1e4382c8e7cd0/base/sparse/sparsematrix.jl

# Copyright (c) 2009-2018: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
# and other contributors:
#
# https://github.com/JuliaLang/julia/contributors
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""
`alignment_std(X, ldims, rdims, rows, cols, cols_if_complete, cols_otherwise, sep)` returns the
alignment for specified parts of array `X`, returning the (left,right) info.
It will look in X's `rows`, `cols` (both lists of indices)
and figure out what's needed to be fully aligned, for example looking all
the way down a column and finding out the maximum size of each element.
Parameter `sep::Integer` is number of spaces to put between elements.
`cols_if_complete` and `cols_otherwise` indicate screen width to use.
Alignment is reported as a vector of (left,right) tuples, one for each
column going across the screen.
"""
function alignment_std(io::IO, X::AbstractVecOrMat, ldims::Vector, rdims::Vector,
        rows::AbstractVector, cols::AbstractVector,
        cols_if_complete::Integer, cols_otherwise::Integer, sep::Integer)
    a = Tuple{Int, Int}[]
    for j in cols # need to go down each column one at a time
        l = r = 0
        for i in rows # plumb down and see what largest element sizes are
            aij = Base.alignment(io, X[mirror_world_index(i, ldims), mirror_world_index(j, rdims)])
            l = max(l, aij[1]) # left characters
            r = max(r, aij[2]) # right characters
        end
        push!(a, (l, r)) # one tuple per column of X, pruned to screen width
        if length(a) > 1 && sum(map(sum,a)) + sep*length(a) >= cols_if_complete
            pop!(a) # remove this latest tuple if we're already beyond screen width
            break
        end
    end
    if 1 < length(a) < length(axes(X,2))
        while sum(map(sum,a)) + sep*length(a) >= cols_otherwise
            pop!(a)
        end
    end
    return a
end

"""
`print_matrix_row_std(io, X, A, ldims, rdims, i, cols, sep)` produces the aligned output for
a single matrix row X[i, cols] where the desired list of columns is given.
The corresponding alignment A is used, and the separation between elements
is specified as string sep.
`print_matrix_row_std` will also respect compact output for elements.
"""
function print_matrix_row_std(io::IO,
        X::AbstractVecOrMat, A::Vector, ldims::Vector, rdims::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString)
    isempty(A) || first(axes(cols,1)) == 1 || throw(DimensionMismatch("indices of cols ($(axes(cols,1))) must start at 1"))
    for k = 1:length(A)
        j = cols[k]
        x = X[mirror_world_index(i, ldims), mirror_world_index(j, rdims)]
        a = Base.alignment(io, x)
        sx = sprint(show, x, context=io)
        l = repeat(" ", A[k][1]-a[1]) # pad on left and right as needed
        r = repeat(" ", A[k][2]-a[2])
        prettysx = Base.replace_in_print_matrix(X,mirror_world_index(i, ldims),mirror_world_index(j, rdims),sx)
        print(io, l, prettysx, r)
        if k < length(A); print(io, sep); end
    end
end

function print_matrix_std(io::IO, X::AbstractVecOrMat, ldims::Vector, rdims::Vector,
                      pre::AbstractString = " ",  # pre-matrix string
                      sep::AbstractString = "  ", # separator between elements
                      post::AbstractString = "",  # post-matrix string
                      hdots::AbstractString = "  \u2026  ",
                      vdots::AbstractString = "\u22ee",
                      ddots::AbstractString = "  \u22f1  ",
                      hmod::Integer = 5, vmod::Integer = 5)
    if !get(io, :limit, false)
        screenheight = screenwidth = typemax(Int)
    else
        sz = displaysize(io)
        screenheight, screenwidth = sz[1] - 4, sz[2]
    end
    screenwidth -= length(pre) + length(post)
    presp = repeat(" ", length(pre))  # indent each row to match pre string
    postsp = ""
    @assert textwidth(hdots) == textwidth(ddots)
    sepsize = length(sep)
    rowsA, colsA = axes(X,1), axes(X,2)
    m, n = length(rowsA), length(colsA)
    # To figure out alignments, only need to look at as many rows as could
    # fit down screen. If screen has at least as many rows as A, look at A.
    # If not, then we only need to look at the first and last chunks of A,
    # each half a screen height in size.
    halfheight = div(screenheight,2)
    if m > screenheight
        rowsA = [rowsA[1:halfheight]; rowsA[m-div(screenheight-1,2)+1:m]]
    end
    # Similarly for columns, only necessary to get alignments for as many
    # columns as could conceivably fit across the screen
    maxpossiblecols = div(screenwidth, 1+sepsize)
    if n > maxpossiblecols
        colsA = [colsA[1:maxpossiblecols]; colsA[(n-maxpossiblecols+1):n]]
    end
    A = alignment_std(io, X, ldims, rdims, rowsA, colsA, screenwidth, screenwidth, sepsize)
    # Nine-slicing is accomplished using print_matrix_row_std repeatedly
    if m <= screenheight # rows fit vertically on screen
        if n <= length(A) # rows and cols fit so just print whole matrix in one piece
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row_std(io, X,A,ldims,rdims,i,colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        else # rows fit down screen but cols don't, so need horizontal ellipsis
            c = div(screenwidth-length(hdots)+1,2)+1  # what goes to right of ellipsis
            Ralign = reverse(alignment_std(io, X, ldims, rdims, rowsA, reverse(colsA), c, c, sepsize)) # alignments for right
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)
            Lalign = alignment_std(io, X, ldims, rdims, rowsA, colsA, c, c, sepsize) # alignments for left of ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row_std(io, X,Lalign,ldims,rdims,i,colsA[1:length(Lalign)],sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)))
                print_matrix_row_std(io, X,Ralign,ldims,rdims,i,n-length(Ralign).+colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        end
    else # rows don't fit so will need vertical ellipsis
        if n <= length(A) # rows don't fit, cols do, so only vertical ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row_std(io, X,A,ldims, rdims,i,colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    Base.print_matrix_vdots(io, vdots,A,sep,vmod,1)
                    println(io, i == last(rowsA) ? post : postsp)
                end
            end
        else # neither rows nor cols fit, so use all 3 kinds of dots
            c = div(screenwidth-length(hdots)+1,2)+1
            Ralign = reverse(alignment_std(io, X, ldims, rdims, rowsA, reverse(colsA), c, c, sepsize))
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)
            Lalign = alignment_std(io, X, ldims, rdims, rowsA, colsA, c, c, sepsize)
            r = mod((length(Ralign)-n+1),vmod) # where to put dots on right half
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row_std(io, X,Lalign,ldims, rdims,i,colsA[1:length(Lalign)],sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)))
                print_matrix_row_std(io, X,Ralign,ldims, rdims,i,n-length(Ralign).+colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    Base.print_matrix_vdots(io, vdots,Lalign,sep,vmod,1)
                    print(io, ddots)
                    Base.print_matrix_vdots(io, vdots,Ralign,sep,vmod,r)
                    println(io, i == last(rowsA) ? post : postsp)
                end
            end
        end
    end
end

function showarray_stdord(io::IO, X::AbstractVecOrMat, ldims::Vector, rdims::Vector, repr::Bool = true; header = true)
    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if !isempty(X)
        punct = (" ", "  ", "")
        print_matrix_std(io, X, ldims, rdims, punct...)
    end
end
showarray_stdord(io::IO, X::Vector, dims::Vector, repr::Bool = true; header = true) = showarray_stdord(io, X, dims, [1], repr; header = header)

function showsparsearray_stdord(io::IO, S::SparseMatrixCSC, ldims::Vector, rdims::Vector)
    if get(io, :limit, false)
        rows = displaysize(io)[1]
        half_screen_rows = div(rows - 8, 2)
    else
        half_screen_rows = typemax(Int)
    end
    pad = ndigits(max(S.m,S.n))
    sep = "\n  "
    if !haskey(io, :compact) && length(axes(S, 2)) > 1
        io = IOContext(io, :compact => true)
    end

    colval = zeros(Int, nnz(S))
    for col in 1:S.n, k in S.colptr[col] : (S.colptr[col+1]-1)
        colval[k] = col
    end
    rowval = S.rowval
    rowval_std = map(x -> mirror_world_index(x, reverse(ldims)), rowval)
    colval_std = map(x -> mirror_world_index(x, reverse(rdims)), colval)
    idx_nzval = map((x,y,z) -> (x,y,z), rowval_std, colval_std, S.nzval)
    sort!(idx_nzval, by = x -> (x[2], x[1]) )
    rowval_std = map(x -> x[1], idx_nzval)
    colval_std = map(x -> x[2], idx_nzval)
    nzval_std = map(x -> x[3], idx_nzval)

    for (k, val) in enumerate(nzval_std)
        if k < half_screen_rows || k > nnz(S)-half_screen_rows
            print(io, sep, '[', rpad(rowval_std[k], pad), ", ", lpad(colval_std[k], pad), "]  =  ")
            show(io, nzval_std[k])
        elseif k == half_screen_rows
            print(io, sep, '\u22ee')
        end
    end
end

end # module


text:
function add_additional_state_provider(obj::FieldPropagator, arg0::FieldAdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (FieldAdditionalStateProvider,), arg0)
end

function add_event_detector(obj::FieldPropagator, arg0::FieldEventDetector)
    return jcall(obj, "addEventDetector", void, (FieldEventDetector,), arg0)
end

function clear_events_detectors(obj::FieldPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function get_additional_state_providers(obj::FieldPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_attitude_provider(obj::FieldPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_events_detectors(obj::FieldPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::FieldPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::FieldPropagator)
    return jcall(obj, "getGeneratedEphemeris", FieldBoundedPropagator, ())
end

function get_initial_state(obj::FieldPropagator)
    return jcall(obj, "getInitialState", FieldSpacecraftState, ())
end

function get_managed_additional_states(obj::FieldPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_max_date(obj::FieldBoundedPropagator)
    return jcall(obj, "getMaxDate", FieldAbsoluteDate, ())
end

function get_min_date(obj::FieldBoundedPropagator)
    return jcall(obj, "getMinDate", FieldAbsoluteDate, ())
end

function get_mode(obj::FieldPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_pv_coordinates(obj::FieldPVCoordinatesProvider, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function is_additional_state_managed(obj::FieldPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::FieldPropagator, arg0::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate,), arg0)
end

function reset_initial_state(obj::FieldPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end

function set_attitude_provider(obj::FieldPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::FieldPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_master_mode(obj::FieldPropagator, arg0::FieldOrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (FieldOrekitStepHandler,), arg0)
end

function set_master_mode(obj::FieldPropagator, arg0::RealFieldElement, arg1::FieldOrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (RealFieldElement, FieldOrekitFixedStepHandler), arg0, arg1)
end

function set_slave_mode(obj::FieldPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end


text:
# # Intro to AbstractGPs: one-dimensional regression
#
# ## Setup
#
# Loading the necessary packages.

using AbstractGPs
using Distributions
using StatsFuns

using Plots
default(; legend=:outertopright, size=(700, 400))

using Random
Random.seed!(42)  # setting the seed for reproducibility of this notebook
#md nothing #hide

# Load toy regression
# [dataset](https://github.com/GPflow/GPflow/blob/7705cee6723f78066981f27954130daaede55dfc/doc/sphinx/notebooks/basics/data/regression_1D.csv)
# taken from GPflow examples.

x = [
    0.8658165855998895,
    0.6661700880180962,
    0.8049218148148531,
    0.7714303440386239,
    0.14790478354654835,
    0.8666105548197428,
    0.007044577166530286,
    0.026331737288148638,
    0.17188596617099916,
    0.8897812990554013,
    0.24323574561119998,
    0.028590102134105955,
]
y = [
    1.5255314337144372,
    3.6434202968230003,
    3.010885733911661,
    3.774442382979625,
    3.3687639483798324,
    1.5506452040608503,
    3.790447985799683,
    3.8689707574953,
    3.4933565751758713,
    1.4284538820635841,
    3.8715350915692364,
    3.7045949061144983,
]
scatter(x, y; xlabel="x", ylabel="y", legend=false)

# We split the observations into train and test data.

x_train = x[1:8]
y_train = y[1:8]
x_test = x[9:end]
y_test = y[9:end]
#md nothing #hide

# We instantiate a Gaussian process with a Matern kernel. The kernel has
# fixed variance and length scale parameters of default value 1.

f = GP(Matern52Kernel())
#md nothing #hide

# We create a finite dimensional projection at the inputs of the training dataset
# observed under Gaussian noise with variance $\sigma^2 = 0.1$, and compute the
# log-likelihood of the outputs of the training dataset.

fx = f(x_train, 0.1)
logpdf(fx, y_train)

# We compute the posterior Gaussian process given the training data, and calculate the
# log-likelihood of the test dataset.

p_fx = posterior(fx, y_train)
logpdf(p_fx(x_test), y_test)

# We plot the posterior Gaussian process (its mean and a ribbon of 2 standard deviations
# around it) on a grid along with the observations.

scatter(
    x_train,
    y_train;
    xlim=(0, 1),
    xlabel="x",
    ylabel="y",
    title="posterior (default parameters)",
    label="Train Data",
)
scatter!(x_test, y_test; label="Test Data")
plot!(0:0.001:1, p_fx; label=false, ribbon_scale=2)

# ## Markov Chain Monte Carlo
#
# Previously we computed the log likelihood of the untuned kernel parameters of the GP.
# We now also perform approximate inference over said kernel parameters using different
# Markov chain Monte Carlo (MCMC) methods. I.e., we approximate the posterior distribution
# of the kernel parameters with samples from a Markov chain.
#
# We define a function which returns the log-likelihood of the data for different variance
# and inverse lengthscale parameters of the Matern kernel. We ensure that these parameters are
# positive with the softplus function
# ```math
# f(x) = \log (1 + \exp x).
# ```

function gp_loglikelihood(x, y)
    function loglikelihood(params)
        kernel =
            softplus(params[1]) * (Matern52Kernel() ∘ ScaleTransform(softplus(params[2])))
        f = GP(kernel)
        fx = f(x, 0.1)
        return logpdf(fx, y)
    end
    return loglikelihood
end

const loglik_train = gp_loglikelihood(x_train, y_train)
#md nothing #hide

# We define a Gaussian prior for the joint distribution of the two transformed kernel
# parameters. We assume that both parameters are independent with mean 0 and variance 1.

logprior(params) = logpdf(MvNormal(2, 1), params)
#md nothing #hide

# ### Hamiltonian Monte Carlo
#
# We start with a Hamiltonian Monte Carlo (HMC) sampler. More precisely, we use the
# [No-U-Turn sampler (NUTS)](http://www.jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf),
# which is provided by the Julia packages
# [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl/) and
# [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl/).
#
# #### AdvancedHMC
#
# We start with performing inference with AdvancedHMC.

using AdvancedHMC
using ForwardDiff

# Set the number of samples to draw and warmup iterations.

n_samples = 2_000
n_adapts = 1_000
#md nothing #hide

# Define a Hamiltonian system of the log joint probability.

logjoint_train(params) = loglik_train(params) + logprior(params)
metric = DiagEuclideanMetric(2)
hamiltonian = Hamiltonian(metric, logjoint_train, ForwardDiff)
#md nothing #hide

# Define a leapfrog solver, with initial step size chosen heuristically.

initial_params = rand(2)
initial_ϵ = find_good_stepsize(hamiltonian, initial_params)
integrator = Leapfrog(initial_ϵ)
#md nothing #hide

# Define an HMC sampler, with the following components:
# - multinomial sampling scheme,
# - generalised No-U-Turn criteria, and
# - windowed adaption for step-size and diagonal mass matrix

proposal = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))
#md nothing #hide

# We draw samples from the posterior distribution of kernel parameters. These samples
# are in the unconstrained space $\mathbb{R}^2$.

samples, _ = sample(
    hamiltonian, proposal, initial_params, n_samples, adaptor, n_adapts; progress=false
)
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample",
    ylabel="counts",
    layout=2,
    title=["variance" "inverse length scale"],
    legend=false,
)
vline!(mean_samples'; linewidth=2)

# We approximate the log-likelihood of the test data using the posterior Gaussian processes
# for kernels with the sampled kernel parameters. We can observe that there is a significant
# improvement over the log-likelihood of the test data with respect to the posterior
# Gaussian process with default kernel parameters of value 1.

function gp_posterior(x, y, p)
    kernel = softplus(p[1]) * (Matern52Kernel() ∘ ScaleTransform(softplus(p[2])))
    f = GP(kernel)
    return posterior(f(x, 0.1), y)
end

mean(logpdf(gp_posterior(x_train, y_train, p)(x_test), y_test) for p in samples)

# We sample 5 functions from each posterior GP given by the final 100 samples of kernel
# parameters.

plt = plot(; xlim=(0, 1), xlabel="x", ylabel="y", title="posterior (AdvancedHMC)")
for (i, p) in enumerate(samples[(end - 100):end])
    sampleplot!(
        plt,
        0:0.02:1,
        gp_posterior(x_train, y_train, p);
        samples=5,
        seriescolor="red",
        label=(i == 1 ? "samples" : nothing),
    )
end
scatter!(plt, x_train, y_train; label="Train Data", markercolor=1)
scatter!(plt, x_test, y_test; label="Test Data", markercolor=2)
plt

# #### DynamicHMC
#
# We repeat the inference with DynamicHMC. DynamicHMC requires us to
# implement the LogDensityProblems interface for `loglik_train`.

using DynamicHMC
using LogDensityProblems

## Log joint density
function LogDensityProblems.logdensity(ℓ::typeof(loglik_train), params)
    return ℓ(params) + logprior(params)
end

## The parameter space is two-dimensional
LogDensityProblems.dimension(::typeof(loglik_train)) = 2

## `loglik_train` does not allow to evaluate derivatives of
## the log-likelihood function
function LogDensityProblems.capabilities(::Type{<:typeof(loglik_train)})
    return LogDensityProblems.LogDensityOrder{0}()
end

# Now we can draw samples from the posterior distribution of kernel parameters with
# DynamicHMC. Again we use [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
# to compute the derivatives of the log joint density with automatic differentiation.

samples =
    mcmc_with_warmup(
        Random.GLOBAL_RNG,
        ADgradient(:ForwardDiff, loglik_train),
        n_samples;
        reporter=NoProgressReport(),
    ).chain
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample",
    ylabel="counts",
    layout=2,
    title=["variance" "inverse length scale"],
    legend=false,
)
vline!(mean_samples'; linewidth=2)

# Again we can observe that there is a significant improvement over the log-likelihood
# of the test data with respect to the posterior Gaussian process with default kernel
# parameters.

mean(logpdf(gp_posterior(x_train, y_train, p)(x_test), y_test) for p in samples)

# We sample a function from the posterior GP for the final 100 samples of kernel
# parameters.

plt = plot(; xlim=(0, 1), xlabel="x", ylabel="y", title="posterior (DynamicHMC)")
scatter!(plt, x_train, y_train; label="Train Data")
scatter!(plt, x_test, y_test; label="Test Data")
for p in samples[(end - 100):end]
    sampleplot!(plt, 0:0.02:1, gp_posterior(x_train, y_train, p); seriescolor="red")
end
plt

# ### Elliptical slice sampling
#
# Instead of HMC, we use
# [elliptical slice sampling](http://proceedings.mlr.press/v9/murray10a/murray10a.pdf)
# which is provided by the Julia package
# [EllipticalSliceSampling.jl](https://github.com/TuringLang/EllipticalSliceSampling.jl/).

using EllipticalSliceSampling

# We draw 2000 samples from the posterior distribution of kernel parameters.

samples = sample(ESSModel(
    MvNormal(2, 1), # Gaussian prior
    loglik_train,
), ESS(), n_samples; progress=false)
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample",
    ylabel="counts",
    layout=2,
    title=["variance" "inverse length scale"],
)
vline!(mean_samples'; layout=2, labels="mean")

# Again we can observe that there is a significant improvement over the log-likelihood
# of the test data with respect to the posterior Gaussian process with default kernel
# parameters.

mean(logpdf(gp_posterior(x_train, y_train, p)(x_test), y_test) for p in samples)

# We sample a function from the posterior GP for the final 100 samples of kernel
# parameters.

plt = plot(;
    xlim=(0, 1), xlabel="x", ylabel="y", title="posterior (EllipticalSliceSampling)"
)
scatter!(plt, x_train, y_train; label="Train Data")
scatter!(plt, x_test, y_test; label="Test Data")
for p in samples[(end - 100):end]
    sampleplot!(plt, 0:0.02:1, gp_posterior(x_train, y_train, p); seriescolor="red")
end
plt

# ## Variational Inference
#
# Sanity check for the Evidence Lower BOund (ELBO) implemented according to
# M. K. Titsias's _Variational learning of inducing variables in sparse Gaussian processes_.

elbo(VFE(f(rand(5))), fx, y_train)

# We use the LBFGS algorithm to maximize the given ELBO. It is provided by the Julia
# package [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl).

using Optim

# We define a function which returns the negative ELBO for different variance and inverse
# lengthscale parameters of the Matern kernel and different pseudo-points. We ensure that
# the kernel parameters are positive with the softplus function
# ```math
# f(x) = \log (1 + \exp x),
# ```
# and that the pseudo-points are in the unit interval $[0,1]$ with the logistic function
# ```math
# f(x) = \frac{1}{1 + \exp{(-x)}}.
# ```

jitter = 1e-6  # "observing" the latent process with some (small) amount of jitter improves numerical stability

function objective_function(x, y)
    function negative_elbo(params)
        kernel =
            softplus(params[1]) * (Matern52Kernel() ∘ ScaleTransform(softplus(params[2])))
        f = GP(kernel)
        fx = f(x, 0.1)
        z = logistic.(params[3:end])
        approx = VFE(f(z, jitter))
        return -elbo(approx, fx, y)
    end
    return negative_elbo
end
#md nothing #hide

# We randomly initialize the kernel parameters and 5 pseudo points, and minimize the
# negative ELBO with the LBFGS algorithm and obtain the following optimal parameters:

x0 = rand(7)
opt = optimize(objective_function(x_train, y_train), x0, LBFGS())

#-

opt.minimizer

# The optimized value of the variance is

softplus(opt.minimizer[1])

# and of the inverse lengthscale is

softplus(opt.minimizer[2])

# We compute the log-likelihood of the test data for the resulting approximate
# posterior. We can observe that there is a significant improvement over the
# log-likelihood with the default kernel parameters of value 1.

opt_kernel =
    softplus(opt.minimizer[1]) *
    (Matern52Kernel() ∘ ScaleTransform(softplus(opt.minimizer[2])))
opt_f = GP(opt_kernel)
opt_fx = opt_f(x_train, 0.1)
ap = posterior(VFE(opt_f(logistic.(opt.minimizer[3:end]), jitter)), opt_fx, y_train)
logpdf(ap(x_test), y_test)

# We visualize the approximate posterior with optimized parameters.

scatter(
    x_train,
    y_train;
    xlim=(0, 1),
    xlabel="x",
    ylabel="y",
    title="posterior (VI with sparse grid)",
    label="Train Data",
)
scatter!(x_test, y_test; label="Test Data")
plot!(0:0.001:1, ap; label=false, ribbon_scale=2)
vline!(logistic.(opt.minimizer[3:end]); label="Pseudo-points")

# ## Exact Gaussian Process Inference
#
# Here we use Type-II MLE to train the hyperparameters of the Gaussian process.
# This means that our loss function is the negative log marginal likelihood.

# We re-calculate the log-likelihood of the test dataset with the
# default kernel parameters of value 1 for the sake of comparison.

logpdf(p_fx(x_test), y_test)

# We define a function which returns the negative log marginal
# likelihood for different variance and inverse lengthscale parameters
# of the Matern kernel and different pseudo-points. We ensure that the
# kernel parameters are positive with the softplus function
# ``f(x) = \log (1 + \exp x)``.

function loss_function(x, y)
    function negativelogmarginallikelihood(params)
        kernel =
            softplus(params[1]) * (Matern52Kernel() ∘ ScaleTransform(softplus(params[2])))
        f = GP(kernel)
        fx = f(x, 0.1)
        return -logpdf(fx, y)
    end
    return negativelogmarginallikelihood
end

#md nothing #hide

# We randomly initialize the kernel parameters, and minimize the
# negative log marginal likelihood with the LBFGS algorithm
# and obtain the following optimal parameters:

θ0 = randn(2)
opt = Optim.optimize(loss_function(x_train, y_train), θ0, LBFGS())

#-

opt.minimizer

# The optimized value of the variance is

softplus(opt.minimizer[1])

# and of the inverse lengthscale is

softplus(opt.minimizer[2])

# We compute the log-likelihood of the test data for the resulting optimized
# posterior. We can observe that there is a significant improvement over the
# log-likelihood with the default kernel parameters of value 1.

opt_kernel =
    softplus(opt.minimizer[1]) *
    (Matern52Kernel() ∘ ScaleTransform(softplus(opt.minimizer[2])))

opt_f = GP(opt_kernel)
opt_fx = opt_f(x_train, 0.1)
opt_p_fx = posterior(opt_fx, y_train)
logpdf(opt_p_fx(x_test), y_test)

# We visualize the posterior with optimized parameters.

scatter(
    x_train,
    y_train;
    xlim=(0, 1),
    xlabel="x",
    ylabel="y",
    title="posterior (optimized parameters)",
    label="Train Data",
)
scatter!(x_test, y_test; label="Test Data")
plot!(0:0.001:1, opt_p_fx; label=false, ribbon_scale=2)

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
"""
    absolute_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)

Absolute distance between two optimality cuts

"""
function absolute_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)
    c̃₁ = vcat(c₁.δQ, c₁.q) ./ num_subproblems(c₁)
    c̃₂ = vcat(c₂.δQ, c₂.q) ./ num_subproblems(c₂)
    return norm(c̃₁-c̃₂) / max(norm(c̃₁), norm(c̃₂))
end
function absolute_distance(c₁::AnyOptimalityCut{T}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return Inf
end
function absolute_distance(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::AnyOptimalityCut{T}) where T <: AbstractFloat
    return Inf
end

"""
    angular_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)

Angular distance between two optimality cuts

"""
function angular_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)
    return 1-abs(c₁.δQ ⋅ c₂.δQ) / (norm(c₁.δQ) * norm(c₂.δQ))
end
function angular_distance(c₁::AnyOptimalityCut{T}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return Inf
end
function angular_distance(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::AnyOptimalityCut{T}) where T <: AbstractFloat
    return Inf
end

"""
    spatioangular_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)

Spatioangular distance between two optimality cuts.

"""
function spatioangular_distance(c₁::AnyOptimalityCut, c₂::AnyOptimalityCut)
    return 1-abs(c₁.δQ ⋅ c₂.δQ) / (norm(c₁.δQ) * norm(c₂.δQ)) +
        (abs(c₁.q / num_subproblems(c₁) - c₂.q / num_subproblems(c₂)) /
        max(abs(c₁.q / num_subproblems(c₁)), abs(c₂.q / num_subproblems(c₂))))
end
function spatioangular_distance(c₁::AnyOptimalityCut{T}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return Inf
end
function spatioangular_distance(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::AnyOptimalityCut{T}) where T <: AbstractFloat
    return Inf
end


text:
module SoftmaxTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: VBSoftmaxOut, VBSoftmaxIn1, VBSoftmaxXi, VBSoftmaxA

#-------------
# Update rules
#-------------

@testset "VBSoftmaxOut" begin
    @test VBSoftmaxOut <: NaiveVariationalRule{Softmax}
    @test outboundType(VBSoftmaxOut) == Message{Categorical}
    @test isApplicable(VBSoftmaxOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBSoftmaxOut(nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[-1.0, 1.0], v=diageye(2)), ProbabilityDistribution(Multivariate, PointMass, m=[1.0, 2.0]), ProbabilityDistribution(Univariate, PointMass, m=1.0)) == Message(Univariate, Categorical, p=[0.11920292202211755, 0.8807970779778824])
end

@testset "VBSoftmaxIn1" begin
    @test VBSoftmaxIn1 <: NaiveVariationalRule{Softmax}
    @test outboundType(VBSoftmaxIn1) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(VBSoftmaxIn1, [ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBSoftmaxIn1(ProbabilityDistribution(Univariate, Categorical, p=[0.8, 0.2]), nothing, ProbabilityDistribution(Multivariate, PointMass, m=[1.0, 2.0]), ProbabilityDistribution(Univariate, PointMass, m=1.0)) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[0.5097033519265977, -0.039212315743835824], w=[0.2310585786300049 0.0; 0.0 0.19039853898894119])
end

@testset "VBSoftmaxXi" begin
    @test VBSoftmaxXi <: NaiveVariationalRule{Softmax}
    @test outboundType(VBSoftmaxXi) == Message{Function}
    @test isApplicable(VBSoftmaxXi, [ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution])

    @test ruleVBSoftmaxXi(ProbabilityDistribution(Univariate, Categorical, p=[0.8, 0.2]), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[-1.0, 1.0], v=diageye(2)), nothing, ProbabilityDistribution(Univariate, PointMass, m=1.0)) == Message(Multivariate, Function, mode=[2.23606797749979, 1.0])
end

@testset "VBSoftmaxA" begin
    @test VBSoftmaxA <: NaiveVariationalRule{Softmax}
    @test outboundType(VBSoftmaxA) == Message{Function}
    @test isApplicable(VBSoftmaxA, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBSoftmaxA(ProbabilityDistribution(Univariate, Bernoulli, p=[0.8, 0.2]), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[-1.0, 1.0], v=diageye(2)), ProbabilityDistribution(Multivariate, PointMass, m=[1.0, 2.0]), nothing) == Message(Univariate, Function, mode=-0.09647491510115125)
end

@testset "averageEnergy" begin
    @test averageEnergy(Softmax, ProbabilityDistribution(Univariate, Categorical, p=[0.8, 0.2]), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[-1.0, 1.0], v=diageye(2)), ProbabilityDistribution(Multivariate, PointMass, m=[1.0, 2.0]), ProbabilityDistribution(Univariate, PointMass, m=1.0)) == 2.907107586326735
end

end # module


text:
push(xs::StaticArrays.StaticVector, x) = vcat(xs, StaticArrays.SVector(x))
pushfirst(xs::StaticArrays.StaticVector, x) = StaticArrays.pushfirst(xs, x)
pushfirst(xs::StaticArrays.StaticVector, y, ys...) =
    foldl(pushfirst, (reverse(ys)..., y), init=xs)
pop(xs::StaticArrays.StaticVector) = (StaticArrays.pop(xs), xs[end])
popfirst(xs::StaticArrays.StaticVector) = (StaticArrays.popfirst(xs), xs[1])
deleteat(xs::StaticArrays.StaticVector, i) = StaticArrays.deleteat(xs, i)
_empty(xs::StaticArrays.StaticVector) =
    foldl(ntuple(identity, length(xs)); init=xs) do xs, _
        StaticArrays.pop(xs)
    end
_setindex(xs::StaticArrays.StaticArray, v, I...) = Base.setindex(xs, v, I...)

StaticArrays.SArray(x::SingletonVector) = StaticArrays.SVector{1}(x)
StaticArrays.SVector(x::SingletonVector) = StaticArrays.SVector{1}(x)
StaticArrays.MArray(x::SingletonVector) = StaticArrays.MVector{1}(x)
StaticArrays.MVector(x::SingletonVector) = StaticArrays.MVector{1}(x)


text:
module MsgPackRpcClientSocks

include("const.jl")

type Socks
  pool :: Array
end

function new()
  Socks(Union(Base.TcpSocket, Base.UdpSocket)[])
end

function show(self::Socks)
  println(self.pool)
  nothing
end

function is_empty(self::Socks)
  if length(self.pool) != 0
    return false
  else
    return true
  end
end

#function reject!(self::Socks)
#end

function get_sock(self, ptr::Int)
  self.pool[ptr]
end

function push!(self::Socks, sock::Union(Base.TcpSocket, Base.UdpSocket))
  enqueue!(self, sock)
end

function pop!(self; or_create = true, host::String = "localhost", port::Int = DEFAULT_PORT_NUMBER)
  try
    sock = pop!(self.pool)
    return sock
  catch
    # TODO: Add error handling
  end

  if or_create == true
    try
      return connect(host, port)
    catch
      # TODO: Add error handling
    end
  end

  nothing
end

function enqueue!(self::Socks, sock::Union(Base.TcpSocket, Base.UdpSocket))
  try
    Base.push!(self.pool, sock)
  catch
    self = Socks({})
  end
  self
end

function dequeue!(self::Socks)
  try
    sock = self.pool[1]
    deleteat!(self.pool, 1)
    return sock
  catch
    return nothing
  end
end

function connect_and_push!(self, host::String = "localhost", port::Int = DEFAULT_PORT_NUMBER)
  Base.push!(self.pool, connect(host, port))
  self
end

function connect_in_port_range_and_push!(self::Socks,
                                         host::String = "localhost",
                                         range::UnitRange = DEFAULT_PORT_NUMBER:DEFAULT_PORT_NUMBER)
  @sync begin
    @async for port = range
      connect_and_push!(self, host, port)
    end
  end
  self
end

function delete(self::Socks, sock::Base.TcpSocket)
  i = 1
  for x in self.pool
    if x == sock
      if isopen(self.pool[i])
        close(self.pool[i])
      end
      deleteat!(self.pool, i)
      break
    end
    i += 1
  end
  self
end

function destroy(self::Socks)
  try
    if is_empty(self)
      return true
    end
    for x in self.pool
      if x == nothing
        continue
      end
      if isopen(x)
        close(x)
      end
    end
  catch
    return false
  end
  true
end

end # module MsgPackRpcClientSocks

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    ZerosBridge{T} <: Bridges.Variable.AbstractBridge

Transforms constrained variables in [`MathOptInterface.Zeros`](@ref) to zeros,
which ends up creating no variables in the underlying model.

The bridged variables are therefore similar to parameters with zero values.
Parameters with non-zero value can be created with constrained variables in
[`MOI.EqualTo`](@ref) by combining a [`VectorizeBridge`](@ref) and this bridge.
The functions cannot be unbridged, given a function, we cannot determine, if
the bridged variables were used.

The dual values cannot be determined by the bridge but they can be determined
by the bridged optimizer using [`MathOptInterface.Utilities.get_fallback`](@ref)
if a `CachingOptimizer` is used (since `ConstraintFunction` cannot be got
as functions cannot be unbridged).
"""
struct ZerosBridge{T} <: AbstractBridge
    n::Int # Number of variables
end

const Zeros{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{ZerosBridge{T},OT}

function bridge_constrained_variable(
    ::Type{ZerosBridge{T}},
    model::MOI.ModelLike,
    set::MOI.Zeros,
) where {T}
    return ZerosBridge{T}(MOI.dimension(set))
end

function supports_constrained_variable(::Type{<:ZerosBridge}, ::Type{MOI.Zeros})
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:ZerosBridge})
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{<:ZerosBridge})
    return Tuple{Type,Type}[]
end

# Attributes, Bridge acting as a model
MOI.get(::ZerosBridge, ::MOI.NumberOfVariables)::Int64 = 0

function MOI.get(bridge::ZerosBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# References
MOI.delete(::MOI.ModelLike, ::ZerosBridge) = nothing

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::ZerosBridge)
    return MOI.Zeros(bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimal,
    bridge::ZerosBridge{T},
) where {T}
    return zeros(T, bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintDual,
    ::ZerosBridge{T},
) where {T}
    return error(
        "Unable to query the dual of a variable bound that was reformulated " *
        "using `ZerosBridge`. This usually arises in conic solvers when a " *
        "variable is fixed to a value. As a work-around, instead of creating " *
        "a fixed variable using variable bounds like `p == 1`, add an affine " *
        "equality constraint like `1 * p == 1` (or `[1 * p - 1,] in Zeros(1)`).",
    )
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.VariablePrimal,
    ::ZerosBridge{T},
    ::MOIB.IndexInVector,
) where {T}
    return zero(T)
end

function MOIB.bridged_function(::ZerosBridge{T}, ::MOIB.IndexInVector) where {T}
    return zero(MOI.ScalarAffineFunction{T})
end

function unbridged_map(::ZerosBridge, ::MOI.VariableIndex, ::MOIB.IndexInVector)
    return nothing
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Test, Qsosed

@testset "Disk UV fraction" begin
    bb = BlackBody(1e4)
    @test spectral_band_fraction_frequency(bb, 1e12, 1e18) ≈ 1 rtol = 1e-8
    @test spectral_band_fraction(bb, 1e-5, 1e2) ≈ 1 rtol = 1e-6

    bh = BlackHole(1e8 * M_SUN, 0.5, 0)
    @test disk_uv_fraction(bh, 1.0) ≈ 0
    uvf1 = disk_uv_fraction(bh, 6.01)
    uvf2 = disk_uv_fraction(bh, 30)
    uvf3 = disk_uv_fraction(bh, 200)
    @test uvf1 < uvf2
    @test uvf3 < uvf2
end

@testset "Radial UV fractions" begin
    bh = BlackHole(1e8 * M_SUN, 0.5, 0.0)
    corona =
        Corona(bh, hard_xray_fraction = 0.02, electron_energy = 100, reprocessing = true)
    warm = Warm(corona, electron_energy = 0.2, photon_index = 2.5)
    rr, uvf = radial_uv_fraction(corona, warm, n_r = 500)
    @test length(rr) == length(uvf) == 500
    for (i, r) in enumerate(rr)
        if r <= corona.radius
            @test uvf[i] == 0.0
            continue
        end
        if r <= warm.radius
            @test uvf[i] < disk_uv_fraction(bh, r)
            continue
        end
        if r > warm.radius
            @test uvf[i] == disk_uv_fraction(bh, r)
        end
    end
end

@testset "X-Ray fraction" begin
    bh = BlackHole(1e8 * M_SUN, 0.5, 0.0)
    corona =
        Corona(bh, hard_xray_fraction = 0.02, electron_energy = 100, reprocessing = true)
    warm = Warm(corona, electron_energy = 0.2, photon_index = 2.5)
    xf = total_xray_fraction(corona, warm)
    @test 0.1 < xf < 0.2
end


text:
"""
# KM3NeT Data Definitions v1.3.1-8-gae7dfb2
https://git.km3net.de/common/km3net-dataformat
"""

module Daqdatatypes
  const DAQSUPERFRAME = 101
  const DAQSUMMARYFRAME = 201
  const DAQTIMESLICE = 1001
  const DAQTIMESLICEL0 = 1002
  const DAQTIMESLICEL1 = 1003
  const DAQTIMESLICEL2 = 1004
  const DAQTIMESLICESN = 1005
  const DAQSUMMARYSLICE = 2001
  const DAQEVENT = 10001
end