text:
import Base:iterate, values
export
PosteriorFactorization,
currentPosteriorFactorization,
setCurrentPosteriorFactorization

global current_posterior_factorization = nothing

mutable struct PosteriorFactorization
    graph::FactorGraph
    posterior_factors::Dict{Symbol, PosteriorFactor}

    # Bookkeeping for faster lookup during scheduling
    edge_to_posterior_factor::Dict{Edge, PosteriorFactor}
    node_edge_to_cluster::Dict{Tuple{FactorNode, Edge}, Cluster}
    deterministic_edges::Set{Edge}

    # Fields for assembly
    energy_counting_numbers::Dict{FactorNode, Int64}
    entropy_counting_numbers::Dict{Region, Int64}

    # Flag to ensure required quantities for FE evaluation are computed
    free_energy_flag::Bool
end

"""
Return currently active `PosteriorFactorization`. Create one if there is none.
"""
function currentPosteriorFactorization()
    current_posterior_factorization === nothing ? PosteriorFactorization() : current_posterior_factorization
end

function setCurrentPosteriorFactorization(pfz::PosteriorFactorization)
    global current_posterior_factorization = pfz
end

"""
Initialize an empty `PosteriorFactorization` for sequential construction
"""
function PosteriorFactorization()
    setCurrentPosteriorFactorization(
        PosteriorFactorization(
            currentGraph(),
            Dict{Symbol, PosteriorFactor}(),
            Dict{Edge, PosteriorFactor}(),
            Dict{Tuple{FactorNode, Edge}, Symbol}(),
            deterministicEdges(currentGraph()),
            Dict{FactorNode, Int64}(),
            Dict{Region, Int64}(),
            false
        )
    )
end

"""
Initialize a `PosteriorFactorization` consisting of a single `PosteriorFactor` for the entire graph
"""
function PosteriorFactorization(fg::FactorGraph)
    pfz = setCurrentPosteriorFactorization(
        PosteriorFactorization(
            fg,
            Dict{Symbol, PosteriorFactor}(),
            Dict{Edge, PosteriorFactor}(),
            Dict{Tuple{FactorNode, Edge}, Symbol}(),
            deterministicEdges(fg),
            Dict{FactorNode, Int64}(),
            Dict{Region, Int64}(),
            false
        )
    )
    PosteriorFactor(fg, pfz=pfz, id=Symbol(""))
    return pfz
end

"""
Construct a `PosteriorFactorization` consisting of one `PosteriorFactor` for each argument
"""
function PosteriorFactorization(args::Vararg{Union{T, Set{T}, Vector{T}} where T<:Variable}; ids=Symbol[])
    pfz = PosteriorFactorization()
    isempty(ids) || (length(ids) == length(args)) || error("Length of ids must match length of posterior factor arguments")
    for (i, arg) in enumerate(args)
        if isempty(ids)
            PosteriorFactor(arg, id=generateId(PosteriorFactor))
        else
            PosteriorFactor(arg, id=ids[i])
        end
    end
    return pfz
end

iterate(pfz::PosteriorFactorization) = iterate(pfz.posterior_factors)
iterate(pfz::PosteriorFactorization, state) = iterate(pfz.posterior_factors, state)
values(pfz::PosteriorFactorization) = values(pfz.posterior_factors)

"""
Populate the target regions of the PosteriorFactor.
The targets depend on the variables of interest, the local posterior
factorization, and whether free energy will be evaluated. At the same
time, fields for fast lookup during scheduling are populated in the
posterior factorization.
"""
function setTargets!(pf::PosteriorFactor, pfz::PosteriorFactorization; target_variables=Set{Variable}(), free_energy=false, external_targets=false)
    # Initialize empty set of target cluster node and edges.
    # We cannot build a Set of Clusters directly, because duplicate Clusters are not removed.
    large_regions = Set{Tuple}()

    # Register quantities of interest with the current posterior factor
    for variable in target_variables
        edge = first(variable.edges)
        if edge in pf.internal_edges
            push!(pf.target_variables, variable)
        end
    end

    # Determine which target regions are required by external posterior factors
    if external_targets
        nodes_connected_to_external_edges = nodesConnectedToExternalEdges(pf)
        for node in nodes_connected_to_external_edges
            isa(node, DeltaFactor) && continue # Skip deterministic nodes

            target_edges = localStochasticInternalEdges(node, pf, pfz) # Find internal edges connected to node
            if length(target_edges) == 1 # Only one internal edge, the marginal for a single Variable is required
                push!(pf.target_variables, target_edges[1].variable)
            elseif length(target_edges) > 1 # Multiple internal edges, register the region for computing the joint marginal
                push!(large_regions, (node, target_edges))
            end
        end
    end

    # Determine which targets are required for evaluating the free energy
    if free_energy
        # Initialize counting numbers
        variable_counting_numbers = Dict{Variable, Number}()
        cluster_counting_numbers = Dict{Tuple, Number}()

        # Iterate over large regions
        nodes_connected_to_internal_edges = nodes(pf.internal_edges)
        for node in nodes_connected_to_internal_edges
            target_edges = localStochasticInternalEdges(node, pf, pfz) # Find stochastic internal edges connected to node
            if isDeltaConstraint(node)
                continue
            elseif !isa(node, DeltaFactor) # Node is stochastic
                if length(target_edges) == 1 # Single internal edge
                    increase!(variable_counting_numbers, target_edges[1].variable, Inf) # For average energy evaluation, make sure to include the edge variable
                elseif length(target_edges) > 1 # Multiple internal edges
                    region = (node, target_edges) # Node and edges for region that will later define the Cluster
                    increase!(cluster_counting_numbers, region, Inf) # Make sure to include the region for average energy evaluation
                end
            elseif isa(node, Equality)
                increase!(variable_counting_numbers, node.i[1].edge.variable, 1) # Increase the counting number for the equality-constrained variable
            else # Node is deterministic and not Equality, requires a joint belief for the inbound variables
                outbound_edge = node.interfaces[1].edge
                inbound_target_edges = filter(edge->edge!=outbound_edge, target_edges) # Collect local internal stochastic inbound edges
                if length(inbound_target_edges) == 1 # Just a single local stochastic inbound edge
                    increase!(variable_counting_numbers, inbound_target_edges[1].variable, 1) # Increase counting number for variable on inbound edge
                elseif length(inbound_target_edges) >= 2 # Multiple stochastic inbound edges require the formation of a large (joint) region
                    region = (node, inbound_target_edges) # Region for node and stochastic inbound edges
                    increase!(cluster_counting_numbers, region, 1) # Increase the counting number for that region
                end
            end
        end

        # Iterate over small regions
        for edge in pf.internal_edges
            increase!(variable_counting_numbers, edge.variable, -(degree(edge) - 1)) # Discount single variables
            if edge in pfz.deterministic_edges
                variable_counting_numbers[edge.variable] = 0 # However, ignore entropies of deterministic beliefs
            end
        end

        # All targets with a non-zero counting number are required for free energy evaluation
        for (variable, cnt) in variable_counting_numbers
            if cnt != 0
                push!(pf.target_variables, variable)
            end
        end
        for (region, cnt) in cluster_counting_numbers
            if cnt != 0
                push!(large_regions, region)
            end
        end
    end

    # Determine breaker initialization for EP sites
    for ep_site in pf.ep_sites
        push!(pf.breaker_interfaces, ultimatePartner(ep_site))
    end

    # Determine which interfaces require breakers based on node-specific properties
    for edge in pf.internal_edges
        requiresBreaker(edge.a) && push!(pf.breaker_interfaces, edge.a)
        requiresBreaker(edge.b) && push!(pf.breaker_interfaces, edge.b)

        # Register internal edges with the posterior factorization for fast lookup during scheduling
        pfz.edge_to_posterior_factor[edge] = pf
    end

    # Create clusters, and register clusters with the posterior factorization for fast lookup during scheduling
    for region in large_regions
        cluster = Cluster(region...) # Use the region definition to construct a Cluster
        push!(pf.target_clusters, cluster)
        for edge in cluster.edges
            pfz.node_edge_to_cluster[(cluster.node, edge)] = cluster
        end
    end

    # Register whether posterior factorization is prepared for free energy evaluation
    pfz.free_energy_flag = free_energy

    return pf
end

"""
Return the local stochastic regions around `node`
"""
function localStochasticRegions(node::FactorNode, pfz::PosteriorFactorization)
    regions = Region[]
    for interface in node.interfaces
        partner = ultimatePartner(interface)
        if !(interface.edge in pfz.deterministic_edges) # If edge is stochastic
            push!(regions, region(node, interface.edge))
        end
    end

    return regions
end

"""
Find stochastic edges that are internal to the posterior factor and connected to node.
This function is used for constructing clusters. Therefore, the edges are returned
in the same order as the node's interfaces.
"""
function localStochasticInternalEdges(node::FactorNode, pf::PosteriorFactor, pfz::PosteriorFactorization)
    local_internal_edges = Edge[]
    for interface in node.interfaces
        if (interface.edge in pf.internal_edges) && !(interface.edge in pfz.deterministic_edges) # Edge is internal to posterior factor and stochastic
            push!(local_internal_edges, interface.edge) # Otherwise include the edge
        end
    end

    return local_internal_edges
end


text:
struct Rackett{MW,T} <: SatLiquidVolumeModel
    mw::MW
    tc::T
    pc::T
    zc::T

end

Rackett(;tc,pc,zc,mw=nothing) = Rackett(tc,pc,zc,mw)

function Rackett(model::ThermoModel)
    mw = molecular_weight(model)
    tc = temperature(model,CriticalPoint())
    pc = pressure(model,CriticalPoint())
    zc = compressibility_factor(model,CriticalPoint())
    return Rackett(tc,pc,zc,mw)
end

function mol_volume_impl(::SingleSatT,model::Rackett,t::T) where T
    _2_7 = T(2//7)
    _1 = T(1.0)
    tc = only(model.tc)
    zc = only(model.zc)
    pc = only(model.pc)
    RGAS*tc/pc*zc^(_1 + (_1 - T/Tc)^(_2_7))
end

function mol_volume_impl(::MultiSatT,model::Rackett,t::T,x) where T
    tc = model.tc
    zc = model.zc
    pc = model.pc
    _2_7 = T(2//7)
    _1 = T(1.0)
    zmix = dot(zc,x)
    tmix = dot(tc,x)
    mw = dot(MW,x)
    tr = t/tmix
    RGAS*tc/pc*zmix^(_1 + (_1 - tr)^(_2_7))
end

ext: jl


text:
function solve(U, queries)

end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        U = parse(Int, readline(stdin))
        queries = Vector{String}(missing, 10^4)
        for i in 1:10^4
            _, R = stdin |> readline |> split
            queries[i] = R
        end

        solve(case, U, queries)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Linear programming using Julia

# - Let's consider a function of two variables $f=f(x,y)$ defined by
#
# $$
#     f(x,y) = 5x + 3y.
# $$
#
# - We would like find $(x,y)=(\hat{x}, \hat{y})$ that maximize the value of f(x,y) with the following contstraits, namely:
#
# $$
# \begin{aligned}
# \text{maximize} &\ f(x,y) \\
# \text{subject to}:&\ 0\leq x \leq 2, \\
#                   &\ 0\leq y \leq 30, \\
#                   &\ 5x + 3y \leq 3.
# \end{aligned}
# $$
#
# - Can we solve this problem using Julia ? Let's find out.

# - JuMP is a domain-specific modeling language for mathematical optimization embedded in Julia.
# - Tulip is an open-source interior-point solver for linear optimization, written in pure Julia.

# ## Load Packages

using JuMP
import Tulip
using LazySets
using Plots
using LaTeXStrings
using Colors

# Models are created with the Model function. The following example adopts `Tulip` as a solver. Note that the version of JuMP should be $\geq 0.21$

]st JuMP

model = Model(Tulip.Optimizer)

# The following commands will create two variables, `x` and `y`, with both lower and upper bounds.

xmin, xmax = 0, 2
ymin, ymax = 0, 30
# you can use `≤` instead of `<=`
@variable(model, xmin ≤ x ≤ xmax)
@variable(model, ymin ≤ y ≤ ymax);

# Next we'll set our objective. Note again the model, so we know which model's objective we are setting! The objective sense, `Max` or `Min`, should be provided as the second argument. Note also that we don't have a multiplication * symbol between 5 and our variable `x` - Julia is smart enough to not need it!

f(x,y)=5x+3y
@objective(model, Max, 5x+3y)

# Adding constraints is a lot like setting the objective. Here we create a less-than-or-equal-to constraint using `<=`, but we can also create equality constraints using `==` and greater-than-or-equal-to constraints with `>=`

@constraint(model, con, 1x + 5y <= 3)

# ## Summarize model

# Let's see the summary of our `model`. The result should be equivalent to the condition:
#
# $$
# \begin{aligned}
# \text{maximize} &\ f(x,y) \\
# \text{subject to}:&\ 0\leq x \leq 2, \\
#                   &\ 0\leq y \leq 30, \\
#                   &\ 5x + 3y \leq 3.
# \end{aligned}
# $$

model

# - It's time to solve the model using `JuMP.optimize!`:

optimize!(model)

# Done!

# ## Output the result

# We can access the solution via `objective_value` and `value` function.

x̂ = value(x)
ŷ = value(y)
@assert f̂ = objective_value(model) ≈ f(x̂, ŷ)
@show x̂, ŷ

# # Visualize result

# - Let's visualize the constraint $x+5y\leq 3$. It can be written as $\langle[1,5],[x,y]\rangle\leq 3$ where $\langle\bullet,\bullet\rangle$ stands for inner product of 2-dimensional Euclidian space. LazySets provides a recipe that visualize the half plane of the form $\langle a,x\rangle \leq b$ via `LazySets.HalfSpace(a,b)`.
# - Note that the version of LazySets should be >= 0.35

p = plot(xlims=[xmin, xmax], ylims=[ymin, ymax])
hs = HalfSpace([1.,5.], 3.)
p = plot!(p, hs)
hp = Hyperplane([1.,5.], 3.)
p = plot!(p, hp,color=:red)

# Where is the solution, say, $(\hat{x},\hat{y})$ ?

scatter!(p, [x̂], [ŷ])

# # Linear programming with multiple constraits

# - Again, let's consider a function of two variables $f(x,y)=5x+3y$.
# - We would like find $(x,y)=(\hat{x}, \hat{y})$ that maximize the value of $f(x,y)$ with the following multiple linear contstraits, namely:
#
# $$
# \begin{aligned}
# \text{maximize} &\ f(x,y)=3x+5y, \\
# \text{subject to}:&\
#     \begin{bmatrix}
#         2 & 2 \\
#         2 & -4 \\
#         -2 & 1 \\
#         0 & -1 \\
#         0 & 1
#     \end{bmatrix}
#     \begin{bmatrix} x \\ y \end{bmatrix}
#     \leq
#     \begin{bmatrix}
#         33 \\
#         8 \\
#         -5 \\
#         -1 \\
#         8
#     \end{bmatrix}.
# \end{aligned}
# $$
#
# - This condition is taken from [MML book](https://mml-book.github.io/) p240 Example 7.5. Note that the original example has a typo. See https://github.com/mml-book/mml-book.github.io/issues/547
# - Can we solve this problem using Julia ? Of course.

# ## Define model and constraits

# +
f(x,y)=5x+3y
model = Model(Tulip.Optimizer)

@variable(model, x)
@variable(model, y)

@objective(model, Max, f(x,y))

X=Float64[
    2   2
    2  -4
    -2  1
    0  -1
    0   1
]

bs = Float64[
    33
    8
    -5 # MML book has typo. I've fixed from 5 to -5
    -1
    8
]

@constraint(model, con, X * [x,y] .≤ bs)
# -

# ## Optimize model

optimize!(model)

f̂=objective_value(model)
x̂=value(x)
ŷ=value(y)
@show x̂, ŷ

# ## Visualize result

# +
p = plot(xlim=[0,16], ylim=[0,10], aspect_ratio=:equal, ticks=0:2:16)
hs = HalfSpace{Float64, Vector{Float64}}[]
hp = Hyperplane{Float64, Vector{Float64}}[]

cs = distinguishable_colors(length(bs))
for (i, b) in zip(1:size(X)[1], bs)
    push!(hs, HalfSpace(X[i, 1:2], b))
    push!(hp, Hyperplane(X[i,1:2], b))
end

ph = HPolyhedron(hs)
plot!(p, convert(HPolytope, ph))
plot!.(Ref(p), hp)
scatter!(p, [x̂], [ŷ], label=L"(\hat{x}, \hat{y})", marker=:star, msize=8)
xs = 0:0.01:16
ys = 0:0.01:10
contour!(p, xs, ys, f)
# -

# # References:
#
# - JuMP quick start: https://jump.dev/JuMP.jl/v0.21.3/quickstart/


text:
mutable struct Style
    color::Colorant
    linewidth::Float64
    fontface::String
    fontsize::Float64
    linecap::Symbol
    linejoin::Symbol
    dashpattern::Array
    background::Colorant
    function Style(color, linewidth, fontface, fontsize, linecap, linejoin, dashpattern, background)
        new(color, linewidth, fontface, fontsize, linecap, linejoin, dashpattern, background)
    end
end

# blank
Style() = Style(colorant"black", 1.0, "JuliaMono-Regular", 20,  :butt, :miter, [5, 0], colorant"white")

# color
Style(col::Colorant) = Style(col,  1.0, "JuliaMono-Regular", 20,  :butt, :miter, [5, 0], colorant"white")

# linewidth
Style(lw::Float64) = Style(colorant"black",  lw, "JuliaMono-Regular", 20,  :butt, :miter, [5, 0], colorant"white")

# color linewidth
Style(col::Colorant, lw::Float64) = Style(col, lw, "JuliaMono-Regular", 20,  :butt, :miter, [5, 0], colorant"white")

# font fontsize
Style(ff::String, fs) = Style(colorant"black",  1.0, ff, fs, :butt, :miter, [5, 0], colorant"white")

# color font fontsize
Style(col::Colorant, ff::String, fs) = Style(col, 1.0, ff, fs, :butt, :miter, [5, 0], colorant"white")

function applystyle(style::Style)
    setcolor(style.color)
    setline(style.linewidth)
    fontface(style.fontface)
    fontsize(style.fontsize)
    setlinecap(style.linecap)
    setlinejoin(style.linejoin)
    setdash(style.dashpattern)
end


text:
function gnlse(pulse::Pulse, fiber::Fiber, betas::Array{Any,1}; dispersion::Bool = true,  fr::Float64=0.0, shock::Bool=false, plasma::Bool=false, nsaves::Int64=1000)
    ω0 = 2pi*c/pulse.λ0
    γ = getGamma(fiber.gas, fiber.p, fiber.Tk, ω0, fiber.diameter)

    # grid
    τ = pulse.τ
    dτ = τ[2] - τ[1]
    convt = length(τ)*dτ
    W = 2pi.*fftshift(fftfreq(length(τ), 1/dτ))
    dW = W[2] - W[1]
    convw = dW/(2pi)

    # define the indices for where the signal after them ignored
    # the limits here are 200 nm - 2000 nm
    # downlim = findfirst(x -> x>(2pi*3e8/downlim - ω0), W)
    # uplim = findfirst(x -> x>(2pi*3e8/uplim - ω0), W)

    #dispersion
    if dispersion == true
        if length(betas) > 0

            β2 = wbeta2(ω0, fiber.gas, fiber.p, fiber.Tk, fiber.diameter/2)
            betas0 = [β2]
            β = zeros(Float64, size(W))
            for i in 1:length(betas0)
                @. β = β + betas0[i]/factorial(i+1)*(W)^(i+1)
            end
            L = 1im*β
        else
            L = @. 1im*L_operator.(W,ω0, fiber.gas, fiber.p, fiber.Tk, fiber.diameter/2) - fiber.α/2
        end

        # L = []
        # arr = readdlm("L_operator.txt")
        # for i in 1:length(arr)
        #     append!(L, parse(Complex{Float64}, arr[i]))
        # end
    else
    # if dispersion == true
    #     β = zeros(Float64, size(W))
    #     for i in 1:length(betas)
    #         @. β = β + betas[i]/factorial(i+1)*(W)^(i+1)
    #     end
    #     L = 1im*β
    # else
        L = 0
    end
    #Raman
    τ1 = 62.5e-15
    τ2 = 77e-15
    RW = sdo(τ, τ1, τ2)

    #windows
    gw(z) = (convt) .* exp.(-0.5 .*((W.+ 2e15 .-  maximum(W)./2)./(maximum(W).*0.5)).^30).* exp.(-z .* L)
    gt(z) =  exp.(-0.5 .*((τ)./(maximum(τ).*0.5)).^30) .* convw

    Uw0 = shifted_ifft(pulse.data) .* convt
    zs = (0.0, fiber.length)
    Z = LinRange(0, fiber.length, nsaves)
    prob = ODEProblem(rhs, Uw0, zs, (L, Progress(floor(Int, Float64(fiber.length*1000)), 2),
                      W, τ, ω0, γ, RW, fiber, pulse, fr, shock, plasma, dW, dτ, convw, convt, gw, gt))
    sol = solve(prob, Vern7(), reltol=1e-4, abstol=1e-14, saveat=fiber.length/nsaves)
    res = zeros(Complex{Float64}, size(sol))
    @simd for i in eachindex(sol.t)
        @inbounds res[:,i] = Complex{Float64}.(sol[:,i] .* exp.(L .* sol.t[i]))
    end
    rest = zeros(Complex{Float64}, size(sol))
    @simd for i in eachindex(sol.t)
        @inbounds rest[:,i] = Complex{Float64}.(shifted_fft(res[:,i])).*convw
    end

    τ, W, rest, res, Z
end

function rhs(dUw, Uw, para, z)
    L = para[1]
    W, τ, ω0, γ, RW, fiber, pulse, fr, shock, plasma, dW, dτ, convw, convt, gw, gt = para[3:end]

    update!(para[2], max(floor(Int, Float64(z*1000)), 1))
    U = shifted_fft(Uw .* exp.(L .* z)) .* convw
    Uw = shifted_ifft(abs2.(U) .* U).*convt
    if fr != 0.0
        covRaman(U, RW, convt, convw, fr)
    end
    if shock == true
        @. Uw *= (1.0 .+ W./ω0)
    end

    if plasma == true
        Et = getElectricField(pulse, fiber)
        #Pl = plasmaTerms(Et, τ, U, fiber.Ip, fiber.p, pulse.λ0, fiber.diameter).* (length(τ)*dτ)
        Pl = geisslerPlasma(Et, τ, U, fiber.Ip, fiber.p, pulse.λ0, fiber.diameter, convw, convt).* (length(τ)*dτ)
        #Pl = fedotovPlasma(Et, τ, U, fiber.Ip, fiber.p, pulse.λ0, fiber.diameter).* (length(τ)*dτ)
        # @. dUw = (1im*γ*uuw + Pl) * exp(-z * L)
        dU =  @. (1im*γ*Uw + Pl)
        dUt = shifted_fft(dU) .* gt(z)
        dUw .= shifted_ifft(dUt).* gw(z)
    else
        # @. dUw = (1im*γ*uuw) * exp(-z * L)
        dU = @. (1im*γ*Uw)
        dUt = shifted_fft(dU) .* gt(z)
        dUw .= shifted_ifft(dUt).* gw(z)
    end

end

#based on soliton order
function gnlse(order::Float64, τ::Array{Float64,1}, betas::Array{Float64,1}, λ0::Float64, FWHM::Float64, fiber::Fiber; dispersion::Bool = true,  fr::Float64=0.0, shock::Bool=false, plasma::Bool=false, nsaves=1000)
    τ0 = FWHM/1.76
    ω0 = 2pi*c/λ0
    γ = getGamma(fiber.gas, fiber.p, fiber.Tk, ω0, fiber.diameter)

    P0 = abs(betas[1])*(order^2)/(γ*τ0^2)
    energy = 2τ0*P0
    pulse = sechPulse(τ, energy, FWHM, λ0)

    # grid
    τ = pulse.τ
    dτ = τ[2] - τ[1]
    convt = length(τ)*dτ
    W = 2pi.*fftshift(fftfreq(length(τ), 1/dτ))
    dW = W[2] - W[1]
    convw = (dW/(2pi))

    #dispersion
    if dispersion == true
        β = zeros(Float64, size(W))
        for i in 1:length(betas)
            @. β = β + betas[i]/factorial(i+1)*(W)^(i+1)
        end
        L = 1im*β
    else
        L = 0
    end

    #Raman
    if fr != Nothing
        γ = γ/(1-fr)
        Ht = Htot(τ)
        Ht = f2(Ht, τ)
        RT = Ht
        RT = RT/trapz(τ, RT)
        RW = shifted_ifft(RT)*convt
    end
    function rhs(dUw, Uw, para, z)
        L = para[1]
        # println(z)
        update!(para[2], max(floor(Int, Float64(z*1000)), 1))
        U = shifted_fft(Uw .* exp.(L .* z)) .* convw
        uu = abs2.(U) .* U
        uuw = shifted_ifft(uu).*convt
        if fr != Nothing
            # println(fr)
            Iw = shifted_ifft(abs2.(U)).*convt
            M = shifted_fft(Iw.*RW) .* U .* convw
            Mw = shifted_ifft(M) .* convt
            @. uuw = (1.0 - fr) .* uuw .+ fr.* Mw
        end
        if shock == true
            @. uuw = uuw + uuw * W/ω0
        end

        if plasma == true
            Et = getElectricField(pulse, fiber)
            Pl = plasmaTerms(Et, τ, U, fiber.Ip, fiber.p, pulse.λ0, fiber.diameter).* (length(τ)*dτ)

            # savefig(plot(τ, abs2.(plasma)), "plas")
            @. dUw =  (1im*γ*uuw + Pl) * exp(-z * L)
        else
            @. dUw = (1im*γ*uuw) * exp(-z * L)
        end

    end
    Uw0 = shifted_ifft(pulse.data) .* convt
    zs = (0.0, fiber.length)

    prob = ODEProblem(rhs, Uw0, zs, (L, Progress(floor(Int, Float64(fiber.length*1000)), 2)))
    sol = solve(prob, Vern7(), reltol=1e-8, abstol=1e-16, saveat=fiber.length/nsaves)
    res = zeros(Complex{Float64}, size(sol))
    for i in eachindex(sol.t)
        res[:,i] = Complex{Float64}.(sol[:,i] .* exp.(L .* sol.t[i]))
    end
    rest = zeros(Complex{Float64}, size(sol))
    for i in eachindex(sol.t)
        rest[:,i] = Complex{Float64}.(shifted_fft(res[:,i])).*convw
    end

    τ, W.+ω0, rest, res
end

function myplot(τ, W, solt, sol, name)
    ldw = 10.0.*log10.(abs2.(sol) .+ 1e-27);

    ldt = abs2.(solt);

    savefig(heatmap(ldw, clim=(maximum(ldw)-50.0, maximum(ldw)), color = :pu_or), "$name_sol")

    savefig(heatmap(ldt, clim=(maximum(ldt)-70.0, maximum(ldt)), color = :pu_or), "$name_solt")
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
patch(x, x̄) = x .- x̄

function state(o, x)
  if isleaf(x)
    return init(o, x)
  else
    x, _ = functor(x)
    return map(x -> state(o, x), x)
  end
end

function _update(o, st, x, x̄s...)
  st, x̄ = apply(o, st, x, x̄s...)
  return st, patch(x, x̄)
end

function update(o, state, x::T, x̄s...) where T
  if all(isnothing, x̄s)
    return state, x
  elseif isleaf(x)
    return _update(o, state, x, x̄s...)
  else
    x̄s = map(x̄ -> functor(typeof(x), x̄)[1], x̄s)
    x, restructure = functor(typeof(x), x)
    xstate = map((state, x, x̄s...) -> update(o, state, x, x̄s...), state, x, x̄s...)
    return map(first, xstate), restructure(map(last, xstate))
  end
end

# default all rules to first order calls
apply(o, state, x, dx, dxs...) = apply(o, state, x, dx)


text:
function initialize_avr!(
    device_states,
    device::PSY.DynamicGenerator{M, S, PSY.AVRFixed, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #In AVRFixed, V_ref is used as Vf
    Vf = get_inner_vars(device)[Vf_var]
    PSY.set_V_ref!(PSY.get_avr(device), Vf)
    #Update Control Refs
    device.ext[CONTROL_REFS][V_ref_index] = Vf
end

function initialize_avr!(
    device_states,
    device::PSY.DynamicGenerator{M, S, PSY.AVRSimple, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #In AVRFixed, V_ref is used as Vf
    Vf0 = get_inner_vars(device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(get_inner_vars(device)[VR_gen_var]^2 + get_inner_vars(device)[VI_gen_var]^2)
    #Solve V_ref
    V_ref = Vm
    #Set Vf state equals to Vf0
    avr_ix = get_local_state_ix(device, PSY.AVRSimple)
    avr_states = @view device_states[avr_ix]
    avr_states[1] = Vf0 #δ
    #Set V_ref
    PSY.set_V_ref!(PSY.get_avr(device), V_ref)
    #Update Control Refs
    device.ext[CONTROL_REFS][V_ref_index] = V_ref
end

function initialize_avr!(
    device_states,
    device::PSY.DynamicGenerator{M, S, PSY.AVRTypeI, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(get_inner_vars(device)[VR_gen_var]^2 + get_inner_vars(device)[VI_gen_var]^2)

    #Get parameters
    avr = PSY.get_avr(device)
    Ka = PSY.get_Ka(avr)
    Ke = PSY.get_Ke(avr)
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)
    #Obtain saturated Vf
    Se_Vf = Ae * exp(Be * abs(Vf0))

    #States of AVRTypeI are Vf, Vr1, Vr2, Vm
    #To solve V_ref, Vr1, Vr2
    function f!(out, x)
        V_ref = x[1]
        Vr1 = x[2]
        Vr2 = x[3]

        out[1] = Vf0 * (Ke + Se_Vf) - Vr1 #16.12c
        out[2] = Ka * (V_ref - Vm - Vr2 - (Kf / Tf) * Vf0) - Vr1 #16.12a
        out[3] = (Kf / Tf) * Vf0 + Vr2 #16.12b
    end
    x0 = [1.0, Vf0, Vf0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        device.ext[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(device, PSY.AVRTypeI)
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vf0 #δ
        avr_states[2] = sol_x0[2] #ω
        avr_states[3] = sol_x0[3]
        avr_states[4] = Vm
    end
end

function initialize_avr!(
    device_states,
    device::PSY.DynamicGenerator{M, S, PSY.AVRTypeII, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(get_inner_vars(device)[VR_gen_var]^2 + get_inner_vars(device)[VI_gen_var]^2)

    #Get parameters
    avr = PSY.get_avr(device)
    K0 = PSY.get_K0(avr)
    T1 = PSY.get_T1(avr)
    T2 = PSY.get_T2(avr)
    T3 = PSY.get_T3(avr)
    T4 = PSY.get_T4(avr)
    Tr = PSY.get_Tr(avr)
    Vr_max = PSY.get_Vr_max(avr)
    Vr_min = PSY.get_Vr_min(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)
    #Obtain saturated Vf
    Se_Vf = Ae * (exp(Be * abs(Vf0)))

    #States of AVRTypeII are Vf, Vr1, Vr2, Vm
    #To solve V_ref, Vr1, Vr2
    function f!(out, x)
        V_ref = x[1]
        Vr1 = x[2]
        Vr2 = x[3]

        Vr = K0 * Vr2 + (T4 / T3) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) #16.21

        out[1] = (Vf0 * (1.0 + Se_Vf) - Vr) #16.18
        out[2] = K0 * (1.0 - (T2 / T1)) * (V_ref - Vm) - Vr1 #16.14
        out[3] = (1.0 - (T4 / T3)) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) - K0 * Vr2  #16.20
    end
    x0 = [1.0, Vf0, Vf0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        device.ext[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(device, PSY.AVRTypeII)
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vf0 #Vf
        avr_states[2] = sol_x0[2] #Vr1
        avr_states[3] = sol_x0[3] #Vr2
        avr_states[4] = Vm #Vm
    end
end


text:
function sub_quote(W::Int, T::Symbol, fast::Bool)::Expr
  vtyp = vtype(W, T)
  instrs = "%res = fneg $(fast_flags(fast)) $vtyp %0\nret $vtyp %res"
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T}}, data(v)))
  end
end

@generated vsub(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, JULIA_TYPES[T], false)
@generated vsub_fast(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, JULIA_TYPES[T], true)

@inline vsub(v) = -v
@inline vsub_fast(v) = Base.FastMath.sub_fast(v)
@inline vsub(v::Vec{<:Any,<:NativeTypes}) = vsub(zero(v), v)
@inline vsub_fast(v::Vec{<:Any,<:UnsignedHW}) = vsub(zero(v), v)
@inline vsub_fast(v::Vec{<:Any,<:NativeTypes}) = vsub_fast(zero(v), v)
@inline vsub(x::NativeTypes) = Base.FastMath.sub_fast(x)
@inline vsub_fast(x::NativeTypes) = Base.FastMath.sub_fast(x)

@inline vinv(v) = inv(v)
@inline vinv(v::AbstractSIMD{W,<:FloatingTypes}) where {W} = vfdiv(one(v), v)
@inline vinv(v::AbstractSIMD{W,<:IntegerTypesHW}) where {W} = inv(float(v))
@inline Base.FastMath.inv_fast(v::AbstractSIMD) = Base.FastMath.div_fast(one(v), v)

@inline vabs(v) = abs(v)
@inline vabs(v::AbstractSIMD{W,<:Unsigned}) where {W} = v
@inline vabs(v::AbstractSIMD{W,<:Signed}) where {W} = ifelse(v > 0, v, -v)

@inline vround(v) = round(v)
@inline vround(v::AbstractSIMD{W,<:Integer}) where {W} = v
@inline vround(v::AbstractSIMD{W,<:Integer}, ::RoundingMode) where {W} = v

function bswap_quote(W::Int, T::Symbol, st::Int)::Expr
  typ = 'i' * string(8st)
  suffix = 'v' * string(W) * typ
  vtyp = "<$W x $typ>"
  decl = "declare $(vtyp) @llvm.bswap.$(suffix)($(vtyp))"
  instrs = """
    %res = call $vtyp @llvm.bswap.$(suffix)($vtyp %0)
    ret $vtyp %res
  """
  ret_type = :(_Vec{$W,$T})
  llvmcall_expr(decl, instrs, ret_type, :(Tuple{$ret_type}), vtyp, [vtyp], [:(data(x))])
end
@generated Base.bswap(x::Vec{W,T}) where {T<:IntegerTypesHW,W} = bswap_quote(W, JULIA_TYPES[T], sizeof(T))
@inline Base.bswap(x::VecUnroll{<:Any,<:Any,<:IntegerTypesHW}) = VecUnroll(fmap(bswap, data(x)))
@inline Base.bswap(x::AbstractSIMDVector{<:Any,<:IntegerTypesHW}) = bswap(Vec(x))
@inline Base.bswap(x::AbstractSIMD{<:Any,Float16}) = reinterpret(Float16, bswap(reinterpret(UInt16, x)))
@inline Base.bswap(x::AbstractSIMD{<:Any,Float32}) = reinterpret(Float32, bswap(reinterpret(UInt32, x)))
@inline Base.bswap(x::AbstractSIMD{<:Any,Float64}) = reinterpret(Float64, bswap(reinterpret(UInt64, x)))


text:
    :(function (o::Creature.Organism,)
        function func_3(var_12::Int64=8428943236158514626,var_13::Int8=-102)
            return var_12
        end
        function func_2(var_10::Int64=8331567932446283170,var_11::Bool=true)
            return var_10
        end
        Creature.stepRight(o)
        local var_7::Int8 = 125
        Creature.stepUp(o)
        Creature.eatUp(o,Int(var_7))
        Creature.stepRight(o)
        Creature.stepRight(o)
        Creature.stepRight(o)
        Creature.eatLeft(o,Int(var_7))
        begin
            local var_8::Int8
            for var_8 = 1:var_7
                Creature.stepLeft(o)
                Creature.stepRight(o)
                local var_9::Int16 = -28033
            end
        end
        Creature.stepRight(o)
	)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
@testset "mysql" begin
    # Some mysql drivers don't support reading usernames and passwords from .odbc.ini
    # files so you have to specify those explicitly in the DSN setup call:
    # dsn = ODBC.DSN("MySQL-test", "root", "")
    dsn = ODBC.DSN("MySQL-test")

    @testset "basic queries" begin
        dbs = ODBC.query(dsn, "show databases")
        ODBC.query(dsn, "use mysql")
        data = ODBC.query(dsn, "select table_name from information_schema.tables")
    end

    @testset "create testdb" begin
        ODBC.execute!(dsn, "drop database if exists testdb")
        ODBC.execute!(dsn, "create database testdb")
        ODBC.execute!(dsn, "use testdb")
    end

    @testset "test1" begin
        ODBC.execute!(dsn, "drop table if exists test1")
        ODBC.execute!(dsn, "create table test1
                            (test_bigint bigint,
                             test_bit bit,
                             test_decimal decimal,
                             test_int int,
                             test_numeric numeric,
                             test_smallint smallint,
                             test_mediumint mediumint,
                             test_tiny_int tinyint,
                             test_float float,
                             test_real double,
                             test_date date,
                             test_datetime datetime,
                             test_timestamp timestamp,
                             test_time time,
                             test_year year,
                             test_char char(1),
                             test_varchar varchar(16),
                             test_binary binary(2),
                             test_varbinary varbinary(16),
                             test_tinyblob tinyblob,
                             test_blob blob,
                             test_mediumblob mediumblob,
                             test_longblob longblob,
                             test_tinytext tinytext,
                             test_text text,
                             test_mediumtext mediumtext,
                             test_longtext longtext
                            )")
        data = ODBC.Query(dsn, "select * from information_schema.columns where table_name = 'test1'")
        ODBC.execute!(dsn, "insert test1 VALUES
                            (1, -- bigint
                             1, -- bit
                             1.0, -- decimal
                             1, -- int
                             1.0, -- numeric
                             1, -- smallint
                             1, -- mediumint
                             1, -- tinyint
                             1.2, -- float
                             1.2, -- double
                             '2016-01-01', -- date
                             '2016-01-01 01:01:01', -- datetime
                             '2016-01-01 01:01:01', -- timestamp
                             '01:01:01', -- time
                             2016, -- year
                             'A', -- char(1)
                             'hey there sailor', -- varchar
                             cast('12' as binary(2)), -- binary
                             NULL, -- varbinary
                             'hey there abraham', -- tinyblob
                             'hey there bill', -- blob
                             'hey there charlie', -- mediumblob
                             'hey there dan', -- longblob
                             'hey there ephraim', -- tinytext
                             'hey there frank', -- text
                             'hey there george', -- mediumtext
                             'hey there hank' -- longtext
                            )")
        data = DataFrame(ODBC.Query(dsn, "select * from test1"))

        @test size(data) == (1,27)
        @test Tables.schema(data).types == (
            Union{Int64, Missing},
            Union{Int8, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int32, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int16, Missing},
            Union{Int32, Missing},
            Union{Int8, Missing},
            Union{Float32, Missing},
            Union{Float64, Missing},
            Union{ODBC.API.SQLDate, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTime, Missing},
            Union{Int16, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{String, Missing}
        )
        @test data[1][1] === Int64(1)
        @test data[2][1] === Int8(1)
        @test data[3][1] === DecFP.Dec64(1)
        @test data[4][1] === Int32(1)
        @test data[5][1] === DecFP.Dec64(1)
        @test data[6][1] === Int16(1)
        @test data[7][1] === Int32(1)
        @test data[8][1] === Int8(1)
        @test data[9][1] === Float32(1.2)
        @test data[10][1] === Float64(1.2)
        @test data[11][1] === ODBC.API.SQLDate(2016,1,1)
        @test data[12][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
        @test data[13][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
        @test data[14][1] === ODBC.API.SQLTime(1,1,1)
        @test data[15][1] === Int16(2016)
        @test string(data[16][1]) == "A"
        @test string(data[17][1]) == "hey there sailor"
        @test data[18][1] == UInt8[0x31,0x32]
        @test ismissing(data[19][1])
        @test data[20][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x61,0x62,0x72,0x61,0x68,0x61,0x6d]
        @test data[21][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x62,0x69,0x6c,0x6c]
        @test data[22][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x63,0x68,0x61,0x72,0x6c,0x69,0x65]
        @test data[23][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x64,0x61,0x6e]
        @test string(data[24][1]) == "hey there ephraim"
        @test string(data[25][1]) == "hey there frank"
        @test string(data[26][1]) == "hey there george"
        @test string(data[27][1]) == "hey there hank"

        ODBC.execute!(dsn, "insert test1 VALUES
                            (2, -- bigint
                             1, -- bit
                             2.0, -- decimal
                             2, -- int
                             2.0, -- numeric
                             2, -- smallint
                             2, -- mediumint
                             2, -- tinyint
                             2.2, -- float
                             2.2, -- double
                             '2016-01-01', -- date
                             '2016-01-01 01:01:01', -- datetime
                             '2016-01-01 01:01:01', -- timestamp
                             '01:01:01', -- time
                             2016, -- year
                             'B', -- char(1)
                             'hey there sailor', -- varchar
                             cast('12' as binary(2)), -- binary
                             NULL, -- varbinary
                             'hey there abraham2', -- tinyblob
                             'hey there bill2', -- blob
                             'hey there charlie2', -- mediumblob
                             'hey there dan2', -- longblob
                             'hey there ephraim2', -- tinytext
                             'hey there frank2', -- text
                             'hey there george2', -- mediumtext
                             'hey there hank2' -- longtext
                            )")
        data = DataFrame(ODBC.Query(dsn, "select * from test1"))
        @test size(data) == (2,27)
        @test data[1][1] === Int64(1)
        @test data[1][2] === Int64(2)

        # @testset "Streaming mysql data to CSV" begin
        #     # Test exporting test1 to CSV
        #     temp_filename = "mysql_test1.csv"
        #     source = ODBC.Query(dsn, "select * from test1")
        #     CSV.write(temp_filename, source)

        #     open(temp_filename) do f
        #         @test readline(f) == (
        #             "test_bigint,test_bit,test_decimal,test_int,test_numeric," *
        #             "test_smallint,test_mediumint,test_tiny_int,test_float,test_real," *
        #             "test_date,test_datetime,test_timestamp,test_time,test_year," *
        #             "test_char,test_varchar,test_binary,test_varbinary,test_tinyblob," *
        #             "test_blob,test_mediumblob,test_longblob,test_tinytext,test_text," *
        #             "test_mediumtext,test_longtext"
        #         )
        #         @test readline(f) == (
        #             "1,1,1.0,1,1.0,1,1,1,1.2,1.2,2016-01-01,2016-01-01T01:01:01,2016-01-01T01:01:01,01:01:01,2016,A,hey there sailor,\"UInt8[0x31, 0x32]\",,\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x61, 0x62, 0x72, 0x61, 0x68, 0x61, 0x6d]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x62, 0x69, 0x6c, 0x6c]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x63, 0x68, 0x61, 0x72, 0x6c, 0x69, 0x65]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x64, 0x61, 0x6e]\",hey there ephraim,hey there frank,hey there george,hey there hank"
        #         )
        #         @test readline(f) == (
        #             "2,1,2.0,2,2.0,2,2,2,2.2,2.2,2016-01-01,2016-01-01T01:01:01,2016-01-01T01:01:01,01:01:01,2016,B,hey there sailor,\"UInt8[0x31, 0x32]\",,\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x61, 0x62, 0x72, 0x61, 0x68, 0x61, 0x6d, 0x32]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x62, 0x69, 0x6c, 0x6c, 0x32]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x63, 0x68, 0x61, 0x72, 0x6c, 0x69, 0x65, 0x32]\",\"UInt8[0x68, 0x65, 0x79, 0x20, 0x74, 0x68, 0x65, 0x72, 0x65, 0x20, 0x64, 0x61, 0x6e, 0x32]\",hey there ephraim2,hey there frank2,hey there george2,hey there hank2"
        #         )
        #     end
        #     rm(temp_filename)
        # end

        ODBC.execute!(dsn, "drop table if exists test1")
    end

    @testset "large query" begin
        ODBC.execute!(dsn, "drop table if exists test2")
        ODBC.execute!(dsn, """
        CREATE TABLE test2
        (
            ID INT NOT NULL PRIMARY KEY,
            first_name VARCHAR(25),
            last_name VARCHAR(25),
            Salary DECIMAL,
            `hourly rate` real,
            hireDate DATE,
            `last clockin` DATETIME
        );""")
        randoms = joinpath(dirname(pathof(ODBC)), "../test/data/randoms.csv")
        ODBC.execute!(dsn, "load data local infile '$randoms' into table test2
                            fields terminated by ',' lines terminated by '\n'
                            (id,first_name,last_name,salary,`hourly rate`,hiredate,`last clockin`)")

        data = DataFrame(ODBC.Query(dsn, "select count(*) from test2"))
        @test size(data) == (1,1)
        @test data[1][1] === 70000

        df = DataFrame(ODBC.Query(dsn, "select * from test2"))
        @test size(df) == (70000,7)
        @test df[1] == collect(1:70000)
        @test df[end][1] === ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0)
    end

    @testset "prepared statement" begin
        ODBC.execute!(dsn, "create table test3 as select * from test2 limit 0")
        ODBC.execute!(dsn, "delete from test3")

        stmt = ODBC.prepare(dsn, "insert into test3 values(?,?,?,?,?,?,?)")
        ODBC.execute!(stmt, [101, "Steve", "McQueen", 1.0, 100.0, Date(2016,1,1), DateTime(2016,1,1)])

        ODBC.execute!(stmt, [102, "Dean", "Martin", 1.5, 10.1, Date(2016,1,2), DateTime(2016,1,2)])

        df = DataFrame(ODBC.Query(dsn, "select * from test3"))
        @test size(df) == (2,7)
        @test df[1][end-1] == 101
        @test df[1][end] == 102
        @test df[2][end-1] == "Steve"
        @test df[2][end] == "Dean"
        @test df[3][end-1] == "McQueen"
        @test df[3][end] == "Martin"
        @test df[4][end-1] == DecFP.Dec64(1)
        @test df[4][end] == DecFP.Dec64(2)
        @test df[5][end-1] == 100.0
        @test df[5][end] == 10.1
        @test df[6][end-1] == ODBC.API.SQLDate(2016,1,1)
        @test df[6][end] == ODBC.API.SQLDate(2016,1,2)
        @test df[7][end-1] == ODBC.API.SQLTimestamp(2016,1,1,0,0,0,0)
        @test df[7][end] == ODBC.API.SQLTimestamp(2016,1,2,0,0,0,0)
        ODBC.execute!(dsn, "drop table if exists test2")
        ODBC.execute!(dsn, "drop table if exists test3")
    end

    ODBC.disconnect!(dsn)
end


text:
struct Memo
    items :: Vector{Any}
end
Memo() = Memo([])

function Base.push!(m::Memo, x)
    push!(m.items, x)
end

function Base.getindex(m::Memo, i::Integer)
    m.items[UInt(i)+UInt(1)]
end

function Base.setindex!(m::Memo, v, i::Integer)
    j = UInt(i) + UInt(1)
    length(m.items) < j && resize!(m.items, j)
    m.items[j] = v
    m
end


text:
using Surrogates
using Test

#1D
obj = x -> sin(x) + sin(x)^2 + sin(x)^3
lb = 0.0
ub = 10.0
x = sample(5,lb,ub,LowDiscrepancySample(2))
y = obj.(x)
p = 3.5
InverseDistance = InverseDistanceSurrogate(x,y,lb,ub,p=2.4)
InverseDistance_kwargs = InverseDistanceSurrogate(x,y,lb,ub)
prediction = InverseDistance(5.0)
add_point!(InverseDistance,5.0,-0.91)
add_point!(InverseDistance,[5.1,5.2],[1.0,2.0])

#ND

lb = [0.0,0.0]
ub = [10.0,10.0]
n = 100
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]^2
y = f.(x)
p = 3.0
InverseDistance = InverseDistanceSurrogate(x,y,p,lb,ub)
prediction = InverseDistance((1.0,2.0))
add_point!(InverseDistance,(5.0,3.4),-0.91)
add_point!(InverseDistance,[(5.1,5.2),(5.3,6.7)],[1.0,2.0])

# Multi-output #98
f  = x -> [x^2, x]
lb = 1.0
ub = 10.0
x  = sample(5, lb, ub, SobolSample())
push!(x, 2.0)
y  = f.(x)
surrogate = InverseDistanceSurrogate(x, y, lb, ub, p = 1.2)
surrogate_kwargs = InverseDistanceSurrogate(x, y, lb, ub)
@test surrogate(2.0) ≈ [4, 2]
surrogate((0.0, 0.0))

f  = x -> [x[1], x[2]^2]
lb = [1.0, 2.0]
ub = [10.0, 8.5]
x  = sample(20, lb, ub, SobolSample())
push!(x, (1.0, 2.0))
y  = f.(x)
surrogate = InverseDistanceSurrogate(x, y, lb, ub, p = 1.2)
surrogate_kwargs = InverseDistanceSurrogate(x,y,lb,ub)
@test surrogate((1.0, 2.0)) ≈ [1, 4]
x_new = (2.0, 2.0)
y_new = f(x_new)
add_point!(surrogate, x_new, y_new)
@test surrogate(x_new) ≈ y_new
surrogate((0.0, 0.0))


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

import REPL.LineEdit
import ReplMaker

ans = nothing

"""
	finished(s)

Examine the buffer in the repl to see if the input is complete
"""
function finished(s)
    str = String(take!(copy(LineEdit.buffer(s))))
    if length(str) == 0
        return false
    elseif occursin(r"(^|\n)[^%\n]*(;|\$)[ ]*(%[^\n]*)?$",str)
        return true
    else
        return false
    end
end

"""
    parserepl(input)

Called on valid input from the repl.
"""
function parserepl(input)
    if !isempty(strip(input))
        global ans = RExpr(input[1:end-1])
        if input[end] == ';'
            display(ans)
        end
    end
end

"""
    replinit(repl)

Start up the Reduce  repl mode.
"""
function repl_init(repl)
    ReplMaker.initrepl(
        parserepl,
        valid_input_checker=finished,
        prompt_text="reduce> ",
        prompt_color=:cyan,
        start_key='}',
        mode_name="reduce",
        repl=repl,
    )
end


text:
module CUSPARSE

using ..CuArrays
using ..CuArrays: libcusparse, unsafe_free!, @argout, @workspace

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

const SparseChar = Char

# core library
include("libcusparse_common.jl")
include("error.jl")
include("libcusparse.jl")

# low-level wrappers
include("array.jl")
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const handles_lock = ReentrantLock()
const created_handles = Dict{Tuple{UInt,Int},cusparseHandle_t}()
const active_handles = Vector{Union{Nothing,cusparseHandle_t}}()

function handle()
    tid = Threads.threadid()
    if @inbounds active_handles[tid] === nothing
        ctx = context()
        key = (objectid(ctx), tid)
        lock(handles_lock) do
            active_handles[tid] = get!(created_handles, key) do
                handle = cusparseCreate()
                atexit(()->CUDAdrv.isvalid(ctx) && cusparseDestroy(handle))
                handle
            end
        end
    end
    @inbounds active_handles[tid]
end

function __init__()
    resize!(active_handles, Threads.nthreads())
    fill!(active_handles, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        # we don't eagerly initialize handles, but do so lazily when requested
        active_handles[tid] = nothing
    end
end

end


text:
##################################################################################################
# Integration with UncertainData.jl (sampling from full supports of the furnishing distributions)
#################################################################################################

uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainValueDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# UncertainIndexDataset
UVX_idx = UncertainIndexDataset(uvals_x)
UVY_idx = UncertainIndexDataset(uvals_y)

# Real-valued vectors
x = resample.(uvals_x);
y = resample.(uvals_y);

# Vectors of uncertain values
uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainValueDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# UncertainIndexDataset
UVX_idx = UncertainIndexDataset(uvals_x)
UVY_idx = UncertainIndexDataset(uvals_y)

# Real-valued vectors
x = resample.(uvals_x);
y = resample.(uvals_y);

jtest = JointDistanceDistributionTest(B = 20)

@test causality(x, y, jtest) isa Vector{T} where T
@test causality(uvals_x, uvals_y, jtest) isa Vector{T} where T
@test causality(x, uvals_y, jtest) isa Vector{T} where T
@test causality(uvals_x, y, jtest) isa Vector{T} where T
@test causality(UVX, UVY, jtest) isa Vector{T} where T
@test causality(x, UVY, jtest) isa Vector{T} where T
@test causality(UVX, y, jtest) isa Vector{T} where T

# 20 intervals, so 20 values should be in the joint distance distribution
@test causality(x, y, jtest) |> length == 20
@test causality(uvals_x, uvals_y, jtest) |> length == 20
@test causality(x, uvals_y, jtest) |> length == 20
@test causality(uvals_x, y, jtest) |> length == 20
@test causality(UVX, UVY, jtest) |> length == 20
@test causality(x, UVY, jtest) |> length == 20
@test causality(UVX, y, jtest) |> length == 20

################################################################
# Integration with UncertainData.jl (with sampling constraints)
################################################################
onevar_constraints = ConstrainedResampling(TruncateStd(1))
twovar_constraints = ConstrainedResampling(TruncateStd(2), TruncateStd(1))

test_jdd = JointDistanceDistributionTest()
@test causality(x, y, test_jdd, twovar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, uvals_y, test_jdd, twovar_constraints) isa Vector{T} where T<:Real
@test causality(x, uvals_y, test_jdd, onevar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, y, test_jdd, onevar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, UVY, test_jdd, twovar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, UVY, test_jdd, twovar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, uvals_y, test_jdd, twovar_constraints) isa Vector{T} where T<:Real
@test causality(x, UVY, test_jdd, onevar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, y, test_jdd, onevar_constraints) isa Vector{T} where T<:Real


text:
function setadbackend(::Val{:forward_diff})
    Base.depwarn("`Turing.setadbackend(:forward_diff)` is deprecated. Please use `Turing.setadbackend(:forwarddiff)` to use `ForwardDiff`.", :setadbackend)
    setadbackend(Val(:forwarddiff))
end

function setadbackend(::Val{:reverse_diff})
    Base.depwarn("`Turing.setadbackend(:reverse_diff)` is deprecated. Please use `Turing.setadbackend(:tracker)` to use `Tracker` or `Turing.setadbackend(:reversediff)` to use `ReverseDiff`. To use `ReverseDiff`, please make sure it is loaded separately with `using ReverseDiff`.",  :setadbackend)
    setadbackend(Val(:tracker))
end


text:
# StdMap

eval(cxxnewfile("StdMap.cxx", """
    #include <map>
    #include <memory>
    #include <string>
    #include <utility>
    #include <vector>

    static_assert(sizeof(bool) == 1, "");
    """))

function generate(::Type{StdMap{K,T}}) where {K,T}
    CK = K == Bool ? "bool" : K == RefStdString ? "std::string" : K == SharedStdString ? "std::shared_ptr<std::string>" : cxxtype[K]
    NK = cxxname(CK)
    CT = T == Bool ? "bool" :
         T == RefStdString ? "std::string" :
         T == SharedStdString ? "std::shared_ptr<std::string>" :
         T == RefStdVector{RefStdString} ? "std::vector<std::string>" :
         T == RefStdVector{SharedStdString} ? "std::vector<std::shared_ptr<std::string>>" :
         T == SharedStdVector{RefStdString} ? "std::shared_ptr<std::vector<std::string>>" :
         T == SharedStdVector{SharedStdString} ? "std::shared_ptr<std::vector<std::shared_ptr<std::string>>>" : cxxtype[T]
    NT = cxxname(CT)

    eval(cxxfunction(FnName(:RefStdMap_new, "std_map_$(NK)_$(NT)_new", libSTL),
                     FnResult(Ptr{StdMap{K,T}}, "std::map<$CK, $CT> *", RefStdMap{K,T}, expr -> :(RefStdMap{$K,$T}($expr))),
                     [FnArg(:keytype, Nothing, "keytype", "void", Type{K}, identity; skip=true),
                      FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     "return new std::map<$CK, $CT>;"))

    eval(cxxfunction(FnName(:RefStdMap_delete, "std_map_$(NK)_$(NT)_delete", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "std::map<$CK, $CT> * restrict", RefStdMap{K,T}, identity)],
                     "delete map;"))

    eval(cxxfunction(FnName(:(Base.copy), "std_map_$(NK)_$(NT)_copy", libSTL),
                     FnResult(Ptr{StdMap{K,T}}, "std::map<$CK, $CT> *", RefStdMap{K,T}, expr -> :(RefStdMap{$K,$T}($expr))),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "const std::map<$CK, $CT> * restrict", RefStdMap{K,T}, identity)],
                     "return new std::map<$CK, $CT>(*map);"))

    # GCStdMap

    eval(cxxfunction(FnName(:GCStdMap_sizeof, "std_map_$(NK)_$(NT)_sizeof", libSTL), FnResult(Csize_t, "std::size_t"),
                     [FnArg(:keytype, Nothing, "keytype", "void", Type{K}, identity; skip=true),
                      FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     "return sizeof(std::map<$CK, $CT>);"))
    push!(selftests, () -> @assert GCStdMap_sizeof(K, T) <= GCStdMap_size)

    eval(cxxfunction(FnName(:GCStdMap_construct, "std_map_$(NK)_$(NT)_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdMap{K,T}}, "ptr", "void *", GCStdMap{K,T}, identity)], "new(ptr) std::map<$CK, $CT>;"))

    eval(cxxfunction(FnName(:GCStdMap_destruct, "std_map_$(NK)_$(NT)_destruct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdMap{K,T}}, "ptr", "std::map<$CK, $CT> * restrict", GCStdMap{K,T}, identity)],
                     "ptr->~map();"))

    eval(cxxfunction(FnName(:GCStdMap_copy_construct, "std_map_$(NK)_$(NT)_copy_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdMap{K,T}}, "ptr", "std::map<$CK, $CT> * restrict", GCStdMap{K,T}, identity),
                      FnArg(:vec, Ptr{StdMap{K,T}}, "vec", "const std::map<$CK, $CT> * restrict", GCStdMap{K,T}, identity)],
                     "new(ptr) std::map<$CK, $CT>(*vec);"))

    # SharedStdMap

    eval(cxxfunction(FnName(:SharedStdMap_sizeof, "std_shared_ptr_std_map_$(NK)_$(NT)_sizeof", libSTL),
                     FnResult(Csize_t, "std::size_t"),
                     [FnArg(:keytype, Nothing, "keytype", "void", Type{K}, identity; skip=true),
                      FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     "return sizeof(std::shared_ptr<std::map<$CK, $CT>>);"))
    push!(selftests, () -> @assert SharedStdMap_sizeof(K, T) <= SharedStdMap_size)

    eval(cxxfunction(FnName(:SharedStdMap_construct, "std_shared_ptr_std_map_$(NK)_$(NT)_placement_new", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdMap{K,T}, expr -> :(pointer_from_objref($expr)))],
                     """
                     auto res = new(ptr) std::shared_ptr<std::map<$CK, $CT>>;
                     *res = std::make_shared<std::map<$CK, $CT>>();
                     """))

    eval(cxxfunction(FnName(:SharedStdMap_destruct, "std_shared_ptr_std_map_$(NK)_$(NT)_placement_delete", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::map<$CK, $CT>> * restrict", SharedStdMap{K,T},
                            expr -> :(pointer_from_objref($expr)))], "ptr->~shared_ptr();"))

    eval(cxxfunction(FnName(:SharedStdMap_copy_construct, "std_shared_ptr_std_map_$(NK)_$(NT)_placement_copy", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdMap{K,T}, expr -> :(pointer_from_objref($expr))),
                      FnArg(:vec, Ptr{StdMap{K,T}}, "vec", "const std::map<$CK, $CT> * restrict", SharedStdMap{K,T}, identity)],
                     """
                     auto res = new(ptr) std::shared_ptr<std::map<$CK, $CT>>;
                     *res = std::make_shared<std::map<$CK, $CT>>(*vec);
                     """))

    eval(cxxfunction(FnName(:SharedStdMap_get, "std_shared_ptr_std_map_$(NK)_$(NT)_get", libSTL),
                     FnResult(Ptr{StdMap}, "std::map<$CK, $CT> *"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::map<$CK, $CT>> * restrict", SharedStdMap{K,T},
                            expr -> :(pointer_from_objref($expr)))], "return ptr->get();"))

    # StdMap

    eval(cxxfunction(FnName(:(Base.length), "std_map_$(NK)_$(NT)_length", libSTL),
                     FnResult(Csize_t, "std::size_t", Int, expr -> :(convert(Int, $expr))),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity)],
                     "return map->size();"))

    eval(cxxfunction(FnName(:(Base.haskey), "std_map_$(NK)_$(NT)_haskey", libSTL),
                     FnResult(Csize_t, "std::size_t", Bool, expr -> :(convert(Bool, $expr))),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity),
                      FnArg(:key, Ptr{K}, "key", "$CK const *", Any, expr -> :(convert_arg(Ptr{$K}, convert($K, $expr))))],
                     "return map->count(*key);"))

    eval(cxxfunction(FnName(:(Base.getindex), "std_map_$(NK)_$(NT)_getindex", libSTL),
                     FnResult(Ptr{T}, "$CT *", T, expr -> :(convert_result($T, $expr))),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "std::map<$CK, $CT> * restrict", StdMap{K,T}, identity),
                      FnArg(:key, Ptr{K}, "key", "$CK const *", Any, expr -> :(convert_arg(Ptr{$K}, convert($K, $expr))))],
                     "return &map->at(*key);"))

    eval(cxxfunction(FnName(:(Base.setindex!), "std_map_$(NK)_$(NT)_setindex_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "std::map<$CK, $CT> * restrict", StdMap{K,T}, identity),
                      FnArg(:elt, Ptr{T}, "elt", "$CT const *", Any, expr -> :(convert_arg(Ptr{$T}, convert($T, $expr)))),
                      FnArg(:key, Ptr{K}, "key", "$CK const *", Any, expr -> :(convert_arg(Ptr{$K}, convert($K, $expr))))],
                     "(*map)[*key] = *elt;"))

    eval(cxxfunction(FnName(:(Base.delete!), "std_map_$(NK)_$(NT)_delete_", libSTL),
                     FnResult(Nothing, "void", StdMap{K,T}, expr -> :map),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "std::map<$CK, $CT> * restrict", StdMap{K,T}, identity),
                      FnArg(:key, Ptr{K}, "key", "$CK const *", Any, expr -> :(convert_arg(Ptr{$K}, convert($K, $expr))))],
                     "map->erase(*key);"))

    eval(cxxfunction(FnName(:(Base.:(==)), "std_map_$(NK)_$(NT)_equals", libSTL), FnResult(Bool, "bool"),
                     [FnArg(:map1, Ptr{StdMap{K,T}}, "map1", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity),
                      FnArg(:map2, Ptr{StdMap{K,T}}, "map2", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity)],
                     "return *map1 == *map2;"))

    eval(cxxfunction(FnName(:cbegin, "std_map_$(NK)_$(NT)_cbegin", libSTL),
                     FnResult(StdMapIterator{K,T}, "std::map<$CK, $CT>::const_iterator"),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity)],
                     """
                     using MI = std::map<$CK, $CT>::const_iterator;
                     static_assert(sizeof(MI) == sizeof(void *), "");
                     return map->cbegin();
                     """))

    eval(cxxfunction(FnName(:cend, "std_map_$(NK)_$(NT)_cend", libSTL),
                     FnResult(StdMapIterator{K,T}, "std::map<$CK, $CT>::const_iterator"),
                     [FnArg(:map, Ptr{StdMap{K,T}}, "map", "const std::map<$CK, $CT> * restrict", StdMap{K,T}, identity)],
                     "return map->cend();"))

    # Iterators

    eval(cxxfunction(FnName(:StdMap_offsetof_pair_first, "std_map_$(NK)_$(NT)_offsetof_pair_first", libSTL),
                     FnResult(Csize_t, "std::size_t"),
                     [FnArg(:keytype, Nothing, "keytype", "void", Type{K}, identity; skip=true),
                      FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     """
                     using P = std::pair<$CK const *, $CT const *>;
                     return offsetof(P, first);
                     """))
    eval(cxxfunction(FnName(:StdMap_offsetof_pair_second, "std_map_$(NK)_$(NT)_offsetof_pair_second", libSTL),
                     FnResult(Csize_t, "std::size_t"),
                     [FnArg(:keytype, Nothing, "keytype", "void", Type{K}, identity; skip=true),
                      FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     """
                     using P = std::pair<$CK const *, $CT const *>;
                     return offsetof(P, second);
                     """))
    push!(selftests, () -> @assert StdMap_offsetof_pair_first(K, T) == fieldoffset(Pair{Ptr{K},Ptr{T}}, 1))
    push!(selftests, () -> @assert StdMap_offsetof_pair_second(K, T) == fieldoffset(Pair{Ptr{K},Ptr{T}}, 2))
    eval(cxxfunction(FnName(:(Base.getindex), "std_map_$(NK)_$(NT)_const_iterator_getindex", libSTL),
                     FnResult(Pair{Ptr{K},Ptr{T}}, "const std::pair<$CK const *, $CT const *>", Pair{K,T},
                              expr -> :(convert_result($K, $expr[1]) => convert_result($T, $expr[2]))),
                     [FnArg(:iter, StdMapIterator{K,T}, "iter", "std::map<$CK, $CT>::const_iterator")],
                     "return std::pair<$CK const *, $CT const *>(&iter->first, &iter->second);"))

    eval(cxxfunction(FnName(:inc, "std_map_$(NK)_$(NT)_const_iterator_inc", libSTL),
                     FnResult(StdMapIterator{K,T}, "std::map<$CK, $CT>::const_iterator"),
                     [FnArg(:iter, StdMapIterator{K,T}, "iter", "std::map<$CK, $CT>::const_iterator")], "return ++iter;"))

    eval(cxxfunction(FnName(:dec, "std_map_$(NK)_$(NT)_const_iterator_dec", libSTL),
                     FnResult(StdMapIterator{K,T}, "std::map<$CK, $CT>::const_iterator"),
                     [FnArg(:iter, StdMapIterator{K,T}, "iter", "std::map<$CK, $CT>::const_iterator")], "return --iter;"))

    eval(cxxfunction(FnName(:(Base.:(==)), "std_map_$(NK)_$(NT)_const_iterator_equals", libSTL), FnResult(Bool, "bool"),
                     [FnArg(:iter1, StdMapIterator{K,T}, "iter1", "std::map<$CK, $CT>::const_iterator"),
                      FnArg(:iter2, StdMapIterator{K,T}, "iter2", "std::map<$CK, $CT>::const_iterator")], "return iter1 == iter2;"))

    return nothing
end

const StdMap_types = filter(T -> !(T <: Complex), value_types) ∪
                     Set([RefStdString, SharedStdString, RefStdVector{RefStdString}, RefStdVector{SharedStdString},
                          SharedStdVector{RefStdString}, SharedStdVector{SharedStdString}])
const StdMap_keys = filter(T -> T <: Union{Integer,StdString}, StdMap_types)
for K in sort!(collect(StdMap_keys); by=string), T in sort!(collect(StdMap_types); by=string)
    generate(StdMap{K,T})
end

Base.copy(vec::GCStdMap{K,T}) where {K,T} = GCStdMap{K,T}(Base.copy, vec)
Base.copy(vec::SharedStdMap{K,T}) where {K,T} = SharedStdMap{K,T}(Base.copy, vec)

free(map::RefStdMap) = RefStdMap_delete(map)

################################################################################

Base.convert(::Type{Dict{K,T}}, map::StdMap{K,T}) where {K,T} = Dict{K,T}(k => v for (k, v) in map)
Base.convert(::Type{Dict}, map::StdMap{K,T}) where {K,T} = convert(Dict{K,T}, map)

Base.isempty(map::StdMap) = length(map) == 0

Base.keytype(::StdMap{K}) where {K} = K
Base.eltype(::StdMap{K,T}) where {K,T} = T

function Base.iterate(map::StdMap)
    iter = cbegin(map)
    iter == cend(map) && return nothing
    return iter[], iter
end
function Base.iterate(map::StdMap, iter::StdMapIterator)
    iter = inc(iter)
    iter == cend(map) && return nothing
    return iter[], iter
end

################################################################################

# GCStdMap{K,T}() where {K,T} = GCStdMap{K,T}(StdMap{K,T}())
#
# GCStdMapIterator{K,T}() where {K,T} = GCStdMapIterator{K,T}(StdMapIterator{K,T}())
#
# free(map::GCStdMap) = free(map.managed)
# free(iter::GCStdMapIterator) = free(iter.managed)
#
# Base.length(map::GCStdMap) = length(map.managed)
# Base.isempty(map::GCStdMap) = isempty(map.managed)
# Base.haskey(map::GCStdMap, key) = haskey(map.managed, key)
# Base.getindex(map::GCStdMap, key) = getindex(map.managed, key)
# Base.setindex!(map::GCStdMap, elt, key) = setindex!(map.managed, elt, key)
# Base.delete!(map::GCStdMap, key) = delete!(map.managed, key)
# Base.size(map::GCStdMap) = size(map.managed)
# Base.keytype(::GCStdMap{K}) where {K} = keytype(StdMap{K})
# Base.eltype(::GCStdMap{K,T}) where {K,T} = eltype(StdMap{K,T})

################################################################################

# cbegin(map::GCStdMap{K,T}) where {K,T} = GCStdMapIterator{K,T}(cbegin(map.managed))
# cend(map::GCStdMap{K,T}) where {K,T} = GCStdMapIterator{K,T}(cend(map.managed))
# Base.:(==)(iter1::GCStdMapIterator{K,T}, iter2::GCStdMapIterator{K,T}) where {K,T} = iter1.managed == iter2.managed
# Base.getindex(iter::GCStdMapIterator) = getindex(iter.managed)
# inc!(iter::GCStdMapIterator) = inc!(iter.managed)
# dec!(iter::GCStdMapIterator) = dec!(iter.managed)
# is_cbegin(iter::GCStdMapIterator{K,T}, map::GCStdMap{K,T}) where {K,T} = is_cbegin(iter.managed, map.managed)
# is_end(iter::GCStdMapIterator{K,T}, map::GCStdMap{K,T}) where {K,T} = is_end(iter.managed, map.managed)
#
# function Base.iterate(map::StdMap{K,T}) where {K,T}
#     iter = GCStdMapIterator{K,T}(cbegin(map))
#     is_cend(iter.managed, map) && return nothing
#     return (iter[], iter)
# end
# function Base.iterate(map::StdMap{K,T}, iter::GCStdMapIterator{K,T}) where {K,T}
#     inc!(iter)
#     is_cend(iter.managed, map) && return nothing
#     return (iter[], iter)
# end
#
# function Base.iterate(map::GCStdMap{K,T}) where {K,T}
#     iter = cbegin(map)
#     is_cend(iter, map) && return nothing
#     return (iter[], iter)
# end
# function Base.iterate(map::GCStdMap{K,T}, iter::GCStdMapIterator{K,T}) where {K,T}
#     inc!(iter)
#     is_cend(iter, map) && return nothing
#     return (iter[], iter)
# end


text:
# PROJECT_BOX    Project a point onto a box (hyper-rectangle).
#
#   project_box(v,l,u) is the projection of v onto
#   the set { x | l <= x <= u }.

function project_box(
	v,
	l,
	u
)
    x = max(l, min(v, u))
    return x
end


text:
using ShirleyRayTracer
using Test

@testset "ShirleyRayTracer.jl" begin
    # Write your tests here.
end


text:
@testset "sdfilewriter" begin

@testset "molblock" begin
    demomol = joinpath(dirname(@__FILE__), "..", "assets", "test", "demo.mol")
    mol = sdftomol(demomol)
    buf = IOBuffer(write=true)
    sdfilewriter(buf, (mol,))
    mol2 = sdftomol(split(String(take!(buf)), "\n"))
    @test nodecount(mol2) == 37
end

end #sdfilewriter


text:
module WRF

#
# Module for reading in WRF outputs and re-packaging selected variables
# for use with the moisture simulation/data assimilation code.
#
#

using Calendar
import Calendar.CalendarTime

using NetCDF

type WRFData

    # source file path
    file_path :: String

    # the variable fields
    fields :: Dict{String,Array{Float64}}

    # the times of recording history
    tm :: Array{CalendarTime}
end

WRFData(file_path) = WRFData(file_path, Dict{String,Array{Float64}}(), Array(CalendarTime,0))

function load_wrf_data(file_path::String, var_list::Array{ASCIIString})

    w = WRFData(file_path)
    nc = NetCDF.open(file_path)

    vl = Array(ASCIIString, 0)
    def_var_list = ["T2", "Q2", "PSFC", "RAINC", "RAINNC"]
    if var_list != nothing
       append!(vl, var_list)
    end
    append!(vl, def_var_list)

    for v in vl
        fld = NetCDF.readvar(nc, v, [1, 1, 1], [-1, -1, -1])
        w.fields[v] = fld
    end

    tm = NetCDF.readvar(nc, "Times", [1, 1], [-1, -1])'
    for i in 1:size(tm,1)
        push!(w.tm, Calendar.parse("yyyy-MM-dd_HH:mm:ss", ascii(squeeze(tm[i,:], 1)), "GMT"))
    end

    # read in grid dimensions
    w.fields["lon"] = squeeze(NetCDF.readvar(nc, "XLONG", [1, 1, 1], [-1, -1, -1])[:,:,1], 3)
    w.fields["lat"] = squeeze(NetCDF.readvar(nc, "XLAT", [1, 1, 1], [-1, -1, -1])[:,:,1], 3)

    # compute derived variables
    compute_rainfall_per_hour(w)
    compute_equilibrium_moisture(w)

    NetCDF.close(nc)

    return w
end

function lat(w::WRFData)
    return w.fields["lat"]
end

function lon(w::WRFData)
    return w.fields["lon"]
end

function times(w::WRFData)
    return w.tm
end

function field(w::WRFData, fname)
    return w.fields[fname]
end

function interpolated_field(w::WRFData, fname)
    F = field(w, fname)
    IF = zeros(Float64, size(F))
    N = size(F,3)
    IF[:,:,2:N] = 0.5 * (F[:,:,1:N-1] + F[:,:,2:N])
    return IF
end

function compute_rainfall_per_hour(w::WRFData)
    rainnc = field(w, "RAINNC")
    rainc = field(w, "RAINC")
    rain = zeros(Float64, size(rainnc))
    rain_old = zeros(Float64, size(rain[:,:,1]))

    for i in 2:length(w.tm)
        dt = (w.tm[i] - w.tm[i-1]).millis / 1000
        rain[:,:,i] = ((rainc[:,:,i] + rainnc[:,:,i]) - rain_old) * 3600.0 / dt
        rain_old[:,:,:] = rainc[:,:,i]
        rain_old += rainnc[:,:,i]
     end

     w.fields["RAIN"] = rain
     delete!(w.fields, "RAINC")
     delete!(w.fields, "RAINNC")
end

function domain_extent(w::WRFData)
    lon = field(w, "LON")
    lat = field(w, "LAT")

    return (min(lon), min(lat), max(lon), max(lat))
end

function compute_equilibrium_moisture(w::WRFData)
    # pull out required fields from WRF
    P, Q, T = field(w, "PSFC"), field(w, "Q2"), field(w, "T2")

    N = size(P,3)

    # Xi stands for X interpolated
    Pi = 0.5 * (P[:,:,1:N-1] + P[:,:,2:N])
    Qi = 0.5 * (Q[:,:,1:N-1] + Q[:,:,2:N])
    Ti = 0.5 * (T[:,:,1:N-1] + T[:,:,2:N])

    # compute saturated vapor pressure
    Pws = exp(54.842763 - 6763.22/Ti - 4.210 * log(Ti) + 0.000367*Ti
              + tanh(0.0415*(Ti - 218.8)) .*
              (53.878 - 1331.22/Ti - 9.44523 * log(Ti) + 0.014025*Ti))

    # compute current water vapor pressure
    Pw = Pi .* Qi ./ (0.622 + (1 - 0.622) * Qi)

    # compute relative humidity in percent
    RH = 100 * Pw ./ Pws

    # drying/wetting fuel equilibrium moisture contents
    Ed = zeros(Float64, size(P))
    Ew = zeros(Float64, size(P))

    Ed[:,:,2:] = 0.924*RH.^0.679 + 0.000499*exp(0.1*RH) + 0.18*(21.1 + 273.15 - Ti).*(1 - exp(-0.115*RH))
    Ew[:,:,2:] = 0.618*RH.^0.753 + 0.000454*exp(0.1*RH) + 0.18*(21.1 + 273.15 - Ti).*(1 - exp(-0.115*RH))

    Ed *= 0.01
    Ew *= 0.01

    w.fields["Ed"] = Ed
    w.fields["Ew"] = Ew

end

function slice_field(w::WRFData, fn::String)
    # remove the time dimension from a static field and save a lot of storage
    w.fields[fn] = squeeze(w.fields[fn][:,:,1], 3)
end

end


text:
@testset "Scan" begin
    function stub_intensity_scan()
        scans = ["scan1d.jld2", "scan2d.jld2", "scan3d.jld2"]
        scans = joinpath.(@__DIR__, scans)
        return load(scans[1]), load(scans[2]), load(scans[3])
    end

    params = ScanParameters(
        medium = Medium(),
        excitation = Excitation(),
        f0 = 15e6,
        hydrophone_id = :Onda_HGL0200_2322,
        preamp_id = :Onda_AH2020_1238_20dB,
    )

    scan1d, scan2d, scan3d = stub_intensity_scan();

    data_1d = IntensityScans.ScanData(scan1d.version, scan1d.axes, scan1d.scope_info, scan1d.time, scan1d.waveform, scan1d.coordinates, scan1d.metrics)
    data_2d = IntensityScans.ScanData(scan2d.version, scan2d.axes, scan2d.scope_info, scan2d.time, scan2d.waveform, scan2d.coordinates, scan2d.metrics)
    data_3d = IntensityScans.ScanData(scan3d.version, scan3d.axes, scan3d.scope_info, scan3d.time, scan3d.waveform, scan3d.coordinates, scan3d.metrics)

    metrics1d = compute_metrics(data_1d, params)
    metrics2d = compute_metrics(data_2d, params)
    metrics3d = compute_metrics(data_3d, params)
end


text:
@enum Color::UInt Black Red
msb(::Type{T}) where {T<:Integer} = bitrotate(one(T), -1)
msb(::T) where {T<:Integer} = msb(T)

"""
    ProtoNode{T}

Immutable struct for representing a tree node.
The main idea here is that each `ProtoNode` is backed be a vector of such nodes held in
the Red-Black tree itself, and as such has a (relatively) stable address (modulo resizing
of the underlying array, which will have to be dealt with).

Mutation of the `ProtoNode` must happen in conjunction with loads and stores to its
underlying array.

**NOTE**: We store the color of the `ProtoNode` implicitly in the most significant bit of
the `parent` field. As such, care must be taken when accessing that field.
"""
struct ProtoNode{T}
    parent::UInt
    left::UInt
    right::UInt
    nchildren::UInt
    key::T
end

"""
    RedBlackTree{T}

A Red-Black tree holding elements of type `T`.
The tree itself is distributed throughout the `nodes` field with the `root` field pointing
to the root of the tree (or set to zero if the tree is empty).

Unused indices are stored in the `free` field and can be used to assign new `ProtoNodes`
to items inserted into the tree.
"""
mutable struct RedBlackTree{T}
    root::UInt
    nodes::Vector{ProtoNode{T}}
    free::Vector{Int}
end

Base.length(tree::RedBlackTree) = length(tree.nodes) - length(tree.free)

function RedBlackTree{T}(len::Integer = 10) where {T}
    # Due to our load-store schenanigans, make sure that we only hold `isbitstypes`.
    isbitstype(T) || error("Can only construct a RedBlackTree for `isbitstype` key types!")

    root = zero(Int)
    nodes = Vector{ProtoNode{T}}(undef, max(one(len), len))
    free = collect(reverse(eachindex(nodes)))
    return RedBlackTree{T}(root, nodes, free)
end

"""
    Base.getindex(tree::RedBlackTree) -> Int

Return a free `index` in `tree` that points to a free `ProtoNode`.

Implementation Details
----------------------
If no such index exists, then the underlying `nodes` vector will be doubled in size to
create more nodes.
"""
function Base.getindex(tree::RedBlackTree)
    nodes, free = tree.nodes, tree.free
    if isempty(free)
        currentlength = length(nodes)
        resize!(nodes, 2 * currentlength)
        for i = (currentlength + 1):length(nodes)
            push!(free, i)
        end
    end
    return pop!(free)
end
Base.getindex(tree::RedBlackTree, index::Integer) = tree.nodes[index].key

"""
    return!(tree::RedBlackTree, index)

Return `index`'s `ProtoNode` back to `tree` for future reuse.
"""
return!(tree::RedBlackTree, index::Integer) = push!(tree.free, index)

function Base.setindex!(tree::RedBlackTree{T}, k0, index::Integer) where {T}
    k = convert(T, k0)
    ptr = Ptr{T}(pointer(tree, index)) + 4 * sizeof(UInt)
    return unsafe_store!(ptr, k)
end

@inline Base.pointer(tree::RedBlackTree, i::Integer) = pointer(tree.nodes, i)

"""
    getparent(tree::RedBlackTree, n) -> Integer

Return the parent index for the the node at `n` in `tree`.
"""
function getparent(tree::RedBlackTree{T}, index::Integer) where {T}
    ptr = pointer(tree, index)
    unmasked = unsafe_load(Ptr{UInt}(ptr))
    return unmasked & ~msb(unmasked)
end

"""
    setparent!(tree::RedBlackTree, n, p)

Assign `p` as `n`'s parent.
"""
function setparent!(tree::RedBlackTree, n::Integer, p::Integer)
    ptr = Ptr{UInt}(pointer(tree, n))
    old = unsafe_load(ptr)
    new = (old & msb(old)) | (p & ~msb(old))
    return unsafe_store!(ptr, new)
end

function getcolor(tree::RedBlackTree{T}, index::Integer) where {T}
    ptr = pointer(tree, index)
    unmasked = unsafe_load(Ptr{UInt}(ptr))
    # N.B.: LLVM is able to elide the bounds check for the enum creation, so we don't
    # need to get fancy with bitcasts or anything like that.
    return Color(bitrotate(unmasked & msb(unmasked), 1))
end

function setcolor!(tree::RedBlackTree, index::Integer, color)
    ptr = Ptr{UInt}(pointer(tree, index))
    old = unsafe_load(ptr)
    new = (old & ~msb(old)) | bitrotate(UInt(color) & one(UInt), -1)
    return unsafe_store!(ptr, new)
end

function getchild(tree::RedBlackTree, index::Integer, direction::Direction)
    ptr = pointer(tree, index) + (1 + UInt(direction)) * sizeof(UInt)
    return unsafe_load(Ptr{UInt}(ptr))
end

leftchild(tree::RedBlackTree, index) = getchild(tree, index, Left)
rightchild(tree::RedBlackTree, index) = getchild(tree, index, Right)

function setchild!(tree::RedBlackTree, index::Integer, direction::Direction, child)
    ptr = pointer(tree, index) + (1 + UInt(direction)) * sizeof(UInt)
    return unsafe_store!(Ptr{UInt}(ptr), convert(UInt, child))
end
setleft!(tree::RedBlackTree, index, child) = setchild!(tree, index, Left, child)
setright!(tree::RedBlackTree, index, child) = setchild!(tree, index, Right, child)

isred(tree::RedBlackTree, n::Integer) = !iszero(n) && (getcolor(tree, n) == Red)

# Sub-tree size tracking
function subtreesize(tree::RedBlackTree, n::Integer)
    ptr = pointer(tree, n) + 3 * sizeof(UInt)
    return unsafe_load(Ptr{UInt}(ptr))
end

function setsubtree!(tree::RedBlackTree, n::Integer, v)
    ptr = pointer(tree, n) + 3 * sizeof(UInt)
    unsafe_store!(Ptr{UInt}(ptr), convert(UInt, v))
end

function setsubtree!(tree::RedBlackTree, n::Integer)
    nchildren = 1
    left = leftchild(tree, n)
    iszero(left) || (nchildren += subtreesize(tree, left))

    right = rightchild(tree, n)
    iszero(right) || (nchildren += subtreesize(tree, right))
    setsubtree!(tree, n, nchildren)
    return nothing
end

#####
##### Convenience Accessors
#####

function child_direction(tree::RedBlackTree, p::Integer, n::Integer)
    return Direction(leftchild(tree, p) != n)
end

child_direction(tree::RedBlackTree, n::Integer) =
    child_direction(tree, getparent(tree, n), n)

function rotateroot!(tree::RedBlackTree, p::Integer, dir::Direction)
    g = getparent(tree, p)
    s = getchild(tree, p, reverse(dir))
    c = getchild(tree, s, dir)
    setchild!(tree, p, reverse(dir), c)
    iszero(c) || setparent!(tree, c, p)

    setchild!(tree, s, dir, p)
    setparent!(tree, p, s)
    setparent!(tree, s, g)
    iszero(g) ? (tree.root = s) : setchild!(tree, g, child_direction(tree, g, p), s)

    # update subtree sizes
    setsubtree!(tree, p )
    setsubtree!(tree, s)
    return s
end

rotateleft!(tree::RedBlackTree, n::Integer) = rotateroot!(tree, n, Left)
rotateright!(tree::RedBlackTree, n::Integer) = rotateroot!(tree, n, Right)

# Now that we have the nitty-gritty accessor implementation out of the way, we can
# get to the actual tree implementation.
@enum Findings::UInt8 NotFoundLeft NotFoundRight Found
function search(tree::RedBlackTree{T}, key::T) where {T}
    n = tree.root
    lastn = n
    direction = Left

    while !iszero(n)
        v = tree[n]
        key == v && return n, Found
        direction = Direction(v < key)
        lastn, n = n, getchild(tree, n, direction)
    end
    return lastn, Base.bitcast(Findings, direction)
end

function Base.haskey(tree::RedBlackTree{T}, k0) where {T}
    _, status = search(tree, convert(T, k0))
    return status == Found
end
Base.in(k, tree::RedBlackTree) = haskey(tree, k)

#####
##### Insertion
#####

function Base.push!(tree::RedBlackTree{T}, k::T) where {T}
    p, status = search(tree, k)
    status == Found && return tree

    # Add 1 to all subtree sizes to the nodes seen along the way.
    g = p
    while !iszero(g)
        setsubtree!(tree, g, subtreesize(tree, g) + 1)
        g = getparent(tree, g)
    end

    # Create a new node and finish inserting
    n = tree[]
    tree[n] = k
    _insert!(tree, n, p, Base.bitcast(Direction, status))
    return tree
end

"""
    _insert!(tree::RedBlackTree, n, p, direction)

Insert node `n` into the tree as the `direction` child under `p`.
"""
function _insert!(tree::RedBlackTree, n::Integer, p::Integer, direction::Direction)
    setcolor!(tree, n, Red)
    setparent!(tree, n, p)
    setleft!(tree, n, 0)
    setright!(tree, n, 0)
    setsubtree!(tree, n, 1)

    # There is no parent.
    # `n` is the new root of the tree and insertion is complete.
    if iszero(p)
        tree.root = n
        return nothing
    end

    # Set `n` as the `direction` child of `p`.
    setchild!(tree, p, direction, n)

    # Insertion procedure
    local g
    while true
        # Case 1: Parent is Black, insertion is complete.
        if getcolor(tree, p) == Black
            return nothing
        end

        # From now on, we know `getcolor(tree, p) == Red`.
        g = getparent(tree, p)
        if iszero(g)
            # Case 4: `p` is the root and Red.
            setcolor!(tree, p, Black)
            return nothing
        end

        # `p` is `Red` and `g` exists.
        direction = child_direction(tree, g, p)

        # `u` is `n`'s uncle.
        u = getchild(tree, g, reverse(direction))
        if iszero(u) || getcolor(tree, u) == Black
            break
        end

        # Case 2: (Parent and Uncle are Red)
        setcolor!(tree, p, Black)
        setcolor!(tree, u, Black)
        setcolor!(tree, g, Red)

        # `g` becomes the new current node.
        # Iterate one black level higher.
        n = g
        p = getparent(tree, n)

        # Case 3 - `p` is NULL and insertion is complete.
        iszero(p) && return nothing
    end

    # Case 5: `p` is Red and `u` is black
    if n == getchild(tree, p, reverse(direction))
        rotateroot!(tree, p, direction)
        n = p
        p = getchild(tree, g, direction)
    end

    # Case 6: `p` is Red, `u` is black, and `n` is the outer grandchild of `g`.
    rotateroot!(tree, g, reverse(direction))
    setcolor!(tree, p, Black)
    setcolor!(tree, g, Red)
    return nothing
end

#####
##### Deletion
#####

function swap_successor!(tree::RedBlackTree, n::Integer, m::Integer)
    # Swap colors
    cₙ = getcolor(tree, n)
    cₘ = getcolor(tree, m)
    setcolor!(tree, m, cₙ)
    setcolor!(tree, n, cₘ)

    # Swap Left Children
    lₙ = leftchild(tree, n)
    lₘ = leftchild(tree, m)
    setleft!(tree, n, lₘ)
    iszero(lₘ) || setparent!(tree, lₘ, n)
    setleft!(tree, m, lₙ)
    iszero(lₙ) || setparent!(tree, lₙ, m)

    # Need to be more careful with right children and parents for the case where
    # `m` is the right child of `n`
    if n == getparent(tree, m)
        ##println("    Successor is child")
        p = getparent(tree, n)
        setparent!(tree, m, p)
        iszero(p) || setchild!(tree, p, child_direction(tree, p, n), m)

        rₘ = rightchild(tree, m)
        setright!(tree, n, rₘ)
        iszero(rₘ) || setparent!(tree, rₘ, n)
        setparent!(tree, n, m)
        setright!(tree, m, n)
    else
        # Swap parents
        #println("    Successor is not child")
        pₙ = getparent(tree, n)
        pₘ = getparent(tree, m)
        setparent!(tree, m, pₙ)
        iszero(pₙ) || setchild!(tree, pₙ, child_direction(tree, pₙ, n), m)
        setparent!(tree, n, pₘ)
        setchild!(tree, pₘ, child_direction(tree, pₘ, m), n)

        # Swap right children
        rₙ = rightchild(tree, n)
        rₘ = rightchild(tree, m)
        setright!(tree, n, rₘ)
        iszero(rₘ) || setparent!(tree, rₘ, n)
        setright!(tree, m, rₙ)
        setparent!(tree, rₙ, m)
    end
    n == tree.root && (tree.root = m)
    setsubtree!(tree, n)
    setsubtree!(tree, m)
    return nothing
end

function Base.findmin(tree::RedBlackTree, n::Integer)
    while true
        m = leftchild(tree, n)
        iszero(m) && return n
        n = m
    end
end

function _delete!(tree::RedBlackTree, n::Integer)
    # Case 1: `n` is the only node in the tree
    if n == tree.root && length(tree) == 1
        tree.root = 0
        return nothing
    end

    # Case 2: Do we have two non-NIL children?
    left = leftchild(tree, n)
    right = rightchild(tree, n)
    if !iszero(left) && !iszero(right)
        replacement = findmin(tree, right)
        swap_successor!(tree, n, replacement)

        # Update `left` and `right`, then fall through to case 3.
        left = leftchild(tree, n)
        right = rightchild(tree, n)
    end

    # Case 3: `n` now has at most one non-NIL child.
    color = getcolor(tree, n)
    # If `n` is Red - then it has no children.
    # Else if `n` is black, is it has a child, it must be a red child.
    # Replace `node` with its child and color the child `Black`.
    p = getparent(tree, n)

    # Decrement subtree size from `n` to the root.
    g = p
    while !iszero(g)
        setsubtree!(tree, g, subtreesize(tree, g) - 1)
        g = getparent(tree, g)
    end

    if color == Red
        setchild!(tree, p, child_direction(tree, p, n), 0)
        return nothing
    elseif iszero(left) && iszero(right)
        _delete_complex!(tree, n)
        return nothing
    else
        m = iszero(left) ? right : left
        setcolor!(tree, m, Black)
        setparent!(tree, m, p)
        if iszero(p)
            tree.root = m
        else
            setchild!(tree, p, child_direction(tree, p, n), m)
        end
        return nothing
    end
end

function _delete_complex!(tree::RedBlackTree, n::Integer)
    p = getparent(tree, n)
    direction = child_direction(tree, p, n)
    setchild!(tree, p, direction, 0)

    local s, d, c
    while true
        # s: sibling
        # d: distant nephew
        # c: close nephew
        s = getchild(tree, p, reverse(direction))
        d = getchild(tree, s, reverse(direction))
        c = getchild(tree, s, direction)

        getcolor(tree, s) == Red && @goto case_d3
        isred(tree, d) && @goto case_d6
        isred(tree, c) && @goto case_d5
        getcolor(tree, p) == Red && @goto case_d4

        # All of p, c, s, d are Black
        setcolor!(tree, s, Red)
        n = p
        p = getparent(tree, n)
        iszero(p) && return nothing
        direction = child_direction(tree, p, n)
    end

    @label case_d3
    rotateroot!(tree, p, direction)
    setcolor!(tree, p, Red)
    setcolor!(tree, s, Black)
    s = c
    d = getchild(tree, s, reverse(direction))
    isred(tree, d) && @goto case_d6
    c = getchild(tree, s, direction)
    isred(tree, c) && @goto case_d5

    @label case_d4
    setcolor!(tree, s, Red)
    setcolor!(tree, p, Black)
    return nothing

    @label case_d5
    rotateroot!(tree, s, reverse(direction))
    setcolor!(tree, s, Red)
    setcolor!(tree, c, Black)
    d = s
    s = c

    @label case_d6
    rotateroot!(tree, p, direction)
    setcolor!(tree, s, getcolor(tree, p))
    setcolor!(tree, p, Black)
    setcolor!(tree, d, Black)
    return nothing
end

function Base.delete!(tree::RedBlackTree{T}, k0) where {T}
    k = convert(T, k0)
    n, status = search(tree, k)
    if status == Found
        _delete!(tree, n)
        push!(tree.free, n)
        return true
    end
    return false
end

#####
##### Count Children
#####

function nchildren(tree::RedBlackTree{T}, k0) where {T}
    k = convert(T, k0)
    count = 0
    n = tree.root
    while !iszero(n)
        v = tree[n]
        if v == k
            right = rightchild(tree, n)
            iszero(right) || (count += subtreesize(tree, right))
            return count
        elseif k < v
            count += 1
            right = rightchild(tree, n)
            iszero(right) || (count += subtreesize(tree, right))
            n = leftchild(tree, n)
        else
            n = rightchild(tree, n)
        end
    end
end

#####
##### misc
#####

function traverse(f::F, tree::RedBlackTree) where {F}
    root = tree.root
    iszero(root) || traverse(f, tree, tree.root)
    return nothing
end

function traverse(f::F, tree::RedBlackTree, i) where {F}
    left = leftchild(tree, i)
    iszero(left) || traverse(f, tree, left)
    f(i)
    right = rightchild(tree, i)
    iszero(right) || traverse(f, tree, right)
    return nothing
end

function istree(tree::RedBlackTree{T}) where {T}
    v = T[]
    seen = Set{Int}()
    traverse(tree) do i
        if in(i, seen)
            error("Tree is broken!")
        end
        push!(seen, i)
        push!(v, tree[i])
    end
    return issorted(v)
end

function no_red_children(tree::RedBlackTree)
    passed = true
    traverse(tree) do i
        if getcolor(tree, i) == Red
            left = leftchild(tree, i)
            right = rightchild(tree, i)
            passed &= (iszero(left) || getcolor(tree, left) == Black)
            passed &= (iszero(right) || getcolor(tree, right) == Black)
        end
    end
    return passed
end

"""
    pathsequal(tree::RedBlackTree)

Check that all paths to all NIL decendents for each node in `tree` passes through the same
number of black nodes.
"""
function pathsequal(tree::RedBlackTree)
    root = tree.root
    iszero(root) && return true
    return pathlength(tree, tree.root) != 0
end

function pathlength(tree::RedBlackTree, i)
    left = leftchild(tree, i)
    right = rightchild(tree, i)
    if iszero(left) && iszero(right)
        return Int(getcolor(tree, i) == Black) + 1
    else
        left_length = iszero(left) ? 1 : pathlength(tree, left)
        right_length = iszero(right) ? 1 : pathlength(tree, right)

        # We've failed - abort!
        if iszero(left_length) || iszero(right_length) || left_length != right_length
            println("Node: $i - Left: $left_length, Right: $right_length")
            return 0
        else
            return left_length + Int(getcolor(tree, i) == Black)
        end
    end
end

function validate(tree::RedBlackTree)
    return istree(tree) && no_red_children(tree) && pathsequal(tree)
end

#####
##### Abstract Trees
#####

struct TreeIndex{T,I}
    tree::RedBlackTree{T}
    index::I
end

AbstractTrees.print_tree(tree::RedBlackTree; kw...) =
    AbstractTrees.print_tree(TreeIndex(tree, tree.root); kw...)

function AbstractTrees.children(ti::TreeIndex)
    left, right = leftchild(ti.tree, ti.index), rightchild(ti.tree, ti.index)
    return [TreeIndex(ti.tree, i) for i in (left, right) if !iszero(i)]
end

function AbstractTrees.printnode(io::IO, ti::TreeIndex)
    metadata = join((ti.index, getcolor(ti.tree, ti.index), subtreesize(ti.tree, ti.index)), " - ")
    key = ti.tree[ti.index]
    print(io, "$key ($metadata)")
end