text:
function plot_parametric_surface(func; n_points=101, axis_title="", axis=nothing, camera=(25, -135), colorfunc=nothing, rel_colorfunc=nothing, colormap="viridis", kwargs...)
    # Plot Data
    xs = ys = range(0, 1, length=n_points)

    F = [func([xi, yi]) for xi in xs, yi in xs]
    fx = [el[1] for el in F]
    fy = [el[2] for el in F]
    fz = [el[3] for el in F]

    if colorfunc !== nothing
        colmap = [colorfunc([xi, yi]) for xi in xs, yi in xs]
        # Check relative colofunc
        if rel_colorfunc !== nothing
            col_min, col_max = get_min_max_grid_value(rel_colorfunc)
            colmap = colmap / col_max
        else
            col_min = minimum(colmap)
            col_max = maximum(colmap)
            colmap = colmap / col_max
        end

        kwargs = (facecolors = get_cmap(colormap)(colmap), kwargs...)
    end

    if axis === nothing
        fig = figure(figsize=(8, 8))
        plot_surface(fx, fy, fz, rcount=100, ccount=100; kwargs...)
        xlabel("x")
        ylabel("y")
        zlabel("z")
        title(axis_title)
        gca().view_init(camera...)
        return fig
    end

    axis.plot_surface(fx, fy, fz, rcount=200, ccount=200; kwargs...)
    axis.set_ylabel("y")
    axis.set_zlabel("z")
    axis.set_xlabel("x")
    axis.set_title(axis_title)
    axis.view_init(camera[1], camera[2])
    return axis
end

function plot_parametric_wireframe(func; n_points=21, axis_title="", axis=nothing, camera=(25, -135), kwargs...)
    # Plot Data
    xs = ys = range(0, 1, length=n_points)

    F = [func([xi, yi]) for xi in xs, yi in xs]
    fx = [el[1] for el in F]
    fy = [el[2] for el in F]
    fz = [el[3] for el in F]

    if axis === nothing
        fig = figure(figsize=(8, 8))
        plot_wireframe(fx, fy, fz; kwargs...)
        xlabel("x")
        ylabel("y")
        zlabel("z")
        title(axis_title)
        gca().view_init(camera...)
        return fig
    end

    axis.plot_wireframe(fx, fy, fz; kwargs...)
    axis.set_ylabel("y")
    axis.set_zlabel("z")
    axis.set_xlabel("x")
    axis.set_title(axis_title)
    axis.view_init(camera[1], camera[2])
    return axis
end

function plot_diffeomorphism(ψ; n_points=15, axis=nothing)
    xs = 0.:(1/n_points):1.

    if axis === nothing
        fig = figure()

        # Plot horizontal lines
        for xi in xs
            line_x, line_y = [xi for yi in xs], [yi for yi in xs]
            warped_line_x = [ψ([xi, yi])[1] for yi in xs]
            warped_line_y = [ψ([xi, yi])[2] for yi in xs]

            # plot(line_x, line_y, c="lightgrey", lw=0.7)
            plot(warped_line_x, warped_line_y, c="k")
        end

        # Plot Warped Lines Vertically
        for yi in xs
            line_x, line_y = [xi for xi in xs], [yi for xi in xs]
            warped_line_x = [ψ([xi, yi])[1] for xi in xs]
            warped_line_y = [ψ([xi, yi])[2] for xi in xs]

            # plot(line_x, line_y, c="lightgrey", lw=0.7)
            plot(warped_line_x, warped_line_y, c="k")
        end

        return fig
    end

        # Plot horizontal lines
    for xi in xs
        line_x, line_y = [xi for yi in xs], [yi for yi in xs]
        warped_line_x = [ψ([xi, yi])[1] for yi in xs]
        warped_line_y = [ψ([xi, yi])[2] for yi in xs]

        # axis.plot(line_x, line_y, c="lightgrey", lw=0.7)
        axis.plot(warped_line_x, warped_line_y, c="k")
    end

    # Plot Warped Lines Vertically
    for yi in xs
        line_x, line_y = [xi for xi in xs], [yi for xi in xs]
        warped_line_x = [ψ([xi, yi])[1] for xi in xs]
        warped_line_y = [ψ([xi, yi])[2] for xi in xs]

        # axis.plot(line_x, line_y, c="lightgrey", lw=0.7)
        axis.plot(warped_line_x, warped_line_y, c="k")
    end
    return axis
end

"""
Retrieve color scale
"""
function get_min_max_grid_value(func, n_points=101)
    xs = 0:(1/n_points):1
    A = [func([xi, yi]) for xi in xs, yi in xs]

    return minimum(A), maximum(A)
end

function add_point(point, axis; kwargs...)
    if length(point) == 2
        axis.scatter([point[1]], [point[2]]; kwargs...)
    elseif length(point) == 3
        axis.scatter([point[1]], [point[2]], [point[3]]; kwargs...)
    else
        println("Dim must be 2 or 3")
    end
end

function plot_contourf(func, axis, n_points=101; kwargs...)
    # Plot Data
    xs = ys = range(0, 1, length=n_points)

    F = [func([xi, yi]) for xi in xs, yi in xs]
    axis.contourf(xs, xs, F; kwargs...)
    return axis
end


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

"""
    RectangleRegion{T,N}

A rectangle region in `N`-dimensional space with coordinates of type `T`.
"""
struct RectangleRegion{T,N} <: AbstractRegion{T,N}
  lowerleft::SVector{N,T}
  upperright::SVector{N,T}
end

RectangleRegion(lowerleft::NTuple{N,T}, upperright::NTuple{N,T}) where {N,T} =
  RectangleRegion{T,N}(lowerleft, upperright)

"""
    center(rectangle)

Return the center of the `rectangle`.
"""
center(r::RectangleRegion{T,N}) where {N,T} = @. (r.lowerleft + r.upperright) / 2

"""
    lowerleft(rectangle)

Return the lower left corner of the `rectangle`.
"""
lowerleft(r::RectangleRegion) = r.lowerleft

"""
    upperright(rectangle)

Return the upper right corner of the `rectangle`.
"""
upperright(r::RectangleRegion) = r.upperright

"""
    side(rectangle, i)

Return the `i`-th side of the `rectangle`.
"""
side(r::RectangleRegion, i::Int) = r.upperright[i] - r.lowerleft[i]

"""
    sides(rectangle)

Return all the sides of the `rectangle` as a tuple.
"""
sides(r::RectangleRegion{T,N}) where {N,T} = SVector{N,T}(ntuple(i->side(r, i), N))

"""
    diagonal(rectangle)

Return the diagonal of the `rectangle`.
"""
diagonal(r::RectangleRegion) = norm(u-l for (l,u) in zip(r.lowerleft, r.upperright))

"""
    volume(rectangle)

Return the volume of the `rectangle`.
"""
volume(r::RectangleRegion) = prod(u-l for (l,u) in zip(r.lowerleft, r.upperright))

ext: jl


text:
# Top level test script for Stan.jl
using ClassicalLaminateTheory
using Base.Test

println("Running tests for ClassicalLaminateTheory-j0.3/4-v1.0.0:")

tests = [
  "test_clt1.jl",
  "test_clt2.jl",
  "test_clt3.jl",
  "test_clt4.jl"
]

println("Running tests:")

for my_test in tests
    println("\n  * $(my_test) *")
    include(my_test)
end

println()


text:
mutable struct SendOptions
  isSSL::Bool
  username::String
  passwd::String
  verbose::Bool
end

function SendOptions(; isSSL::Bool = false, username::AbstractString = "",
                     passwd::AbstractString = "", verbose::Bool = false, kwargs...)
  kwargs = Dict(kwargs)
  if get(kwargs, :blocking, nothing) ≠ nothing
    @warn "options `blocking` is deprecated, blocking behaviour is default now, " *
          "use `@async send(...)` for non-blocking style."
    pop!(kwargs, :blocking)
  end
  length(keys(kwargs)) ≠ 0 && throw(MethodError("got unsupported keyword arguments"))
  SendOptions(isSSL, String(username), String(passwd), verbose)
end

function Base.show(io::IO, o::SendOptions)
  println(io, "SSL:      ", o.isSSL)
  print(  io, "verbose:  ", o.verbose)
  !isempty(o.username) && print(io, "\nusername: ", o.username)
end

mutable struct SendResponse
  body::IOBuffer
  code::Int
  total_time::Float64

  SendResponse() = new(IOBuffer(), 0, 0.0)
end

function Base.show(io::IO, o::SendResponse)
  println(io, "Return Code: ", o.code)
  println(io, "Time:        ", o.total_time)
  print(io,   "Response:    ", String(take!(o.body)))
end

mutable struct ReadData{T<:IO}
  typ::Symbol
  src::T
  str::AbstractString
  offset::Csize_t
  sz::Csize_t
end

ReadData() = ReadData{IOBuffer}(:undefined, IOBuffer(), "", 0, 0)
ReadData(io::T) where {T<:IO} = ReadData{T}(:io, io, "", 0, 0)
ReadData(io::IOBuffer)        = ReadData{IOBuffer}(:io, io, "", 0, io.size)

mutable struct ConnContext
  curl::Ptr{CURL}  # CURL handle
  url::String
  rd::ReadData
  resp::SendResponse
  options::SendOptions
  close_ostream::Bool
  bytes_recd::Int
  finalizer::Vector{Function}
end

function ConnContext(; curl = curl_easy_init(),
                     url::String = "",
                     rd::ReadData = ReadData(),
                     resp::SendResponse = SendResponse(),
                     options::SendOptions = SendOptions())
  curl == C_NULL && throw("curl_easy_init() failed")

  ctxt = ConnContext(curl, url, rd, resp, options, false, 0, Function[])

  @ce_curl curl_easy_setopt curl CURLOPT_URL url
  @ce_curl curl_easy_setopt curl CURLOPT_WRITEFUNCTION c_curl_write_cb
  @ce_curl curl_easy_setopt curl CURLOPT_WRITEDATA ctxt
  @ce_curl curl_easy_setopt curl CURLOPT_READFUNCTION c_curl_read_cb
  @ce_curl curl_easy_setopt curl CURLOPT_READDATA ctxt
  @ce_curl curl_easy_setopt curl CURLOPT_UPLOAD 1

  if options.isSSL
    @ce_curl curl_easy_setopt curl CURLOPT_USE_SSL CURLUSESSL_ALL
  end

  if !isempty(options.username)
    @ce_curl curl_easy_setopt curl CURLOPT_USERNAME options.username
    @ce_curl curl_easy_setopt curl CURLOPT_PASSWORD options.passwd
  end

  if options.verbose
    @ce_curl curl_easy_setopt curl CURLOPT_VERBOSE 1
  end

  ctxt
end

function setopt!(ctxt::ConnContext, opt, val)
  @ce_curl curl_easy_setopt ctxt.curl opt val
  ctxt
end

setmail_from!(ctxt::ConnContext, from::String) =
  setopt!(ctxt, CURLOPT_MAIL_FROM, from)

function setmail_rcpt!(ctxt::ConnContext, R::Vector{String})
  R′ = foldl(curl_slist_append, R, init = C_NULL)
  R′ == C_NULL && error("mail rcpts invalid")
  setopt!(ctxt, CURLOPT_MAIL_RCPT, R′)
  push!(ctxt.finalizer, () -> curl_slist_free_all(R′))
  ctxt
end

function connect(ctxt::ConnContext)
  @ce_curl curl_easy_perform ctxt.curl
  ctxt
end

cleanup!(::Nothing) = nothing
function cleanup!(ctxt::ConnContext)
  curl = ctxt.curl
  curl ≠ C_NULL && curl_easy_cleanup(curl)

  for f ∈ ctxt.finalizer
    f()
  end
  empty!(ctxt.finalizer)

  if ctxt.close_ostream
    close(ctxt.resp.body)
    ctxt.resp.body = nothing
    ctxt.close_ostream = false
  end
end

function getresponse!(ctxt::ConnContext)
  code = Array{Int}(undef, 1)
  @ce_curl curl_easy_getinfo ctxt.curl CURLINFO_RESPONSE_CODE code

  total_time = Array{Float64}(undef, 1)
  @ce_curl curl_easy_getinfo ctxt.curl CURLINFO_TOTAL_TIME total_time

  ctxt.resp.code = code[1]
  ctxt.resp.total_time = total_time[1]
end


text:
module Inference

using ..Core, ..Utilities
using DynamicPPL: Metadata, _tail, VarInfo, TypedVarInfo,
    islinked, invlink!, getlogp, tonamedtuple, VarName, getsym, vectorize,
    settrans!, _getvns, getdist, split_var_str, CACHERESET, AbstractSampler,
    Model, runmodel!, Sampler, SampleFromPrior, SampleFromUniform,
    Selector, AbstractSamplerState, DefaultContext, PriorContext,
    LikelihoodContext, MiniBatchContext, set_flag!, unset_flag!
using Distributions, Libtask, Bijectors
using ProgressMeter, LinearAlgebra
using ..Turing: PROGRESS, NamedDist, NoDist, Turing
using StatsFuns: logsumexp
using Random: GLOBAL_RNG, AbstractRNG, randexp
using AbstractMCMC, DynamicPPL

import MCMCChains: Chains
import AdvancedHMC; const AHMC = AdvancedHMC
import ..Core: getchunksize, getADtype
import AbstractMCMC: AbstractTransition, sample, step!, sample_init!,
    transitions_init, sample_end!, AbstractSampler, transition_type,
    callback, init_callback, AbstractCallback, psample
import DynamicPPL: tilde, dot_tilde, getspace, get_matching_type

export  InferenceAlgorithm,
        Hamiltonian,
        AbstractGibbs,
        GibbsComponent,
        StaticHamiltonian,
        AdaptiveHamiltonian,
        SampleFromUniform,
        SampleFromPrior,
        MH,
        ESS,
        Gibbs,      # classic sampling
        HMC,
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,       # Hamiltonian-like sampling
        DynamicNUTS,
        IS,
        SMC,
        CSMC,
        PG,
        PIMH,
        PMMH,
        IPMCMC,  # particle-based sampling
        assume,
        observe,
        step,
        WelfordVar,
        WelfordCovar,
        NaiveCovar,
        get_var,
        get_covar,
        add_sample!,
        reset!,
        step!,
        resume

#######################
# Sampler abstraction #
#######################
abstract type AbstractAdapter end
abstract type InferenceAlgorithm end
abstract type ParticleInference <: InferenceAlgorithm end
abstract type Hamiltonian{AD} <: InferenceAlgorithm end
abstract type StaticHamiltonian{AD} <: Hamiltonian{AD} end
abstract type AdaptiveHamiltonian{AD} <: Hamiltonian{AD} end

getchunksize(::Type{<:Hamiltonian{AD}}) where AD = getchunksize(AD)
getADtype(::Hamiltonian{AD}) where AD = AD

"""
    mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)

Decide if a proposal ``x'`` with log probability ``\\log p(x') = logp_proposal`` and
log proposal ratio ``\\log k(x', x) - \\log k(x, x') = log_proposal_ratio`` in a
Metropolis-Hastings algorithm with Markov kernel ``k(x_t, x_{t+1})`` and current state
``x`` with log probability ``\\log p(x) = logp_current`` is accepted by evaluating the
Metropolis-Hastings acceptance criterion
```math
\\log U \\leq \\log p(x') - \\log p(x) + \\log k(x', x) - \\log k(x, x')
```
for a uniform random number ``U \\in [0, 1)``.
"""
function mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)
    # replacing log(rand()) with -randexp() yields test errors
    return log(rand()) + logp_current ≤ logp_proposal + log_proposal_ratio
end

######################
# Default Transition #
######################

struct Transition{T, F<:AbstractFloat} <: AbstractTransition
    θ  :: T
    lp :: F
end

function Transition(spl::Sampler, nt::NamedTuple=NamedTuple())
    theta = merge(tonamedtuple(spl.state.vi), nt)
    lp = getlogp(spl.state.vi)
    return Transition{typeof(theta), typeof(lp)}(theta, lp)
end

function additional_parameters(::Type{<:Transition})
    return [:lp]
end

##########################################
# Internal variable names for MCMCChains #
##########################################

const TURING_INTERNAL_VARS = (internals = [
    "elapsed",
    "eval_num",
    "lf_eps",
    "lp",
    "weight",
    "le",
    "acceptance_rate",
    "hamiltonian_energy",
    "hamiltonian_energy_error",
    "max_hamiltonian_energy_error",
    "is_accept",
    "log_density",
    "n_steps",
    "numerical_error",
    "step_size",
    "nom_step_size",
    "tree_depth",
    "is_adapt",
],)

#########################################
# Default definitions for the interface #
#########################################

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    kwargs...
)
    return sample(rng, model, Sampler(alg, model), N; progress=PROGRESS[], kwargs...)
end

function AbstractMCMC.sample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    resume_from=nothing,
    kwargs...
)
    if resume_from === nothing
        return sample(model, Sampler(alg, model), N; progress=PROGRESS[], kwargs...)
    else
        return resume(resume_from, N)
    end
end

function AbstractMCMC.psample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer,
    n_chains::Integer;
    kwargs...
)
    return psample(GLOBAL_RNG, model, alg, N, n_chains; progress=false, kwargs...)
end

function AbstractMCMC.psample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer,
    n_chains::Integer;
    kwargs...
)
    return psample(rng, model, Sampler(alg, model), N, n_chains; progress=false, kwargs...)
end

function AbstractMCMC.sample_init!(
    ::AbstractRNG,
    model::Model,
    spl::Sampler,
    N::Integer;
    kwargs...
)
    # Resume the sampler.
    set_resume!(spl; kwargs...)

    # Set the parameters to a starting value.
    initialize_parameters!(spl; kwargs...)
end

function AbstractMCMC.sample_end!(
    ::AbstractRNG,
    ::Model,
    ::AbstractSampler,
    ::Integer,
    ::Vector{<:AbstractTransition};
    kwargs...
)
    # Silence the default API function.
end

function initialize_parameters!(
    spl::AbstractSampler;
    init_theta::Union{Nothing,Vector}=nothing,
    verbose::Bool=false,
    kwargs...
)
    # Get `init_theta`
    if init_theta !== nothing
        verbose && @info "Using passed-in initial variable values" init_theta
        # Convert individual numbers to length 1 vector; `ismissing(v)` is needed as `size(missing)` is undefined`
        init_theta = [ismissing(v) || size(v) == () ? [v] : v for v in init_theta]
        # Flatten `init_theta`
        init_theta_flat = foldl(vcat, map(vec, init_theta))
        # Create a mask to inidicate which values are not missing
        theta_mask = map(x -> !ismissing(x), init_theta_flat)
        # Get all values
        theta = spl.state.vi[spl]
        @assert length(theta) == length(init_theta_flat) "Provided initial value doesn't match the dimension of the model"
        # Update those which are provided (i.e. not missing)
        theta[theta_mask] .= init_theta_flat[theta_mask]
        # Update in `vi`
        spl.state.vi[spl] = theta
    end
end

##########################
# Chain making utilities #
##########################

function _params_to_array(ts::Vector{<:AbstractTransition}, spl::Sampler)
    names_set = Set{String}()
    # Extract the parameter names and values from each transition.
    dicts = map(ts) do t
        nms, vs = flatten_namedtuple(t.θ)
        for nm in nms
            push!(names_set, nm)
        end
        # Convert the names and values to a single dictionary.
        return Dict(nms[j] => vs[j] for j in 1:length(vs))
    end
    names = collect(names_set)
    vals = [get(dicts[i], key, missing) for i in eachindex(dicts),
        (j, key) in enumerate(names)]

    return names, vals
end

function flatten_namedtuple(nt::NamedTuple)
    names_vals = mapreduce(vcat, keys(nt)) do k
        v = nt[k]
        if length(v) == 1
            return [(string(k), v)]
        else
            return mapreduce(vcat, zip(v[1], v[2])) do (vnval, vn)
                return collect(FlattenIterator(vn, vnval))
            end
        end
    end
    return [vn[1] for vn in names_vals], [vn[2] for vn in names_vals]
end

function get_transition_extras(ts::Vector{<:AbstractTransition})
    # Get the extra field names from the sampler state type.
    # This handles things like :lp or :weight.
    extra_params = additional_parameters(eltype(ts))

    # Get the values of the extra parameters.
    local extra_names
    all_vals = []

    # Iterate through each transition.
    for t in ts
        extra_names = String[]
        vals = []

        # Iterate through each of the additional field names
        # in the struct.
        for p in extra_params
            # Check whether the field contains a NamedTuple,
            # in which case we need to iterate through each
            # key/value pair.
            prop = getproperty(t, p)
            if prop isa NamedTuple
                for (k, v) in pairs(prop)
                    push!(extra_names, string(k))
                    push!(vals, v)
                end
            else
                push!(extra_names, string(p))
                push!(vals, prop)
            end
        end
        push!(all_vals, vals)
    end

    # Convert the vector-of-vectors to a matrix.
    valmat = [all_vals[i][j] for i in 1:length(ts), j in 1:length(all_vals[1])]

    return extra_names, valmat
end

# Default Chains constructor.
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Sampler,
    N::Integer,
    ts::Vector{<:AbstractTransition};
    discard_adapt::Bool=true,
    save_state=true,
    kwargs...
)
    # Check if we have adaptation samples.
    if discard_adapt && :n_adapts in fieldnames(typeof(spl.alg))
        ts = ts[(spl.alg.n_adapts+1):end]
    end

    # Convert transitions to array format.
    # Also retrieve the variable names.
    nms, vals = _params_to_array(ts, spl)

    # Get the values of the extra parameters in each Transition struct.
    extra_params, extra_values = get_transition_extras(ts)

    # Extract names & construct param array.
    nms = [nms; extra_params]
    parray = hcat(vals, extra_values)

    # If the state field has average_logevidence or final_logevidence, grab that.
    le = missing
    if :average_logevidence in fieldnames(typeof(spl.state))
        le = getproperty(spl.state, :average_logevidence)
    elseif :final_logevidence in fieldnames(typeof(spl.state))
        le = getproperty(spl.state, :final_logevidence)
    end

    # Check whether to invlink! the varinfo
    if islinked(spl.state.vi, spl)
        invlink!(spl.state.vi, spl)
    end

    # Set up the info tuple.
    if save_state
        info = (range = rng, model = model, spl = spl, vi = spl.state.vi)
    else
        info = NamedTuple()
    end

    # Chain construction.
    return Chains(
        parray,
        string.(nms),
        deepcopy(TURING_INTERNAL_VARS);
        evidence=le,
        info=info,
        sorted=true
    )
end

function save(c::Chains, spl::AbstractSampler, model, vi, samples)
    nt = NamedTuple{(:spl, :model, :vi, :samples)}((spl, model, deepcopy(vi), samples))
    return setinfo(c, merge(nt, c.info))
end

function resume(c::Chains, n_iter::Int; kwargs...)
    @assert !isempty(c.info) "[Turing] cannot resume from a chain without state info"

    # Sample a new chain.
    newchain = sample(
        c.info[:range],
        c.info[:model],
        c.info[:spl],
        n_iter;
        resume_from=c,
        reuse_spl_n=n_iter,
        kwargs...
    )

    # Stick the new samples at the end of the old chain.
    return vcat(c, newchain)
end

function set_resume!(
    s::Sampler;
    resume_from::Union{Chains, Nothing}=nothing,
    kwargs...
)
    # If we're resuming, grab the sampler info.
    if resume_from !== nothing
        s = resume_from.info[:spl]
    end
end

#########################
# Default sampler state #
#########################

"""
A blank `AbstractSamplerState` that contains only `VarInfo` information.
"""
mutable struct SamplerState{VIType<:VarInfo} <: AbstractSamplerState
    vi :: VIType
end

#######################################
# Concrete algorithm implementations. #
#######################################

include("ess.jl")
include("hmc.jl")
include("mh.jl")
include("is.jl")
include("AdvancedSMC.jl")
include("gibbs.jl")
include("../contrib/inference/sghmc.jl")
include("../contrib/inference/AdvancedSMCExtensions.jl")

################
# Typing tools #
################

for alg in (:SMC, :PG, :PMMH, :IPMCMC, :MH, :IS, :ESS, :Gibbs)
    @eval getspace(::$alg{space}) where {space} = space
end
for alg in (:HMC, :HMCDA, :NUTS, :SGLD, :SGHMC)
    @eval getspace(::$alg{<:Any, space}) where {space} = space
end

floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
floatof(::Type) = Real # fallback if type inference failed

function get_matching_type(
    spl::AbstractSampler,
    vi::VarInfo,
    ::Type{T},
) where {T}
    return T
end
function get_matching_type(
    spl::AbstractSampler,
    vi::VarInfo,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{<:Union{Missing, AbstractFloat}},
)
    return Union{Missing, floatof(eltype(vi, spl))}
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{TV},
) where {T, N, TV <: Array{T, N}}
    return Array{get_matching_type(spl, vi, T), N}
end
function get_matching_type(
    spl::Sampler{<:Union{PG, SMC}},
    vi::VarInfo,
    ::Type{TV},
) where {T, N, TV <: Array{T, N}}
    return TArray{T, N}
end

## Fallback functions

alg_str(spl::Sampler) = string(nameof(typeof(spl.alg)))
transition_type(spl::Sampler) = typeof(Transition(spl))

# utility funcs for querying sampler information
require_gradient(spl::Sampler) = false
require_particles(spl::Sampler) = false

_getindex(x, inds::Tuple) = _getindex(x[first(inds)...], Base.tail(inds))
_getindex(x, inds::Tuple{}) = x

# assume
function tilde(ctx::DefaultContext, sampler, right, vn::VarName, _, vi)
    return _tilde(sampler, right, vn, vi)
end
function tilde(ctx::PriorContext, sampler, right, vn::VarName, inds, vi)
    if ctx.vars !== nothing
        vi[vn] = vectorize(right, _getindex(getfield(ctx.vars, getsym(vn)), inds))
        settrans!(vi, false, vn)
    end
    return _tilde(sampler, right, vn, vi)
end
function tilde(ctx::LikelihoodContext, sampler, right, vn::VarName, inds, vi)
    if ctx.vars !== nothing
        vi[vn] = vectorize(right, _getindex(getfield(ctx.vars, getsym(vn)), inds))
        settrans!(vi, false, vn)
    end
    return _tilde(sampler, NoDist(right), vn, vi)
end
function tilde(ctx::MiniBatchContext, sampler, right, left::VarName, inds, vi)
    return tilde(ctx.ctx, sampler, right, left, inds, vi)
end

function _tilde(sampler, right, vn::VarName, vi)
    return Turing.assume(sampler, right, vn, vi)
end
function _tilde(sampler, right::NamedDist, vn::VarName, vi)
    name = right.name
    if name isa String
        sym_str, inds = split_var_str(name, String)
        sym = Symbol(sym_str)
        vn = VarName{sym}(inds)
    elseif name isa Symbol
        vn = VarName{name}("")
    elseif name isa VarName
        vn = name
    else
        throw("Unsupported variable name. Please use either a string, symbol or VarName.")
    end
    return _tilde(sampler, right.dist, vn, vi)
end

# observe
function tilde(ctx::DefaultContext, sampler, right, left, vi)
    return _tilde(sampler, right, left, vi)
end
function tilde(ctx::PriorContext, sampler, right, left, vi)
    return 0
end
function tilde(ctx::LikelihoodContext, sampler, right, left, vi)
    return _tilde(sampler, right, left, vi)
end
function tilde(ctx::MiniBatchContext, sampler, right, left, vi)
    return ctx.loglike_scalar * tilde(ctx.ctx, sampler, right, left, vi)
end

_tilde(sampler, right, left, vi) = Turing.observe(sampler, right, left, vi)

function assume(spl::Sampler, dist)
    error("Turing.assume: unmanaged inference algorithm: $(typeof(spl))")
end

function observe(spl::Sampler, weight)
    error("Turing.observe: unmanaged inference algorithm: $(typeof(spl))")
end

function assume(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dist::Distribution,
    vn::VarName,
    vi::VarInfo,
)
    if haskey(vi, vn)
        if is_flagged(vi, vn, "del")
            unset_flag!(vi, vn, "del")
            r = spl isa SampleFromUniform ? init(dist) : rand(dist)
            vi[vn] = vectorize(dist, r)
            setorder!(vi, vn, vi.num_produce)
        else
        r = vi[vn]
        end
    else
        r = isa(spl, SampleFromUniform) ? init(dist) : rand(dist)
        push!(vi, vn, r, dist, spl)
    end
    # NOTE: The importance weight is not correctly computed here because
    #       r is genereated from some uniform distribution which is different from the prior
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))

    return r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function observe(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dist::Distribution,
    value,
    vi::VarInfo,
)
    vi.num_produce += one(vi.num_produce)
    return logpdf(dist, value)
end

# .~ functions

# assume
function dot_tilde(ctx::DefaultContext, sampler, right, left, vn::VarName, _, vi)
    vns, dist = get_vns_and_dist(right, left, vn)
    return _dot_tilde(sampler, dist, left, vns, vi)
end
function dot_tilde(
    ctx::LikelihoodContext,
    sampler,
    right,
    left,
    vn::VarName,
    inds,
    vi,
)
    if ctx.vars !== nothing
        var = _getindex(getfield(ctx.vars, getsym(vn)), inds)
        vns, dist = get_vns_and_dist(right, var, vn)
        set_val!(vi, vns, dist, var)
        settrans!.(Ref(vi), false, vns)
    else
        vns, dist = get_vns_and_dist(right, left, vn)
    end
    return _dot_tilde(sampler, NoDist(dist), left, vns, vi)
end
function dot_tilde(ctx::MiniBatchContext, sampler, right, left, vn::VarName, inds, vi)
    return dot_tilde(ctx.ctx, sampler, right, left, vn, inds, vi)
end
function dot_tilde(
    ctx::PriorContext,
    sampler,
    right,
    left,
    vn::VarName,
    inds,
    vi,
)
    if ctx.vars !== nothing
        var = _getindex(getfield(ctx.vars, getsym(vn)), inds)
        vns, dist = get_vns_and_dist(right, var, vn)
        set_val!(vi, vns, dist, var)
        settrans!.(Ref(vi), false, vns)
    else
        vns, dist = get_vns_and_dist(right, left, vn)
    end
    return _dot_tilde(sampler, dist, left, vns, vi)
end

function get_vns_and_dist(dist::NamedDist, var, vn::VarName)
    name = dist.name
    if name isa String
        sym_str, inds = split_var_str(name, String)
        sym = Symbol(sym_str)
        vn = VarName{sym}(inds)
    elseif name isa Symbol
        vn = VarName{name}("")
    elseif name isa VarName
        vn = name
    else
        throw("Unsupported variable name. Please use either a string, symbol or VarName.")
    end
    return get_vns_and_dist(dist.dist, var, vn)
end
function get_vns_and_dist(dist::MultivariateDistribution, var::AbstractMatrix, vn::VarName)
    getvn = i -> VarName(vn, vn.indexing * "[Colon(),$i]")
    return getvn.(1:size(var, 2)), dist
end
function get_vns_and_dist(
    dist::Union{Distribution, AbstractArray{<:Distribution}},
    var::AbstractArray,
    vn::VarName
)
    getvn = ind -> VarName(vn, vn.indexing * "[" * join(Tuple(ind), ",") * "]")
    return getvn.(CartesianIndices(var)), dist
end

function _dot_tilde(sampler, right, left, vns::AbstractArray{<:VarName}, vi)
    return dot_assume(sampler, right, vns, left, vi)
end

# Ambiguity error when not sure to use Distributions convention or Julia broadcasting semantics
function _dot_tilde(
    sampler::AbstractSampler,
    right::Union{MultivariateDistribution, AbstractVector{<:MultivariateDistribution}},
    left::AbstractMatrix{>:AbstractVector},
    vn::AbstractVector{<:VarName},
    vi::VarInfo,
)
    throw(ambiguity_error_msg())
end

function dot_assume(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dist::MultivariateDistribution,
    vns::AbstractVector{<:VarName},
    var::AbstractMatrix,
    vi::VarInfo,
)
    @assert length(dist) == size(var, 1)
    r = get_and_set_val!(vi, vns, dist, spl)
    lp = sum(logpdf_with_trans(dist, r, istrans(vi, vns[1])))
    var .= r
    return var, lp
end
function dot_assume(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    vns::AbstractArray{<:VarName},
    var::AbstractArray,
    vi::VarInfo,
)
    r = get_and_set_val!(vi, vns, dists, spl)
    # Make sure `r` is not a matrix for multivariate distributions
    lp = sum(logpdf_with_trans.(dists, r, istrans(vi, vns[1])))
    var .= r
    return var, lp
end
function dot_assume(
    spl::Sampler,
    ::Any,
    ::AbstractArray{<:VarName},
    ::Any,
    ::VarInfo
)
    error("[Turing] $(alg_str(spl)) doesn't support vectorizing assume statement")
end

function get_and_set_val!(
    vi::VarInfo,
    vns::AbstractVector{<:VarName},
    dist::MultivariateDistribution,
    spl::AbstractSampler,
)
    n = length(vns)
    if haskey(vi, vns[1])
        if is_flagged(vi, vns[1], "del")
            unset_flag!(vi, vns[1], "del")
            r = spl isa SampleFromUniform ? init(dist, n) : rand(dist, n)
            for i in 1:n
                vn = vns[i]
                vi[vn] = vectorize(dist, r[:, i])
                setorder!(vi, vn, vi.num_produce)
            end
        else
        r = vi[vns]
        end
    else
        r = spl isa SampleFromUniform ? init(dist, n) : rand(dist, n)
        for i in 1:n
            push!(vi, vns[i], r[:,i], dist, spl)
        end
    end
    return r
end
function get_and_set_val!(
    vi::VarInfo,
    vns::AbstractArray{<:VarName},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    spl::AbstractSampler,
)
    if haskey(vi, vns[1])
        if is_flagged(vi, vns[1], "del")
            unset_flag!(vi, vns[1], "del")
            f = (vn, dist) -> spl isa SampleFromUniform ? init(dist) : rand(dist)
            r = f.(vns, dists)
            for i in eachindex(vns)
                vn = vns[i]
                dist = dists isa AbstractArray ? dists[i] : dists
                vi[vn] = vectorize(dist, r[i])
                setorder!(vi, vn, vi.num_produce)
            end
        else
        r = reshape(vi[vec(vns)], size(vns))
        end
    else
        f = (vn, dist) -> spl isa SampleFromUniform ? init(dist) : rand(dist)
        r = f.(vns, dists)
        push!.(Ref(vi), vns, r, dists, Ref(spl))
    end
    return r
end

function set_val!(
    vi::VarInfo,
    vns::AbstractVector{<:VarName},
    dist::MultivariateDistribution,
    val::AbstractMatrix,
)
    @assert size(val, 2) == length(vns)
    foreach(enumerate(vns)) do (i, vn)
        vi[vn] = val[:,i]
    end
    return val
end
function set_val!(
    vi::VarInfo,
    vns::AbstractArray{<:VarName},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    val::AbstractArray,
)
    @assert size(val) == size(vns)
    foreach(CartesianIndices(val)) do ind
        dist = dists isa AbstractArray ? dists[ind] : dists
        vi[vns[ind]] = vectorize(dist, val[ind])
    end
    return val
end

# observe
function dot_tilde(ctx::DefaultContext, sampler, right, left, vi)
    return _dot_tilde(sampler, right, left, vi)
end
function dot_tilde(ctx::PriorContext, sampler, right, left, vi)
    return 0
end
function dot_tilde(ctx::LikelihoodContext, sampler, right, left, vi)
    return _dot_tilde(sampler, right, left, vi)
end
function dot_tilde(ctx::MiniBatchContext, sampler, right, left, vi)
    return ctx.loglike_scalar * dot_tilde(ctx.ctx, sampler, right, left, left, vi)
end

function _dot_tilde(sampler, right, left::AbstractArray, vi)
    return dot_observe(sampler, right, left, vi)
end
# Ambiguity error when not sure to use Distributions convention or Julia broadcasting semantics
function _dot_tilde(
    sampler::AbstractSampler,
    right::Union{MultivariateDistribution, AbstractVector{<:MultivariateDistribution}},
    left::AbstractMatrix{>:AbstractVector},
    vi::VarInfo,
)
    throw(ambiguity_error_msg())
end

function dot_observe(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dist::MultivariateDistribution,
    value::AbstractMatrix,
    vi::VarInfo,
)
    vi.num_produce += one(vi.num_produce)
    Turing.DEBUG && @debug "dist = $dist"
    Turing.DEBUG && @debug "value = $value"
    return sum(logpdf(dist, value))
end
function dot_observe(
    spl::Union{SampleFromPrior, SampleFromUniform},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    value::AbstractArray,
    vi::VarInfo,
)
    vi.num_produce += one(vi.num_produce)
    Turing.DEBUG && @debug "dists = $dists"
    Turing.DEBUG && @debug "value = $value"
    return sum(logpdf.(dists, value))
end
function dot_observe(
    spl::Sampler,
    ::Any,
    ::Any,
    ::VarInfo,
)
    error("[Turing] $(alg_str(spl)) doesn't support vectorizing observe statement")
end

##############
# Utilities  #
##############

getspace(spl::Sampler) = getspace(spl.alg)
function ambiguity_error_msg()
    return "Ambiguous `lhs .~ rhs` or `@. lhs ~ rhs` syntax. The broadcasting can either be
    column-wise following the convention of Distributions.jl or element-wise following
    Julia's general broadcasting semantics. Please make sure that the element type of `lhs`
    is not a supertype of the support type of `AbstractVector` to eliminate ambiguity."
end

end # module


text:
#########
# Rules #
#########

abstract type Rule end

mutable struct Terminal <: Rule
  name::AbstractString
  value::AbstractString
  action
end
Terminal(name::AbstractString, value) = Terminal(name, string(value), no_action)

mutable struct ReferencedRule <: Rule
  name::AbstractString
  symbol::Symbol
  action
end
ReferencedRule(name::AbstractString, symbol::Symbol) = ReferencedRule(name, symbol, no_action)

mutable struct AndRule <: Rule
  name::AbstractString
  values::Array{Rule}
  action
end
AndRule(name::AbstractString, values::Array{Rule}) = AndRule(name, values, no_action)
AndRule(name::AbstractString, left::Rule, right::Rule) = AndRule(name, [left, right], no_action)

mutable struct OrRule <: Rule
  name::AbstractString
  values::Array{Rule}
  action
end
OrRule(name::AbstractString, rules::Array) = OrRule(name,rules,liftchild)
OrRule(name::AbstractString, left::OrRule, right::OrRule) = OrRule(name, append!(left.values, right.values), liftchild)
OrRule(name::AbstractString, left::OrRule, right::Rule) = OrRule(name, push!(left.values, right), liftchild)
OrRule(name::AbstractString, left::Rule, right::OrRule) = OrRule(name, [left, right], liftchild)
OrRule(name::AbstractString, left::Rule, right::Rule) = OrRule(name, [left, right], liftchild)

mutable struct OneOrMoreRule <: Rule
  name::AbstractString
  value::Rule
  action
end

mutable struct ZeroOrMoreRule <: Rule
  name::AbstractString
  value::Rule
  action
end

mutable struct MultipleRule <: Rule
  name::AbstractString
  value::Rule
  minCount::Int
  maxCount::Int
  action
end
MultipleRule(name::AbstractString, value::Rule, minCount::Int, maxCount::Int) = MultipleRule(name, value, minCount, maxCount, no_action)

mutable struct RegexRule <: Rule
  name::AbstractString
  value::Regex
  action
end
RegexRule(name::AbstractString, value::Regex) = RegexRule(name, value, no_action)

mutable struct OptionalRule <: Rule
  name::AbstractString
  value::Rule
  action
end

mutable struct LookAheadRule <: Rule
    name::AbstractString
    value::Rule
    action
end

mutable struct SuppressRule <: Rule
  name::AbstractString
  value::Rule
  action
end

mutable struct ListRule <: Rule
  name::AbstractString
  entry::Rule
  delim::Rule
  min::Int
  action
end
ListRule(name::AbstractString, entry::Rule, delim::Rule, min::Int=1) = ListRule(name, entry, delim, min, no_action)

mutable struct NotRule <: Rule
  name
  entry
  action
end
NotRule(name::AbstractString, entry::Rule) = NotRule(name, entry, no_action)

mutable struct EndOfFileRule <: Rule
  name::AbstractString
  action
end

"""
    EmptyRule(name,action)
is a rule that never consumes
"""
mutable struct EmptyRule <: Rule
  name
  action
end

mutable struct IntegerRule <: Rule
  name::AbstractString
  action
end

mutable struct FloatRule <: Rule
  name::AbstractString
  action
end

################
# Constructors #
################

for rule in [FloatRule, IntegerRule, EmptyRule, EndOfFileRule]
  eval(Meta.parse("$rule(name::AbstractString) = $rule(name, no_action)"))
end

for rule in [SuppressRule, LookAheadRule, OptionalRule, ZeroOrMoreRule, OneOrMoreRule]
  eval(Meta.parse("$rule(name::AbstractString, value::Rule) = $rule(name, value, no_action)"))
end

for rule in subtypes(Rule)
  eval(Meta.parse("$rule(args...) = $rule(\"\", args...)"))
end

########
# show #
########

function showRule(io::IO,name::AbstractString, def::AbstractString, action::AbstractString)
  namestring = ""
  if name != ""
    namestring = "$name: "
  end
  actionstring = ""
  if action != string(no_action)
    actionstring = " { $action }"
  end
  print(io, "$namestring$def$actionstring")
end
function showRuleDouble(io::IO,name::AbstractString, def::AbstractString, action::AbstractString)
  namestring = ""
  if name != ""
    namestring = "$name: "
  end
  actionstring = ""
  parenth = ""
  if action != string(no_action)
    parenth = "("
    actionstring = ") { $action }"
  end
  print(io, "$namestring$parenth$def$actionstring")
end

show(io::IO, t::Terminal) = showRule(io, t.name, "'$(t.value)'", string(t.action))
show(io::IO, rule::ReferencedRule) = showRule(io, rule.name, "$(rule.symbol)", string(rule.action))
show(io::IO, rule::OneOrMoreRule) = showRule(io, rule.name, "+($(string(rule.value)))", string(rule.action));
show(io::IO, rule::OptionalRule) = showRule(io, rule.name, "?($(string(rule.value)))", string(rule.action));
show(io::IO, rule::ZeroOrMoreRule) = showRule(io, rule.name, "*($(string(rule.value)))", string(rule.action));
show(io::IO, rule::MultipleRule) = showRule(io, rule.name, "($(rule.value)){$(rule.minCount), $(rule.maxCount)}", string(rule.value));
show(io::IO, rule::RegexRule) = showRule(io, rule.name, "r($(rule.value.pattern))", string(rule.action))
show(io::IO, rule::SuppressRule) = showRule(io, rule.name, "-($(string(rule.value)))", string(rule.action))
show(io::IO, rule::NotRule) = showRule(io, rule.name, "!($(rule.entry))", string(rule.action))
function show(io::IO, rule::AndRule)
  values = [string(r) for r in rule.values]
  joinedValues = join(values, " & ")
  showRuleDouble(io, rule.name, joinedValues, string(rule.action))
end
function show(io::IO, rule::OrRule)
  values = [string(r) for r in rule.values]
  joinedValues = join(values, " | ")
  showRuleDouble(io,rule.name, joinedValues, string(rule.action))
end


text:
"""
```
macro dimension(symb, abbr, name)
```

Creates new dimensions. `name` will be used like an identifier in the type
parameter for a [`Unitful.Dimension`](@ref) object. `symb` will be a symbol
defined in the namespace from which this macro is called that is bound to a
[`Unitful.Dimensions`](@ref) object. For most intents and purposes it is this
object that the user would manipulate in doing dimensional analysis. The symbol
is not exported.

This macro extends [`Unitful.abbr`](@ref) to display the new dimension in an
abbreviated format using the string `abbr`.

Finally, type aliases are created that allow the user to dispatch on
[`Unitful.Quantity`](@ref) and [`Unitful.Units`](@ref) objects of the newly
defined dimension. The type alias for quantities is simply given by `name`,
and the type alias for units is given by `name*"Unit"`, e.g. `LengthUnit`.

Usage example: `@dimension 𝐋 "L" Length` (see `src/Defaults.jl`.)
"""
macro dimension(symb, abbr, name)
    s = Symbol(symb)
    x = Expr(:quote, name)
    uname = Symbol(name,"Unit")
    esc(quote
        Unitful.abbr(::Unitful.Dimension{$x}) = $abbr
        const $s = Unitful.Dimensions{(Unitful.Dimension{$x}(1),)}()
        typealias $(name) DimensionedQuantity{typeof($s)}
        typealias $(uname) DimensionedUnits{typeof($s)}
    end)
end

"""
```
macro derived_dimension(name, dims)
```

Creates type aliases to allow dispatch on [`Unitful.Quantity`](@ref) and
[`Unitful.Units`](@ref) objects of a derived dimension, like area, which is just
length squared. The type aliases are not exported.

`dims` is a [`Unitful.Dimensions`](@ref) object.

Usage examples:

- `@derived_dimension Area 𝐋^2` gives `Area` and `AreaUnit` type aliases
- `@derived_dimension Speed 𝐋/𝐓` gives `Speed` and `SpeedUnit` type aliases
"""
macro derived_dimension(name, dims)
    uname = Symbol(name,"Unit")
    esc(quote
        typealias ($name) DimensionedQuantity{typeof($dims)}
        typealias ($uname) DimensionedUnits{typeof($dims)}
    end)
end

"""
```
macro refunit(symb, name, abbr, dimension, tf)
```

Define a reference unit, typically SI. Rather than define
conversion factors between each and every unit of a given dimension, conversion
factors are given between each unit and a reference unit, defined by this macro.

This macro extends [`Unitful.abbr`](@ref) so that the reference unit can be
displayed in an abbreviated format. If `tf == true`, this macro generates symbols
for every power of ten of the unit, using the standard SI prefixes. A `dimension`
must be given ([`Unitful.Dimensions`](@ref) object) that specifies the dimension
of the reference unit.

Usage example: `@refunit m "m" Meter 𝐋 true`

This example will generate `km`, `m`, `cm`, ...
"""
macro refunit(symb, abbr, name, dimension, tf)
    x = Expr(:quote, name)
    esc(quote
        Unitful.abbr(::Unitful.Unit{$x}) = $abbr
        Unitful.dimension(y::Unitful.Unit{$x}) = $dimension^y.power
        Unitful.basefactor(y::Unitful.Unit{$x}) = (1.0, 1)
        if $tf
            Unitful.@prefixed_unit_symbols $symb $name
        else
            Unitful.@unit_symbols $symb $name
        end
    end)
end

"""
```
macro preferunit(unit)
```

This macro specifies the default unit for promotion for a given dimension,
which is inferred from the given unit.

Usage example: `@preferunit kg`
"""
macro preferunit(unit)
    quote
        dim = dimension($unit)
        if length(typeof(dim).parameters[1]) > 1
            error("@prefer can only be used with a unit that has a pure ",
            "dimension, like 𝐋 or 𝐓 but not 𝐋/𝐓.")
        end
        if length(typeof(dim).parameters[1]) == 1 &&
            typeof(dim).parameters[1][1].power != 1
            error("@prefer cannot handle powers of pure dimensions except 1. ",
            "For instance, it should not be used with units of dimension 𝐋^2.")
        end
        Unitful.dim2refunits(y::typeof(typeof(dim).parameters[1][1])) =
            $unit^y.power
    end
end

# Generated to force a concrete result type.
@generated function dim2refunits(x::Dimensions)
    dim = x.parameters[1]
    y = mapreduce(dim2refunits, *, NoUnits, dim)
    :($y)
end

"""
```
macro unit(symb,abbr,name,equals,tf)
```

Define a unit. Rather than specifying a dimension like in [`@refunit`](@ref),
`equals` should be a [`Unitful.Quantity`](@ref) equal to one of the unit being
defined. If `tf == true`, symbols will be made for each power-of-ten prefix.

Usage example: `@unit mi "mi" Mile (201168//125)*m false`

This example will *not* generate `kmi` (kilomiles).
"""
macro unit(symb,abbr,name,equals,tf)
    # name is a symbol
    # abbr is a string
    x = Expr(:quote, name)
    quote
        inex, ex = Unitful.basefactor(Unitful.unit($(esc(equals))))
        t = Unitful.tensfactor(Unitful.unit($(esc(equals))))
        eq = ($(esc(equals)))/Unitful.unit($(esc(equals)))
        Base.isa(eq, Base.Integer) || Base.isa(eq, Base.Rational) ?
             (ex *= eq) : (inex *= eq)
        Unitful.abbr(::Unitful.Unit{$(esc(x))}) = $abbr
        Unitful.dimension(y::Unitful.Unit{$(esc(x))}) =
            Unitful.dimension($(esc(equals)))^y.power
        Unitful.basefactor(y::Unitful.Unit{$(esc(x))}) =
            Unitful.basefactorhelper(inex, ex, t, y.power)
        if $tf
            Unitful.@prefixed_unit_symbols $(esc(symb)) $(esc(name))
        else
            Unitful.@unit_symbols $(esc(symb)) $(esc(name))
        end
    end
end

"""
```
macro prefixed_unit_symbols(symb,name)
```

Not called directly by the user. Given a unit symbol and a unit's name,
will define units for each possible SI power-of-ten prefix on that unit.

Example: `@prefixed_unit_symbols m Meter` results in nm, cm, m, km, ...
all getting defined in the calling namespace.
"""
macro prefixed_unit_symbols(symb,name)
    expr = Expr(:block)

    z = Expr(:quote, name)
    for (k,v) in prefixdict
        s = Symbol(v,symb)
        u = Unitful.Unit{name}(k,1//1)
        ea = esc(quote
            const $s = Unitful.Units{($u,),typeof(dimension($u))}()
        end)
        push!(expr.args, ea)
    end

    expr
end

"""
```
macro unit_symbols(symb,name)
```

Not called directly by the user. Given a unit symbol and a unit's name,
will define units without SI power-of-ten prefixes.

Example: `@unit_symbols ft Foot` results in `ft` getting defined but not `kft`.
"""
macro unit_symbols(symb,name)
    s = Symbol(symb)
    z = Expr(:quote, name)
    u = Unitful.Unit{name}(0,1//1)
    esc(quote
        const $s = Unitful.Units{($u,),typeof(dimension($u))}()
    end)
end

"""
```
defaults()
```

Includes the file `deps/Defaults.jl` from the Unitful package. This results in
common units and dimensions being generated in the `Unitful` module.
"""
function defaults()
    defpath = joinpath(Pkg.dir("Unitful"),"deps","Defaults.jl")
    include(defpath)
end

"""
```
macro u_str(unit)
```

String macro to easily recall units, dimensions, or quantities defined in the
Unitful module, which does not export such things to avoid namespace pollution.
Note that for now, what goes inside must be parsable as a valid Julia expression.
In other words, u"N m" will fail if you intended to write u"N*m".

Examples:

```jldoctest
julia> 1.0u"m/s"
1.0 m s^-1

julia> 1.0u"N*m"
1.0 m N

julia> typeof(1.0u"m/s")
Quantity{Float64, Dimensions:{𝐋 𝐓^-1}, Units:{m s^-1}}

julia> u"ħ"
1.0545718001391127e-34 J s
```
"""
macro u_str(unit)
    ex = parse(unit)
    replace_value(ex)
end

const allowed_funcs = [:*, :/, :^, :sqrt, :√, :+, :-, ://]
function replace_value(ex::Expr)
    ex.head != :call && error("$(ex.head) != :call")
    ex.args[1] in allowed_funcs ||
        error("""$(ex.args[1]) is not a valid function call when parsing a unit.
         Only the following functions are allowed: $allowed_funcs""")
    for i=2:length(ex.args)
        if typeof(ex.args[i])==Symbol || typeof(ex.args[i])==Expr
            ex.args[i]=replace_value(ex.args[i])
        end
    end
    ex
end

replace_value(sym::Symbol) = :(ustrcheck($sym))
ustrcheck(x::Unitlike) = x
ustrcheck(x::Quantity) = x
ustrcheck(x) = error("Unexpected symbol in unit macro.")

"""
```
basefactor(x::Unit)
```

Specifies conversion factors to reference units.
It returns a tuple. The first value is any irrational part of the conversion,
and the second value is a rational component. This segregation permits exact
conversions within unit systems that have no rational conversion to the
reference units.
"""
function basefactor end

"""
```
dimension(x::Unit)
```

Returns a [`Unitful.Dimensions`](@ref) object describing the given unit `x`.
"""
function dimension end


text: @deprecate showlattice viz


text:
@testset "Search: types, methods" begin
    INDEX_TYPE = UInt32
    for coarse_quantizer in [:naive, :hnsw]
        ivfadc = build_index_random_data(;
                    coarse_quantizer=coarse_quantizer,
                    index_type=INDEX_TYPE)

        # Single query
        K = 3  # number of neighbors
        query = rand(nrows)
        idxs, dists = knn_search(ivfadc, query, K, w=2)
        @test idxs isa Vector{INDEX_TYPE}
        @test dists isa Vector{eltype(query)}
        @test_throws AssertionError knn_search(ivfadc, query, 0)
        @test_throws AssertionError knn_search(ivfadc, query, 1, w=0)

        # Multiple queries
        queries = [rand(nrows) for _ in 1:10]
        idxs, dists = knn_search(ivfadc, queries, K, w=2)
        @test idxs isa Vector{Vector{INDEX_TYPE}}
        @test dists isa Vector{Vector{eltype(query)}}
    end
end

@testset "Search: results" begin
    data = [0    0     0     1    1  1    1    1  20  20    20    20    20;
            0.1  0.11  0.12  8  10  15  14  16  5  5.1  5.2  5.4  5.5]
    for coarse_quantizer in [:naive, :hnsw]
        ivfadc = IVFADCIndex(data, kc=3, k=8, m=2,
                             coarse_quantizer=coarse_quantizer)
        points = [[1.0,10.0], [0.0, 0.0], [20, 5.0]]

        neighbors_w1 = [[5, 4, 7, 6, 8],
                        [1, 2, 3],
                        [9, 10, 11, 12, 13]]
        for (point, result) in zip(points, neighbors_w1)
            neighbors = Int.(knn_search(ivfadc, point, 5, w=1)[1]) .+1
            @test isempty(setdiff(neighbors, result))
        end
        neighbors_w2 = [[5, 4, 7, 6, 8],
                        [1, 2, 3, 4, 5],
                        [9, 10, 11, 12, 13]]
        for (point, result) in zip(points, neighbors_w2)
            neighbors = Int.(knn_search(ivfadc, point, 5, w=2)[1]) .+1
            @test isempty(setdiff(neighbors, result))
        end
    end
end


text:
#__precompile__(true)

module ImagineAnalyses

using AxisArrays, Interpolations, Unitful, IntervalSets
using Reexport
@reexport using CachedCalls, ImagineInterface, UnitAliases

import UnitAliases.HasTimeUnits
import CachedCalls.CachedCall
import Unitful: μm, s, Hz
using DSP
using Statistics

export dumb_bench,
        get_cycles,
        largest_cycle_diff,
        largest_cycle_std,
        mean_cycle,

        find_circular,
        get_circular,
        mon_delay,
        pulse_timings,
        flash_cam_cycs,
        max_exp,
        framerate_efficiency

#This will be a bit of a hodgepodge of analyses one can do on ImagineSignals and acquired images.  May be worth organizing differently if it gets too big.

#this may do something like measure the lag between sig1 and sig2
dumb_bench(sig1::ImagineSignal, sig2::ImagineSignal) = 0.0

#Define new constructors for CachedCalls on a per-function basis when the return type is known (warning: risky when function isn't type stable)
CachedCall(f::typeof(dumb_bench), args...) = CachedCall("dumb benchmark", f, (args...,), Float64)

include("consistency.jl")
include("exposure_placement.jl")

#To implement:

#lag(sig1::ImagineSignal, sig2::ImagineSignal)
    #lag(sig1::AxisArray, sig2::AxisArray)
#offset(sig1::ImagineSignal, sig2::ImagineSignal)
#gain(sig_out::ImagineSignal, sig_in::ImagineSignal)
#corr(sig1::ImagineSignal, sig2::ImagineSignal)
#pulse_count(sig::ImagineSignal)
#range(sig::ImagineSignal)
#signal_consistency(sig::ImagineSignal, cycle_dur::HasTimeUnits) #consistency of a periodic 1d timeseries. may change the interface

#The below use images
#frame_count(camsig::ImagineSignal, img::AbstractImage{T,N}) where {T,N}
#frame_order(camsig::ImagineSignal, img::AbstractImage{T,N}) where {T,N}  #returns an integer sequence of frame numbers
#bidi_image_consistency(img::AbstractImage{T,4}) where {T} #returns each pair of ssds in a bidi sequence
#image_consistency(img::AbstractImage{T,N}, imgtemplate::AbstractImage{T,N-1}) where {T, N} #returns ssd between imgtemplate and each image in the img sequence

end # module


text:
using DelimitedFiles

input = readdlm(ARGS[1], ',', String, '\n');
parseDirection(element::String) = (h=element[1], δ=parse(Int, element[2:end]))
wire1 = input[1, :];
wire2 = input[2, :];
if (wire2[end] == "")
    pop!(wire2)
end
wire1 = parseDirection.(wire1);
wire2 = parseDirection.(wire2);

function explore(turns)
    position = 0
    path = Dict()
    inc = Dict('U' => im, 'D' => -im, 'R' => 1, 'L' => -1)
    steps = 0
    for t in turns
        for d in 1:t.δ
            position = position + inc[t.h] # Not my idea but it is clever. Saves us from switching.
            steps = steps + 1
            if !(position in keys(path))
                path[position] = steps
            end
        end
    end
    return path
end

path1 = explore(wire1)
path2 = explore(wire2)
crosses = collect(intersect(keys(path1), keys(path2)))

manhattan(p) = abs(real(p)) + abs(imag(p))
distances = map(manhattan, crosses)
println("Part 1: $(findmin(distances)[1])")

steps = map(point -> path1[point] + path2[point], crosses)
println("Part 2: $(findmin(steps)[1])")

# Whew, this one was pretty difficult! I had wanted to use tricky linear algebra.
# Glad I didn't -- the procedural method would have been the only feasible way to
# get the steps in part 2.


text:
@testset "Sets" begin
    # No argument.
    @testset "$(S)" for S in [
        CP.AbsoluteValue,
    ]
        @test isbitstype(S)

        @test S() == S()

        s = S()
        @test typeof(copy(s)) <: S
        @test copy(s) === s # Same object.

        @test MOI.dimension(S()) == 2
    end

    # Just a dimension.
    @testset "$(S)" for S in [
        CP.AllEqual,
        CP.AllDifferent,
        CP.SymmetricAllDifferent,
        CP.Membership,
        CP.ElementVariableArray,
        CP.CountDistinct,
        CP.Inverse,
        CP.Contiguity,
        CP.Circuit,
        CP.CircuitPath,
        CP.LexicographicallyLessThan,
        CP.LexicographicallyGreaterThan,
        CP.Sort,
        CP.SortPermutation,
        CP.MinimumAmong,
        CP.MaximumAmong,
        CP.ArgumentMinimumAmong,
        CP.ArgumentMaximumAmong,
        CP.Increasing,
        CP.Decreasing,
    ]
        @test isbitstype(S)

        @test S(2) == S(2)
        @test S(2) != S(3)
        @test S(3) != S(2)

        s = S(2)
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        if S in [
            CP.AllEqual,
            CP.AllDifferent,
            CP.SymmetricAllDifferent,
            CP.Membership,
            CP.Contiguity,
            CP.Circuit,
            CP.Increasing,
            CP.Decreasing,
        ]
            @test MOI.dimension(S(2)) == 2
            @test MOI.dimension(S(3)) == 3
        elseif S in [
            CP.CountDistinct,
            CP.MinimumAmong,
            CP.MaximumAmong,
            CP.ArgumentMinimumAmong,
            CP.ArgumentMaximumAmong,
        ]
            @test MOI.dimension(S(2)) == 2 + 1
            @test MOI.dimension(S(3)) == 3 + 1
        elseif S == CP.ElementVariableArray
            @test MOI.dimension(S(2)) == 2 + 2
            @test MOI.dimension(S(3)) == 3 + 2
        elseif S in [
            CP.Inverse,
            CP.CircuitPath,
            CP.LexicographicallyLessThan,
            CP.LexicographicallyGreaterThan,
            CP.Sort,
        ]
            @test MOI.dimension(S(2)) == 2 * 2
            @test MOI.dimension(S(3)) == 3 * 2
        elseif S == CP.SortPermutation
            @test MOI.dimension(S(2)) == 2 * 3
            @test MOI.dimension(S(3)) == 3 * 3
        else
            error("$(S) not implemented")
        end
    end

    # Two arguments: first a dimension, then a templated constant.
    @testset "$(S)" for S in [
        CP.MinimumDistance,
        CP.MaximumDistance,
    ]
        @test isbitstype(S{Int})

        @test S(2, 0) == S(2, 0)
        @test S(2, 0) != S(2, 1)
        @test S(2, 0) != S(3, 0)
        @test S(3, 0) != S(2, 0)

        s = S(2, 0)
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        if S in [
            CP.MaximumDistance,
            CP.MinimumDistance,
        ]
            @test MOI.dimension(S(2, 0)) == 2
            @test MOI.dimension(S(3, 4)) == 3
        elseif S == CP.Count
            @test MOI.dimension(S(2, 0)) == 2 + 1
            @test MOI.dimension(S(3, 4)) == 3 + 1
        else
            error("$(S) not implemented")
        end
    end

    # Two arguments: two dimensions.
    @testset "$(S)" for S in [
        CP.LexicographicallyLessThan,
        CP.LexicographicallyGreaterThan,
        CP.DoublyLexicographicallyLessThan,
        CP.DoublyLexicographicallyGreaterThan,
    ]
        @test isbitstype(S)

        @test S(2, 2) == S(2, 2)
        @test S(2, 2) != S(2, 1)
        @test S(2, 2) != S(3, 2)
        @test S(3, 2) != S(2, 2)

        s = S(2, 2)
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        if S == CP.LexicographicallyLessThan || S == CP.LexicographicallyGreaterThan || S == CP.DoublyLexicographicallyLessThan || S == CP.DoublyLexicographicallyGreaterThan
            @test MOI.dimension(S(2, 2)) == 2 * 2
            @test MOI.dimension(S(3, 4)) == 3 * 4
        else
            error("$(S) not implemented")
        end
    end

    @testset "AllDifferentExceptConstants" begin
        # Convenience constructor for one value.
        @test CP.AllDifferentExceptConstant(2, 0) == CP.AllDifferentExceptConstant(2, 0)
        @test CP.AllDifferentExceptConstant(2, 0) != CP.AllDifferentExceptConstant(2, 1)
        @test CP.AllDifferentExceptConstant(2, 0) != CP.AllDifferentExceptConstant(3, 0)
        @test CP.AllDifferentExceptConstant(3, 0) != CP.AllDifferentExceptConstant(2, 0)

        s = CP.AllDifferentExceptConstant(2, 0)
        @test typeof(copy(s)) <: CP.AllDifferentExceptConstants
        @test copy(s) == s

        @test MOI.dimension(CP.AllDifferentExceptConstant(2, 0)) == 2
        @test MOI.dimension(CP.AllDifferentExceptConstant(3, 4)) == 3

        # Usual constructor.
        @test CP.AllDifferentExceptConstants(2, Set([0, 1])) == CP.AllDifferentExceptConstants(2, Set([0, 1]))
        @test CP.AllDifferentExceptConstants(2, Set([0, 1])) != CP.AllDifferentExceptConstants(2, Set([2, 3]))
        @test CP.AllDifferentExceptConstants(2, Set([0, 1])) != CP.AllDifferentExceptConstants(3, Set([0, 1]))
        @test CP.AllDifferentExceptConstants(3, Set([0, 1])) != CP.AllDifferentExceptConstants(2, Set([0, 1]))

        s = CP.AllDifferentExceptConstants(2, Set([0, 1]))
        @test typeof(copy(s)) <: CP.AllDifferentExceptConstants
        @test copy(s) == s

        @test MOI.dimension(CP.AllDifferentExceptConstants(2, Set([0, 1]))) == 2
        @test MOI.dimension(CP.AllDifferentExceptConstants(3, Set([2, 3]))) == 3
    end

    @testset "Count{…}" begin
        @test isbitstype(CP.Count{MOI.EqualTo{Int}})

        # Default constructor: MOI.EqualTo.
        @test CP.Count(2, 0) == CP.Count(2, 0)
        @test CP.Count(2, 0) != CP.Count(2, 1)
        @test CP.Count(2, 0) != CP.Count(3, 0)
        @test CP.Count(3, 0) != CP.Count(2, 0)

        s = CP.Count(2, 0)
        @test typeof(copy(s)) <: CP.Count
        @test copy(s) == s

        @test MOI.dimension(CP.Count(2, 0)) == 2 + 1
        @test MOI.dimension(CP.Count(3, 4)) == 3 + 1

        # Directly give a MOI.EqualTo object.
        @test CP.Count(2, MOI.EqualTo(0)) == CP.Count(2, MOI.EqualTo(0))
        @test CP.Count(2, MOI.EqualTo(0)) != CP.Count(2, MOI.EqualTo(1))
        @test CP.Count(2, MOI.EqualTo(0)) != CP.Count(3, MOI.EqualTo(0))
        @test CP.Count(3, MOI.EqualTo(0)) != CP.Count(2, MOI.EqualTo(0))

        s = CP.Count(2, MOI.EqualTo(0))
        @test typeof(copy(s)) <: CP.Count
        @test copy(s) == s

        @test MOI.dimension(CP.Count(2, MOI.EqualTo(0))) == 2 + 1
        @test MOI.dimension(CP.Count(3, MOI.EqualTo(4))) == 3 + 1

        # Other sets.
        @test CP.Count(2, CP.DifferentFrom(0)) == CP.Count(2, CP.DifferentFrom(0))
        @test CP.Count(2, CP.DifferentFrom(0)) != CP.Count(2, CP.DifferentFrom(1))
        @test CP.Count(2, CP.DifferentFrom(0)) != CP.Count(3, CP.DifferentFrom(0))
        @test CP.Count(3, CP.DifferentFrom(0)) != CP.Count(2, CP.DifferentFrom(0))

        s = CP.Count(2, CP.DifferentFrom(0))
        @test typeof(copy(s)) <: CP.Count
        @test copy(s) == s

        @test MOI.dimension(CP.Count(2, CP.DifferentFrom(0))) == 2 + 1
        @test MOI.dimension(CP.Count(3, CP.DifferentFrom(4))) == 3 + 1
    end

    @testset "CountCompare" begin
        @test isbitstype(CP.CountCompare)

        @test CP.CountCompare(2) == CP.CountCompare(2)
        @test CP.CountCompare(2) != CP.CountCompare(3)
        @test CP.CountCompare(3) != CP.CountCompare(2)

        s = CP.CountCompare(2)
        @test typeof(copy(s)) <: CP.CountCompare
        @test copy(s) == s

        @test MOI.dimension(CP.CountCompare(2)) == 2 * 2 + 1
        @test MOI.dimension(CP.CountCompare(3)) == 3 * 2 + 1
    end

    @testset "SlidingSum" begin
        @test CP.SlidingSum(1, 2, 3, 4) == CP.SlidingSum(1, 2, 3, 4)
        @test CP.SlidingSum(1, 2, 3, 4) != CP.SlidingSum(0, 2, 3, 4)
        @test CP.SlidingSum(1, 2, 3, 4) != CP.SlidingSum(1, 4, 2, 4)
        @test CP.SlidingSum(1, 2, 3, 4) != CP.SlidingSum(1, 2, 3, 5)

        s = CP.SlidingSum(1, 2, 3, 4)
        @test typeof(copy(s)) <: CP.SlidingSum
        @test copy(s) == s

        @test MOI.dimension(CP.SlidingSum(1, 2, 3, 4)) == 4
        @test MOI.dimension(CP.SlidingSum(1, 2, 3, 40)) == 40
    end

    @testset "Element" begin
        set1 = [1, 2, 3]
        set2 = [1, 2, 3, 4]
        @test CP.Element(set1) == CP.Element(set1)
        @test CP.Element(set1) != CP.Element(set2)
        @test CP.Element(set2) != CP.Element(set1)

        s = CP.Element(set1)
        @test typeof(copy(s)) <: CP.Element
        @test copy(s) == s

        @test MOI.dimension(CP.Element(set1)) == 2
        @test MOI.dimension(CP.Element(set2)) == 2
    end

    @testset "$(S)" for S in [CP.Domain, CP.AntiDomain]
        set1 = Set([1, 2, 3])
        set2 = Set([1, 2, 3, 4])
        @test S(set1) == S(set1)
        @test S(set1) != S(set2)
        @test S(set2) != S(set1)
        @test S(Set([0, 1])) == S(Set([0, 1]))
        @test S(Set([0, 1])) != S(Set([0, 2]))

        s = S(set1)
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(S(set1)) == 1
        @test MOI.dimension(S(set2)) == 1
    end

    @testset "$(S)" for S in [CP.VectorDomain, CP.VectorAntiDomain]
        set1 = Set([[1, 2, 3], [2, 3, 4]])
        set2 = Set([[1, 2, 3], [3, 4, 5]])
        @test S(3, set1) == S(3, set1)
        @test S(3, set1) != S(3, set2)
        @test S(3, set2) != S(3, set1)

        s = S(3, set1)
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(S(3, set1)) == 3
        @test MOI.dimension(S(3, set2)) == 3
    end

    @testset "Strictly{$(Ssub)}" for Ssub in [MOI.LessThan, MOI.GreaterThan]
        @test CP.Strictly(Ssub(1)) == CP.Strictly(Ssub(1))
        @test CP.Strictly(Ssub(1)) != CP.Strictly(Ssub(2))
        @test CP.Strictly(Ssub(2)) == CP.Strictly(Ssub(2))
        @test CP.Strictly(Ssub(2)) != CP.Strictly(Ssub(1))

        s = CP.Strictly(Ssub(1))
        @test typeof(copy(s)) <: CP.Strictly{Ssub{Int}}
        @test copy(s) == s

        @test MOI.dimension(CP.Strictly(Ssub(3))) == 1
        @test MOI.dimension(CP.Strictly(Ssub(1))) == 1
        @test MOI.dimension(CP.Strictly(Ssub(3))) == MOI.dimension(Ssub(3))

        @test MOI.constant(CP.Strictly(Ssub(3))) == 3
        @test MOIU.shift_constant(CP.Strictly(Ssub(3)), 1) ==
              CP.Strictly(Ssub(4))
    end

    @testset "Strictly{$(Ssub), $(T)}, one argument" for Ssub in [
        CP.LexicographicallyLessThan,
        CP.LexicographicallyGreaterThan,
        CP.Increasing,
        CP.Decreasing,
    ], T in [Int, Float64]
        @test CP.Strictly{Ssub, T}(Ssub(1)) == CP.Strictly{Ssub, T}(Ssub(1))
        @test CP.Strictly{Ssub, T}(Ssub(1)) != CP.Strictly{Ssub, T}(Ssub(2))
        @test CP.Strictly{Ssub, T}(Ssub(2)) == CP.Strictly{Ssub, T}(Ssub(2))
        @test CP.Strictly{Ssub, T}(Ssub(2)) != CP.Strictly{Ssub, T}(Ssub(1))

        s = CP.Strictly{Ssub, T}(Ssub(1))
        @test typeof(copy(s)) <: CP.Strictly{Ssub, T}
        @test copy(s) == s

        if Ssub in
           [CP.LexicographicallyLessThan, CP.LexicographicallyGreaterThan]
            @test MOI.dimension(CP.Strictly(Ssub(3))) == 3 * 2
            @test MOI.dimension(CP.Strictly(Ssub(1))) == 1 * 2
        elseif Ssub in [CP.Increasing, CP.Decreasing]
            @test MOI.dimension(CP.Strictly(Ssub(3))) == 3
            @test MOI.dimension(CP.Strictly(Ssub(1))) == 1
        else
            error("$(Ssub) not implemented")
        end

        if Ssub in [MOI.LessThan, MOI.GreaterThan]
            @test MOI.constant(CP.Strictly(Ssub(3))) == 3
            @test MOIU.shift_constant(CP.Strictly(Ssub(3)), 1) ==
                  CP.Strictly(Ssub(4))
        end
    end

    @testset "Strictly{$(Ssub), $(T)}, two arguments" for Ssub in [
        CP.LexicographicallyLessThan,
        CP.LexicographicallyGreaterThan,
        CP.DoublyLexicographicallyLessThan,
        CP.DoublyLexicographicallyGreaterThan,
    ], T in [Int, Float64]
        @test CP.Strictly{Ssub, T}(Ssub(1, 4)) == CP.Strictly{Ssub, T}(Ssub(1, 4))
        @test CP.Strictly{Ssub, T}(Ssub(1, 4)) != CP.Strictly{Ssub, T}(Ssub(2, 4))
        @test CP.Strictly{Ssub, T}(Ssub(2, 4)) == CP.Strictly{Ssub, T}(Ssub(2, 4))
        @test CP.Strictly{Ssub, T}(Ssub(2, 4)) != CP.Strictly{Ssub, T}(Ssub(1, 4))

        s = CP.Strictly{Ssub, T}(Ssub(1, 4))
        @test typeof(copy(s)) <: CP.Strictly{Ssub, T}
        @test copy(s) == s

        if Ssub in
           [CP.LexicographicallyLessThan, CP.LexicographicallyGreaterThan,
           CP.DoublyLexicographicallyLessThan,
           CP.DoublyLexicographicallyGreaterThan,]
            @test MOI.dimension(CP.Strictly(Ssub(3, 4))) == 3 * 4
            @test MOI.dimension(CP.Strictly(Ssub(1, 4))) == 1 * 4
        else
            error("$(Ssub) not implemented")
        end
    end

    @testset "DifferentFrom" begin
        @test isbitstype(CP.DifferentFrom{Float64})

        @test CP.DifferentFrom(1) == CP.DifferentFrom(1)
        @test CP.DifferentFrom(1) != CP.DifferentFrom(2)
        @test CP.DifferentFrom(2) == CP.DifferentFrom(2)
        @test CP.DifferentFrom(2) != CP.DifferentFrom(1)

        s = CP.DifferentFrom(1)
        @test typeof(copy(s)) <: CP.DifferentFrom
        @test copy(s) == s

        @test MOI.constant(CP.DifferentFrom(3)) == 3
        @test MOI.dimension(CP.DifferentFrom(3)) == 1
        @test MOIU.shift_constant(CP.DifferentFrom(3), 1) == CP.DifferentFrom(4)
    end

    @testset "$(S)" for S in [CP.Reification, CP.Negation]
        @test S(MOI.EqualTo(0.0)) == S(MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(0.0)) != S(MOI.EqualTo(1.0))
        @test S(MOI.EqualTo(1.0)) != S(MOI.EqualTo(0.0))
        @test S(MOI.GreaterThan(0.0)) != S(MOI.EqualTo(1.0))
        @test S(MOI.EqualTo(1.0)) != S(MOI.GreaterThan(0.0))

        s = S(MOI.EqualTo(0.0))
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        if S == CP.Reification
            @test MOI.dimension(S(MOI.EqualTo(0.0))) ==
                  1 + MOI.dimension(MOI.EqualTo(0.0))
            @test MOI.dimension(S(MOI.GreaterThan(0.0))) ==
                  1 + MOI.dimension(MOI.GreaterThan(0.0))
        elseif S == CP.Negation
            @test MOI.dimension(S(MOI.EqualTo(0.0))) ==
                  MOI.dimension(MOI.EqualTo(0.0))
            @test MOI.dimension(S(MOI.GreaterThan(0.0))) ==
                  MOI.dimension(MOI.GreaterThan(0.0))
        else
            error("$(S) not implemented")
        end
    end

    @testset "$(S)" for S in [CP.Equivalence, CP.EquivalenceNot, CP.Implication]
        @test S(MOI.EqualTo(0.0), MOI.EqualTo(0.0)) ==
              S(MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(1.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(1.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(1.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(1.0), MOI.EqualTo(0.0)) !=
              S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0))

        s = S(MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(S(MOI.EqualTo(0.0), MOI.EqualTo(0.0))) ==
              2 * MOI.dimension(MOI.EqualTo(0.0))
        @test MOI.dimension(S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0))) ==
              2 * MOI.dimension(MOI.GreaterThan(0.0))
    end

    @testset "$(S)" for S in [CP.IfThenElse]
        @test S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) ==
              S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test S(MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
              S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))

        s = S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(
            S(MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)),
        ) == 3 * MOI.dimension(MOI.EqualTo(0.0))
        @test MOI.dimension(
            S(MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)),
        ) == 3 * MOI.dimension(MOI.GreaterThan(0.0))
    end

    @testset "$(S)" for S in [CP.Conjunction, CP.Disjunction, CP.ExclusiveDisjunction]
        # Ensure that tuples can be compared.
        @test (MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) ==
            (MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test (MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
                (MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test (MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
                (MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test (MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
                (MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test (MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)) !=
                (MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))
        @test (MOI.EqualTo(1), MOI.EqualTo(0), MOI.EqualTo(0)) !=
                (MOI.GreaterThan(0), MOI.EqualTo(0), MOI.EqualTo(0))
        @test (CP.Domain(Set([0, 1])),) == (CP.Domain(Set([0, 1])),)

        # Then, wrap the tuples in conjunctions and disjunctions.
        @test S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))) ==
              S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))) !=
              S((MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test S((MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))) !=
              S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test S((MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))) !=
              S((MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test S((MOI.EqualTo(1.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))) !=
              S((MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test S((MOI.EqualTo(1), MOI.EqualTo(0), MOI.EqualTo(0))) !=
              S((MOI.GreaterThan(0), MOI.EqualTo(0), MOI.EqualTo(0)))
        @test S((CP.Domain(Set([0, 1])), MOI.EqualTo(0))) ==
              S((CP.Domain(Set([0, 1])), MOI.EqualTo(0)))
        @test S((CP.Domain(Set([0, 1])),)) == S((CP.Domain(Set([0, 1])),))

        s = S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0)))
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(
            S((MOI.EqualTo(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))),
        ) == 3 * MOI.dimension(MOI.EqualTo(0.0))
        @test MOI.dimension(
            S((MOI.GreaterThan(0.0), MOI.EqualTo(0.0), MOI.EqualTo(0.0))),
        ) == 3 * MOI.dimension(MOI.GreaterThan(0.0))
    end

    @testset "$(S)" for S in [CP.True, CP.False]
        @test isbitstype(S)

        @test S() == S()

        s = S()
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        @test MOI.dimension(S()) == 0
    end

    @testset "ValuePrecedence" begin
        @test CP.ValuePrecedence(1, 4, 20) == CP.ValuePrecedence(1, 4, 20)
        @test CP.ValuePrecedence(1, 4, 20) != CP.ValuePrecedence(1, 4, 25)
        @test CP.ValuePrecedence(1, 4, 20) != CP.ValuePrecedence(1, 5, 20)
        @test CP.ValuePrecedence(1, 4, 20) != CP.ValuePrecedence(2, 4, 20)

        s = CP.ValuePrecedence(1, 4, 20)
        @test typeof(copy(s)) <: CP.ValuePrecedence
        @test copy(s) == s

        @test MOI.dimension(CP.ValuePrecedence(1, 4, 20)) == 20
        @test MOI.dimension(CP.ValuePrecedence(1, 5, 20)) == 20
        @test MOI.dimension(CP.ValuePrecedence(2, 4, 20)) == 20
    end

    @testset "$(S)" for S in [CP.WeightedCircuit, CP.WeightedCircuitPath]
        @test S(3, [1 2 3; 4 5 6; 7 8 9]) == S(3, [1 2 3; 4 5 6; 7 8 9])
        @test S(3, [1 2 3; 4 5 6; 7 8 9]) != S(2, [1 2; 3 4])
        @test S(2, [1 2; 3 4]) != S(3, [1 2 3; 4 5 6; 7 8 9])

        s = S(3, [1 2 3; 4 5 6; 7 8 9])
        @test typeof(copy(s)) <: S
        @test copy(s) == s

        if S == CP.WeightedCircuit
            @test MOI.dimension(S(3, [1 2 3; 4 5 6; 7 8 9])) == 3 + 1
        elseif S == CP.WeightedCircuitPath
            @test MOI.dimension(S(3, [1 2 3; 4 5 6; 7 8 9])) == 3 * 2 + 1
        else
            error("$(S) not implemented")
        end
    end

    @testset "GlobalCardinality family" begin
        @testset "Constructor" begin
            for CVCT in [CP.OPEN_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES]
                @test_throws ErrorException CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CVCT, Int}(4, Int[], 5)
                @test_throws ErrorException CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CVCT, Int}(4, [4, 5], -1)
            end
        end

        @testset "GlobalCardinality{FIXED_COUNTED_VALUES, $(CVCT), Int}" for CVCT in [CP.OPEN_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES]
            @test CP.GlobalCardinality{CVCT}(2, [2, 4]) == CP.GlobalCardinality{CVCT}(2, [2, 4])
            @test CP.GlobalCardinality{CVCT}(2, [2, 4]) != CP.GlobalCardinality{CVCT}(3, [2, 4])
            @test CP.GlobalCardinality{CVCT}(3, [2, 4]) != CP.GlobalCardinality{CVCT}(2, [2, 4])
            @test CP.GlobalCardinality{CVCT}(2, [2, 4]) != CP.GlobalCardinality{CVCT}(2, [3, 5])

            s = CP.GlobalCardinality{CVCT}(2, [2, 4])
            @test typeof(copy(s)) <: CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CVCT, Int}
            @test copy(s) == s

            @test MOI.dimension(CP.GlobalCardinality{CVCT}(2, [2, 4])) == 2 + 2
            @test MOI.dimension(CP.GlobalCardinality{CVCT}(3, [2, 4, 6, 8])) == 3 + 4

            if CVCT == CP.OPEN_COUNTED_VALUES
                @test CP.GlobalCardinality{CVCT}(2, [2, 4]) == CP.GlobalCardinality(2, [2, 4])
            end
        end

        @testset "GlobalCardinality{VARIABLE_COUNTED_VALUES, $(CVCT), Int}" for CVCT in [CP.OPEN_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES]
            @test CP.GlobalCardinality{CVCT, Int}(2, 2) == CP.GlobalCardinality{CVCT, Int}(2, 2)
            @test CP.GlobalCardinality{CVCT, Int}(2, 2) != CP.GlobalCardinality{CVCT, Int}(3, 3)
            @test CP.GlobalCardinality{CVCT, Int}(3, 2) != CP.GlobalCardinality{CVCT, Int}(2, 2)
            @test CP.GlobalCardinality{CVCT, Int}(2, 2) != CP.GlobalCardinality{CVCT, Int}(2, 3)

            s = CP.GlobalCardinality{CVCT, Int}(2, 2)
            @test typeof(copy(s)) <: CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CVCT, Int}
            @test copy(s) == s

            @test MOI.dimension(CP.GlobalCardinality{CVCT, Int}(2, 2)) == 2 + 2 * 2
            @test MOI.dimension(CP.GlobalCardinality{CVCT, Int}(3, 4)) == 3 + 2 * 4

            if CVCT == CP.OPEN_COUNTED_VALUES
                @test CP.GlobalCardinality{CVCT, Int}(2, 3) == CP.GlobalCardinality{Int}(2, 3)
            end
        end
    end

    @testset "CumulativeResource family" begin
        @testset "CumulativeResource{$(CRDT)}" for CRDT in [
            CP.NO_DEADLINE_CUMULATIVE_RESOURCE,
            CP.VARIABLE_DEADLINE_CUMULATIVE_RESOURCE,
        ]
            @test CP.CumulativeResource{CRDT}(2) == CP.CumulativeResource{CRDT}(2)
            @test CP.CumulativeResource{CRDT}(2) != CP.CumulativeResource{CRDT}(3)
            @test CP.CumulativeResource{CRDT}(3) != CP.CumulativeResource{CRDT}(2)

            s = CP.CumulativeResource{CRDT}(2)
            @test typeof(copy(s)) <: CP.CumulativeResource{CRDT}
            @test copy(s) == s

            if CRDT == CP.NO_DEADLINE_CUMULATIVE_RESOURCE
                @test MOI.dimension(CP.CumulativeResource{CRDT}(2)) == 2 * 3 + 1
                @test MOI.dimension(CP.CumulativeResource{CRDT}(3)) == 3 * 3 + 1
            elseif CRDT == CP.VARIABLE_DEADLINE_CUMULATIVE_RESOURCE
                @test MOI.dimension(CP.CumulativeResource{CRDT}(2)) == 2 * 4 + 1
                @test MOI.dimension(CP.CumulativeResource{CRDT}(3)) == 3 * 4 + 1
            else
                @assert false
            end
        end

        @test CP.CumulativeResource(2) == CP.CumulativeResource(2)
        @test CP.CumulativeResource(2) != CP.CumulativeResource(3)
        @test CP.CumulativeResource(3) != CP.CumulativeResource(2)

        s = CP.CumulativeResource(2)
        @test typeof(copy(s)) <: CP.CumulativeResource{CP.NO_DEADLINE_CUMULATIVE_RESOURCE}
        @test copy(s) == s

        @test MOI.dimension(CP.CumulativeResource(2)) == 2 * 3 + 1
        @test MOI.dimension(CP.CumulativeResource(3)) == 3 * 3 + 1
    end

    @testset "BinPacking family" begin
        @testset "BinPacking{NO_CAPACITY_BINPACKING}" begin
            @test_throws AssertionError CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2, 3])

            @test CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2]) == CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(2, 2, [1, 2]) != CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 3, [1, 2, 3]) != CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2])

            s = CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test typeof(copy(s)) <: CP.BinPacking{CP.NO_CAPACITY_BINPACKING}
            @test copy(s) == s

            @test MOI.dimension(CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [1, 2])) == 3

            @test_throws AssertionError CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [-1, 2])
            @test_throws AssertionError CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(0, 2, [1, 2])
            @test_throws AssertionError CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 0, [1, 2])
        end

        @testset "BinPacking{FIXED_CAPACITY_BINPACKING}" begin
            @test_throws AssertionError CP.BinPacking(
                1,
                2,
                [1, 2, 3],
                [4],
            )
            @test_throws AssertionError CP.BinPacking(
                1,
                2,
                [1, 2],
                [3, 4],
            )
            @test_throws AssertionError CP.BinPacking(
                1,
                2,
                [1, 2, 3],
                [4, 5],
            )

            @test CP.BinPacking(1, 2, [1, 2], [4]) ==
                CP.BinPacking(1, 2, [1, 2], [4])
            @test CP.BinPacking(2, 2, [1, 2], [4, 5]) !=
                CP.BinPacking(1, 2, [1, 2], [4])
            @test CP.BinPacking(1, 3, [1, 2, 3], [4]) !=
                CP.BinPacking(1, 2, [1, 2], [4])

            s = CP.BinPacking(1, 2, [1, 2], [4])
            @test typeof(copy(s)) <: CP.BinPacking{CP.FIXED_CAPACITY_BINPACKING, Int}
            @test copy(s) == s

            @test MOI.dimension(CP.BinPacking(1, 2, [1, 2], [4])) == 3

            @test_throws AssertionError CP.BinPacking(1, 2, [-1, 2], [4])
            @test_throws AssertionError CP.BinPacking(1, 2, [1, 2], [-4])
            @test_throws AssertionError CP.BinPacking(0, 2, [1, 2], [4])
            @test_throws AssertionError CP.BinPacking(1, 0, [1, 2], [4])
        end

        @testset "BinPacking{VARIABLE_CAPACITY_BINPACKING}" begin
            @test_throws AssertionError CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(
                1,
                2,
                [1, 2, 3],
            )

            @test CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2]) ==
                CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(2, 2, [1, 2]) !=
                CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 3, [1, 2, 3]) !=
                CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2])

            s = CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2])
            @test typeof(copy(s)) <: CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}
            @test copy(s) == s

            @test MOI.dimension(CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [1, 2])) == 4

            @test_throws AssertionError CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 2, [-1, 2])
            @test_throws AssertionError CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(0, 2, [1, 2])
            @test_throws AssertionError CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING}(1, 0, [1, 2])
        end
    end

    @testset "Knapsack family" begin
        # @testset "Easy constructor" begin
        #     @test_throws UndefKeywordError CP.Knapsack(capacity=1)

        #     @test CP.Knapsack(weights=[1, 2]) == CP.Knapsack([1, 2])
        #     @test CP.Knapsack(weights=[1, 2]) == CP.Knapsack{VARIABLE_CAPACITY_KNAPSACK, UNVALUED_KNAPSACK, Int}([1, 2])
        # end

        @testset "Knapsack{FIXED_CAPACITY_KNAPSACK, UNVALUED_KNAPSACK}" begin
            @test CP.Knapsack([1, 2, 3], 3) == CP.Knapsack([1, 2, 3], 3)
            @test CP.Knapsack([1, 2, 3], 3) != CP.Knapsack([1, 2, 3], 4)
            @test CP.Knapsack([1, 2, 3], 4) != CP.Knapsack([1, 2, 3], 3)

            s = CP.Knapsack([1, 2, 3], 3)
            @test typeof(copy(s)) <: CP.Knapsack
            @test copy(s) == s

            @test MOI.dimension(CP.Knapsack([1, 2, 3], 3)) == 3

            @test_throws AssertionError CP.Knapsack([-1, 2, 3], 3)
            @test_throws AssertionError CP.Knapsack([1, 2, 3], -3)
        end

        @testset "Knapsack{VARIABLE_CAPACITY_KNAPSACK, UNVALUED_KNAPSACK}" begin
            @test CP.Knapsack([1, 2, 3]) == CP.Knapsack([1, 2, 3])

            s = CP.Knapsack([1, 2, 3])
            @test typeof(copy(s)) <: CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK}
            @test copy(s) == s

            @test MOI.dimension(CP.Knapsack([1, 2, 3])) == 3 + 1

            @test_throws AssertionError CP.Knapsack([-1, 2, 3])
        end

        @testset "Knapsack{FIXED_CAPACITY_KNAPSACK, VALUED_KNAPSACK}" begin
            @test CP.Knapsack([1, 2, 3], 3, [1, 2, 3]) == CP.Knapsack([1, 2, 3], 3, [1, 2, 3])
            @test CP.Knapsack([1, 2, 3], 3, [1, 2, 3]) != CP.Knapsack([1, 2, 3], 4, [1, 2, 3])
            @test CP.Knapsack([1, 2, 3], 4, [1, 2, 3]) != CP.Knapsack([1, 2, 3], 3, [1, 2, 3])

            s = CP.Knapsack([1, 2, 3], 3, [1, 2, 3])
            @test typeof(copy(s)) <: CP.Knapsack
            @test copy(s) == s

            @test MOI.dimension(CP.Knapsack([1, 2, 3], 3, [1, 2, 3])) == 3 + 1

            @test_throws AssertionError CP.Knapsack([-1, 2, 3], 3, [1, 2, 3])
            @test_throws AssertionError CP.Knapsack([1, 2, 3], 3, [-1, 2, 3])
            @test_throws AssertionError CP.Knapsack([1, 2, 3], -3, [1, 2, 3])
        end

        @testset "Knapsack{VARIABLE_CAPACITY_KNAPSACK, VALUED_KNAPSACK}" begin
            @test CP.Knapsack([1, 2, 3], [1, 2, 3]) == CP.Knapsack([1, 2, 3], [1, 2, 3])

            s = CP.Knapsack([1, 2, 3], [1, 2, 3])
            @test typeof(copy(s)) <: CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.VALUED_KNAPSACK}
            @test copy(s) == s

            @test MOI.dimension(CP.Knapsack([1, 2, 3], [1, 2, 3])) == 3 + 2

            @test_throws AssertionError CP.Knapsack([-1, 2, 3], [1, 2, 3])
            @test_throws AssertionError CP.Knapsack([1, 2, 3], [-1, 2, 3])
        end
    end

    @testset "NonOverlappingOrthotopes family" begin
        for NOOCT in [
            CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES,
            CP.CONDITIONAL_NONVERLAPPING_ORTHOTOPES,
        ]
            @test isbitstype(CP.NonOverlappingOrthotopes{NOOCT})

            @test CP.NonOverlappingOrthotopes{NOOCT}(2, 2) == CP.NonOverlappingOrthotopes{NOOCT}(2, 2)
            @test CP.NonOverlappingOrthotopes{NOOCT}(2, 2) != CP.NonOverlappingOrthotopes{NOOCT}(2, 1)
            @test CP.NonOverlappingOrthotopes{NOOCT}(2, 2) != CP.NonOverlappingOrthotopes{NOOCT}(3, 2)
            @test CP.NonOverlappingOrthotopes{NOOCT}(3, 2) != CP.NonOverlappingOrthotopes{NOOCT}(2, 2)

            s = CP.NonOverlappingOrthotopes{NOOCT}(2, 2)
            @test typeof(copy(s)) <: CP.NonOverlappingOrthotopes{NOOCT}
            @test copy(s) == s

            if NOOCT == CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES
                @test MOI.dimension(CP.NonOverlappingOrthotopes{NOOCT}(2, 2)) == 3 * 2 * 2
                @test MOI.dimension(CP.NonOverlappingOrthotopes{NOOCT}(3, 4)) == 3 * 3 * 4
            elseif NOOCT == CP.CONDITIONAL_NONVERLAPPING_ORTHOTOPES
                @test MOI.dimension(CP.NonOverlappingOrthotopes{NOOCT}(2, 2)) == 3 * 2 * 2 + 2
                @test MOI.dimension(CP.NonOverlappingOrthotopes{NOOCT}(3, 4)) == 3 * 3 * 4 + 3
            else
                error("$(NOOCT) not implemented")
            end

            if NOOCT == CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES
                @test CP.NonOverlappingOrthotopes{NOOCT}(2, 2) == CP.NonOverlappingOrthotopes(2, 2)
            end
        end
    end
end