text:
import StatsFuns: log2π

@node GammaShapeRate Stochastic [ out, (α, aliases = [ shape ]), (β, aliases = [ rate ]) ]

@average_energy GammaShapeRate (q_out::Any, q_α::PointMass, q_β::Any) = begin
    loggammamean(q_α) - mean(q_α) * logmean(q_β) - (mean(q_α) - 1.0) * logmean(q_out) + mean(q_β) * mean(q_out)
end

@average_energy GammaShapeRate (q_out::Any, q_α::GammaDistributionsFamily, q_β::Any) = begin
    loggammamean(q_α) - mean(q_α) * logmean(q_β) - (mean(q_α) - 1.0) * logmean(q_out) + mean(q_β) * mean(q_out)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# Defines the NonlinearExpression, BufferedNonlinearFunction used in
# constructing relaxations of nonlinear functions along with a number of
# helper functions including an Evaluator structure and: set_node_flag!
# set_node!, set_reference_point!, retrieve_node, prior_eval
# copy_subexpression_value!, eliminate_fixed_variables!
#############################################################################

"""
$(TYPEDEF)

Stores a general quadratic function with a buffer.
"""
mutable struct NonlinearExpression{V} <: AbstractEAGOConstraint

    "List of nodes in nonlinear expression"
    nd::Vector{NodeData}
    "Adjacency Matrix for the expression"
    adj::SparseMatrixCSC{Bool,Int64}
    const_values::Vector{Float64}

    setstorage::Vector{V}
    numberstorage::Vector{Float64}
    isnumber::Vector{Bool}
    value::V
    value_available::Bool

    tp1storage::Vector{Float64}
    tp2storage::Vector{Float64}
    tp3storage::Vector{Float64}
    tp4storage::Vector{Float64}
    tpdict::Dict{Int64,Tuple{Int64,Int64,Int64,Int64}}

    # sparsity of constraint + indices in node to reference
    grad_sparsity::Vector{Int64}       # indices of variables in the problem space (size = np)
    reverse_sparsity::Vector{Int64}

    # role in problem
    dependent_variable_count::Int64
    dependent_subexpression_count::Int64
    dependent_subexpressions::Vector{Int64}
    linearity::JuMP._Derivatives.Linearity

    # buffer for subgradients
    cv_grad_buffer::Vector{Float64}
    cc_grad_buffer::Vector{Float64}
end

"""
$(TYPEDEF)

Stores a general nonlinear function with a buffer represented by the sum of a tape
and a scalar affine function.
"""
mutable struct BufferedNonlinearFunction{V} <: AbstractEAGOConstraint
    expr::NonlinearExpression{V}
    saf::SAF
    lower_bound::Float64
    upper_bound::Float64
    last_relax_convex::Bool
    last_relax_concave::Bool
    last_past_reverse::Bool
    has_value::Bool
end

###
### Constructor definitions
###
function NonlinearExpression!(sub::JuMP._SubexpressionStorage, sub_sparsity_dict::Dict{Int64,Vector{Int64}},
                             subexpr_indx::Int64, subexpr_linearity::Vector{JuMP._Derivatives.Linearity},
                             tag::T) where T
    nd = copy(sub.nd)
    adj = copy(sub.adj)
    const_values = copy(sub.const_values)

    lenx = length(nd)
    numberstorage = zeros(lenx)
    isnumber = fill(false, lenx)

    # creates tiepoint storage and populates the sparsity of subexpression
    # not including dependent subexpressions
    tpdict = Dict{Int64,Tuple{Int64,Int64,Int64,Int64}}()
    tp1_count = 0
    tp2_count = 0
    for i = 1:lenx
        # tiepoint storage part
        node = @inbounds nd[i]
        op = node.index
        if double_tp(op)
            tp1_count += 1
            tp2_count += 1
            tpdict[i] = (tp1_count, tp1_count, tp2_count, tp2_count)
        elseif single_tp(op)
            tp1_count += 1
            tpdict[i] = (tp1_count, tp1_count, -1, -1)
        end
    end
    tp1storage = zeros(tp1_count)
    tp2storage = zeros(tp1_count)
    tp3storage = zeros(tp2_count)
    tp4storage = zeros(tp2_count)

    linearity = JuMP._Derivatives.classify_linearity(nd, adj, subexpr_linearity)

    # counts variables in subexpression and collects list of subexpressions
    dependent_subexpressions = Int64[]
    variable_dict = Dict{Int,Bool}()
    for (i,node) in enumerate(nd)
        if node.nodetype === JuMP._Derivatives.VARIABLE
            indx = node.index
            if !haskey(variable_dict, indx)
                variable_dict[indx] = true
            end
        end
        if node.nodetype === JuMP._Derivatives.VALUE
            indx = node.index
            numberstorage[i] = const_values[indx]
        end
        if node.nodetype === JuMP._Derivatives.SUBEXPRESSION
            indx = node.index
            push!(dependent_subexpressions, indx)
        end
    end
    grad_sparsity = collect(keys(variable_dict))
    unique!(dependent_subexpressions)
    sort!(dependent_subexpressions)

    # adds sparsity of dependent subexpressions to sparsity of current subexpression
    dependent_subexpression_count = length(dependent_subexpressions)
    for i = 1:dependent_subexpression_count
        dependent_indx = dependent_subexpressions[i]
        append!(grad_sparsity, sub_sparsity_dict[dependent_indx])
    end
    unique!(grad_sparsity)
    sort!(grad_sparsity)

    # sets subexpression sparsity dictionary
    sub_sparsity_dict[subexpr_indx] = copy(grad_sparsity)

    reverse_sparsity_length = grad_sparsity[end]
    reverse_sparsity = zeros(Int64, reverse_sparsity_length)
    current_grad_sparsity = grad_sparsity[1]
    current_grad_sparsity_count = 1
    for i = 1:reverse_sparsity_length
        if i == current_grad_sparsity
            reverse_sparsity[i] = current_grad_sparsity_count
            current_grad_sparsity_count += 1
            if current_grad_sparsity_count <= length(grad_sparsity)
                current_grad_sparsity = grad_sparsity[current_grad_sparsity_count]
            else
                break
            end
        end
    end

    dependent_variable_count = length(grad_sparsity)
    N = dependent_variable_count

    cv_grad_buffer = zeros(dependent_variable_count)
    cc_grad_buffer = zeros(dependent_variable_count)

    setstorage = fill(MC{N,T}(Interval(-Inf, Inf)), lenx)
    subexpression = NonlinearExpression{MC{N,T}}(nd, adj, const_values, setstorage, numberstorage,
                                                 isnumber, zero(MC{N,T}), false,
                                                 tp1storage, tp2storage,
                                                 tp3storage, tp4storage, tpdict, grad_sparsity,
                                                 reverse_sparsity, dependent_variable_count,
                                                 dependent_subexpression_count,
                                                 dependent_subexpressions, JuMP._Derivatives.CONSTANT,
                                                 cv_grad_buffer, cc_grad_buffer)
    return subexpression
end

function NonlinearExpression()
    return NonlinearExpression{MC{1,NS}}(JuMP.NodeData[], spzeros(Bool, 1), Float64[],
                                         MC{1,NS}[], Float64[], Bool[], zero(MC{1,NS}), false,
                                         Float64[], Float64[], Float64[], Float64[],
                                         Dict{Int64,Tuple{Int64,Int64,Int64,Int64}}(),
                                         Int64[], Int64[], 0, 0, Int64[], JuMP._Derivatives.CONSTANT,
                                         Float64[], Float64[])
end

function BufferedNonlinearFunction(func::JuMP._FunctionStorage, bnds::MOI.NLPBoundsPair,
                                   sub_sparsity_dict::Dict{Int64,Vector{Int64}},
                                   subexpr_linearity::Vector{JuMP._Derivatives.Linearity},
                                   tag::T) where T <: RelaxTag

    nd = copy(func.nd)
    adj = copy(func.adj)
    const_values = copy(func.const_values)

    # sorted by JuMP, _FunctionStorage but amalagation of the sparsity of the
    # function and the nonlinear expressions is not necessarily
    # so we
    grad_sparsity = copy(func.grad_sparsity)
    for nd in func.nd
        if nd.nodetype === JuMP._Derivatives.SUBEXPRESSION
            append!(grad_sparsity, sub_sparsity_dict[nd.index])
        end
    end
    unique!(grad_sparsity)
    sort!(grad_sparsity)

    reverse_sparsity_length = grad_sparsity[end]
    reverse_sparsity = zeros(Int64, reverse_sparsity_length)
    current_grad_sparsity = grad_sparsity[1]
    current_grad_sparsity_count = 1
    for i = 1:reverse_sparsity_length
        if i == current_grad_sparsity
            reverse_sparsity[i] = current_grad_sparsity_count
            current_grad_sparsity_count += 1
            if current_grad_sparsity_count <= length(grad_sparsity)
                current_grad_sparsity = grad_sparsity[current_grad_sparsity_count]
            else
                break
            end
        end
    end

    N = length(grad_sparsity)
    dependent_variable_count = length(grad_sparsity)

    lenx = length(nd)
    setstorage = fill(MC{N,T}(Interval(-Inf, Inf)), lenx)
    numberstorage = zeros(lenx)
    isnumber = fill(false, lenx)

    tpdict = Dict{Int64,Tuple{Int64,Int64,Int64,Int64}}()
    tp1_count = 0
    tp2_count = 0
    for i = 1:lenx
        node = @inbounds nd[i]
        op = node.index
        if double_tp(op)
            tp1_count += 1
            tp2_count += 1
            tpdict[i] = (tp1_count, tp1_count, tp2_count, tp2_count)
        elseif single_tp(op)
            tp1_count += 1
            tpdict[i] = (tp1_count, tp1_count, -1, -1)
        end
    end

    for (i,node) in enumerate(nd)
        if node.nodetype === JuMP._Derivatives.VALUE
            indx = node.index
            numberstorage[i] = const_values[indx]
        end
    end

    tp1storage = zeros(tp1_count)
    tp2storage = zeros(tp1_count)
    tp3storage = zeros(tp2_count)
    tp4storage = zeros(tp2_count)

    dependent_subexpressions = copy(func.dependent_subexpressions)
    dependent_subexpression_count = length(dependent_subexpressions)

    linearity = JuMP._Derivatives.classify_linearity(nd, adj, subexpr_linearity)

    cv_grad_buffer = zeros(dependent_variable_count)
    cc_grad_buffer = zeros(dependent_variable_count)

    expression = NonlinearExpression{MC{N,T}}(nd, adj, const_values, setstorage, numberstorage,
                                              isnumber, zero(MC{N,T}), false,
                                              tp1storage, tp2storage, tp3storage, tp4storage,
                                              tpdict, grad_sparsity, reverse_sparsity,
                                              dependent_variable_count, dependent_subexpression_count,
                                              dependent_subexpressions, JuMP._Derivatives.CONSTANT,
                                              cv_grad_buffer, cc_grad_buffer)

    saf = SAF(SAT[SAT(0.0, VI(i)) for i = 1:length(grad_sparsity)], 0.0)

    lower_bound = bnds.lower
    upper_bound = bnds.upper

    last_relax_convex = false
    last_relax_concave = false
    last_past_reverse = false
    has_value = false

    return BufferedNonlinearFunction{MC{N,T}}(expression, saf, lower_bound, upper_bound,
                                              last_relax_convex, last_relax_concave,
                                              last_past_reverse, has_value)
end

function BufferedNonlinearFunction()
    return BufferedNonlinearFunction{MC{1,NS}}(NonlinearExpression(), SAF(SAT[], 0.0),
                                               -Inf, Inf, false, false, false, false)
end

function set_intersect_value!(expr::NonlinearExpression{V}, value) where V
    if !expr.isnumber[1]
        expr.value = expr.setstorage[1] ∩ value
        expr.setstorage[1] = expr.value
    end

    return nothing
end

function set_node_flag!(f::BufferedNonlinearFunction{V}) where V
    f.last_relax_convex = false
    f.last_relax_concave = false

    return nothing
end

###
### Defines evaluator storage structure
###
"""
    Evaluator

MOI.AbstractNLPEvaluator for calculating relaxations of nonlinear terms.

Checks that the resulting value should be a number...

$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct Evaluator{N,T} <: MOI.AbstractNLPEvaluator

    user_operators::JuMP._Derivatives.UserOperatorRegistry = JuMP._Derivatives.UserOperatorRegistry()
    has_user_mv_operator::Bool = false
    num_mv_buffer::Vector{T} = T[]
    parameter_values::Vector{T} = T[]

    current_node::NodeBB{N,T} = NodeBB{N,T}()
    treat_x_as_number = Bool[]
    lower_variable_bound::Vector{T} = T[]
    upper_variable_bound::Vector{T} = T[]
    x::Vector{T} = T[]
    node_to_variable_map::Vector{Int} = Int[]
    variable_to_node_map::Vector{Int} = Int[]

    variable_count::Int = 0
    node_count::Int = 0

    "Context used to guard against domain violations & branch on these violations if necessary"
    subgrad_tighten::Bool = false
    reverse_subgrad_tighten::Bool = false

    subexpressions::Vector{NonlinearExpression} = NonlinearExpression[]
    subexpressions_eval::Vector{Bool} = Bool[]

    is_post::Bool = false
    is_intersect::Bool = false
    is_first_eval::Bool = false
    interval_intersect::Bool = false

    subgrad_tol::T = 1E-10
end

"""
$(FUNCTIONNAME)

Sets the current node in the Evaluator structure.
"""
function set_node!(evaluator::Evaluator, n::NodeBB)

    evaluator.current_node = NodeBB(n)
    node_to_variable_map = evaluator.node_to_variable_map
    node_lower_bounds = n.lower_variable_bound
    node_upper_bounds = n.upper_variable_bound
    eval_lower_bounds = evaluator.lower_variable_bound
    eval_upper_bounds = evaluator.upper_variable_bound

    for i = 1:length(evaluator.current_node)
        full_variable_index = node_to_variable_map[i]
        eval_lower_bounds[full_variable_index] = node_lower_bounds[i]
        eval_upper_bounds[full_variable_index] = node_upper_bounds[i]
    end
    fill!(evaluator.subexpressions_eval, false)
    evaluator.is_first_eval = true

    return nothing
end

function retrieve_node(d::Evaluator)
    cn = d.current_node
    node_to_variable_map = d.node_to_variable_map

    return NodeBB(copy(d.lower_variable_bound[node_to_variable_map]),
                  copy(d.upper_variable_bound[node_to_variable_map]),
                  cn.lower_bound, cn.upper_bound, cn.depth, cn.id)
end

function retrieve_x!(out::Vector{Float64}, d::Evaluator)
    x = d.x
    node_to_variable_map = d.node_to_variable_map
    for i in 1:length(node_to_variable_map)
        vindx = node_to_variable_map[i]
        out[i] = x[vindx]
    end

    return nothing
end

# Returns false if subexpression has been evaluated at current reference point
prior_eval(d::Evaluator, i::Int64) = d.subexpressions_eval[i]

#=
Assumes the sparsities are sorted...
=#
function copy_subexpression_value!(k::Int, op::Int, subexpression::NonlinearExpression{MC{N1,T}},
                                   numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N2,T}},
                                   cv_buffer::Vector{Float64}, cc_buffer::Vector{Float64},
                                   func_sparsity::Vector{Int64}) where {N1, N2, T <: RelaxTag}

    # fill cv_grad/cc_grad buffers
    sub_sparsity = subexpression.grad_sparsity
    sset = subexpression.setstorage[1]
    fill!(cv_buffer, 0.0)
    fill!(cc_buffer, 0.0)

    sub_sparsity_count = 1
    subs_index = @inbounds sub_sparsity[1]
    for i = 1:length(func_sparsity)
        func_index = @inbounds func_sparsity[i]
        if func_index === subs_index
            @inbounds cv_buffer[i] = sset.cv_grad[sub_sparsity_count]
            @inbounds cc_buffer[i] = sset.cc_grad[sub_sparsity_count]
            sub_sparsity_count += 1
            subs_index = @inbounds sub_sparsity[sub_sparsity_count]
        end
    end

    cv_grad = SVector(cv_buffer)
    cc_grad = SVector(cc_buffer)

    setstorage[k] = MC{N1,T}(sset.cv, sset.cc, sset.Intv, cv_grad, cc_grad, sset.cnst)

    return nothing
end

function eliminate_fixed_variables!(f::NonlinearExpression{V}, v::Vector{VariableInfo}) where V
    num_constants = length(f.const_values)
    indx_to_const_loc = Dict{Int,Int}()
    for i = 1:length(expr.nd)
        nd = @inbounds expr.nd[i]
        # Assumes MOI Variable have been eliminated previously...
        if nd.nodetype === JuMP._Derivatives.VARIABLE
            indx = nd.index
            if v[indx].is_fixed
                if haskey(indx_to_const_loc, indx)
                    const_loc = indx_to_const_loc[indx]
                    expr.nd[i] = NodeData(JuMP._Derivatives.VALUE, const_loc, nd.parent)
                else
                    push!(const_values, v[indx].lower_bound)
                    num_constants += 1
                    indx_to_const_loc[indx] = num_constants
                    expr.nd[i] = NodeData(nd.nodetype, num_constants, nd.parent)
                end
                f.isnumber[i] = true
            end
        end
    end

    return nothing
end

function eliminate_fixed_variables!(f::BufferedNonlinearFunction{V}, v::Vector{VariableInfo}) where V
    eliminate_fixed_variables!(f.expr, v)
end

function set_evaluator_flags!(d, is_post, is_intersect, is_first_eval, interval_intersect)

    d.is_post = is_post
    d.is_intersect = is_intersect
    d.is_first_eval = is_first_eval
    d.interval_intersect = interval_intersect

    return nothing
end

ext: jl


text:
using PyCall
@pyimport matplotlib.pyplot as plt

srand(0)
include("generate_data.jl")

M = 100
N = 100

# A "blurring" matrix
my_matrix = zeros((M, N))
# Note: Julia uses column-major order
for(j in 1:N)
	for(i in 1:M)
		my_matrix[i, j] = exp(-abs(i - j)/5.0)
	end
end

# A true signal that we want to recover
x = zeros(N)
x[Int(N/2):Int(N/2+5)] = 1.0

# A dataset
y = generate_data(x, my_matrix)

# Do regularised inversion
include("reginv.jl")
x_fit = reginv(y, my_matrix, 1.0)

# Plot the data and the inferred signal using matplotlib
plt.plot(x, "ro-", label="True signal")
plt.plot(y, "bo-", label="Blurred, noisy data")
plt.plot(x_fit, "go-", label="MAP")
plt.legend(loc="upper right")
plt.show()

max_stars_repo_name: eggplantbren/MeanImage

max_issues_repo_name: eggplantbren/MeanImage

max_forks_repo_name: eggplantbren/MeanImage


text:
mutable struct DescriptorSetLayoutCreateInfo
    mBindings::Vector{vk.VkDescriptorSetLayoutBinding}

    mHandleRef::Ref{vk.VkDescriptorSetLayoutCreateInfo}

    function DescriptorSetLayoutCreateInfo()
        this = new()
        this.mBindings = Vector{vk.VkDescriptorSetLayoutBinding}()
        return this
    end
end

function addBinding(this::DescriptorSetLayoutCreateInfo,
                    type::vk.VkDescriptorType,
                   stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL_GRAPHICS,
                   count::UInt32 = 0,
                   samplers::Ptr{vk.VkSampler} = Ptr{vk.VkSampler}(C_NULL))
    push!(this.mBindings, vk.VkDescriptorSetLayoutBinding(
                                length(this.mBindings), #binding::UInt32
                                type, #descriptorType::VkDescriptorType
                                count, #descriptorCount::UInt32
                                stage, #stageFlags::VkShaderStageFlags
                                samplers, #pImmutableSamplers::Ptr{VkSampler}
                            ))
end

function addSampler(this::DescriptorSetLayoutCreateInfo,
                   stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL_GRAPHICS,
                   count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_SAMPLER, stage, count)
end

function addUniformBuffer(this::DescriptorSetLayoutCreateInfo,
                         stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL,
                         count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER, stage, count)
end

function addCombinedImageSampler(this::DescriptorSetLayoutCreateInfo,
                                stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL_GRAPHICS,
                                count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, stage, count)
end

function addSampledImage(this::DescriptorSetLayoutCreateInfo,
                        stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL_GRAPHICS,
                        count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_SAMPLED_IMAGE, stage, count)
end

function addStorageBuffer(this::DescriptorSetLayoutCreateInfo,
                         stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL,
                         count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, stage, count)
end

function addStorageImage(this::DescriptorSetLayoutCreateInfo,
                        stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL,
                        count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_STORAGE_IMAGE, stage, count)
end

function addAccelerationStructure(this::DescriptorSetLayoutCreateInfo,
                                 stage::vk.VkShaderStageFlagBits = vk.VkShaderStageFlagBits(
                                        vk.VK_SHADER_STAGE_RAYGEN_BIT_NV |
                                        vk.VK_SHADER_STAGE_ANY_HIT_BIT_NV |
                                        vk.VK_SHADER_STAGE_CLOSEST_HIT_BIT_NV |
                                        vk.VK_SHADER_STAGE_MISS_BIT_NV),
                                 count::UInt32 = UInt32(1))
    addBinding(this, vk.VK_DESCRIPTOR_TYPE_ACCELERATION_STRUCTURE_NV, stage, count)
end

function handleRef(this::DescriptorSetLayoutCreateInfo)::Ref{vk.VkDescriptorSetLayoutCreateInfo}
    this.mHandleRef = Ref(vk.VkDescriptorSetLayoutCreateInfo(
        vk.VK_STRUCTURE_TYPE_DESCRIPTOR_SET_LAYOUT_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        0, #flags::VkDescriptorSetLayoutCreateFlags
        length(this.mBindings), #bindingCount::UInt32
        pointer(this.mBindings) #pBindings::Ptr{VkDescriptorSetLayoutBinding}
    ))
    return this.mHandleRef
end


text:
using Reexport

@reexport module CausalityTests
    import CausalityToolsBase: causality
    import NearestNeighbors
    import StatsBase
    import Distances
    import CausalityToolsBase
    import CausalityToolsBase: CausalityTest, RectangularBinning
    import CrossMappings: crossmap, convergentcrossmap
    import TimeseriesSurrogates: randomshuffle
    import TransferEntropy: transferentropy, TransferOperatorGrid, VisitationFrequency
    import UncertainData
    import UncertainData:
        AbstractUncertainValue,
        AbstractUncertainValueDataset,
        AbstractUncertainIndexValueDataset,
        resample,
        ConstrainedResampling

    resample(v::Vector{Real}) = v

    """
        causality(source::AbstractVector, target::AbstractVector, test::CausalityTest)

    Test for a causal influence from `source` to `target` using the provided causality `test`.

    ## Examples

    ```julia
    x, y = rand(300), rand(300)

    # Define some causality tests and apply them to `x` and `y`.
    test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
    test_cm = CrossMappingTest(n_reps = 10)
    test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
    test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
    test_jdd = JointDistanceDistributionTest()
    test_jddt = JointDistanceDistributionTTest()
    predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
    test_pa = PredictiveAsymmetryTest(predictive_test = predtest)

    causality(x, y, test_ccm)
    causality(x, y, test_cm)
    causality(x, y, test_vf)
    causality(x, y, test_tog)
    causality(x, y, test_jdd)
    causality(x, y, test_jddt)
    causality(x, y, test_pa)
    ```
    """
    function causality(x::AbstractVector, y::AbstractVector, test::CausalityTest) end

    function summarise(test::CausalityTest)
        _type = typeof(test)

        strs = ["$fn = $(test.:($fn))" for fn in fieldnames(_type)]
        return "$_type" * "(" * join(strs, ", ") * ")"
    end

    Base.show(io::IO, test::CausalityTest) = print(io, summarise(test))

    ################################################################
    # Distance based causality tests
    ################################################################
    include("tests_distance_based/DistanceBasedCausalityTest.jl")

    # S-measure test
    # ---------------------
    import ..s_measure

    include("tests_distance_based/SMeasureTest.jl")

    # Joint distances tests
    # ---------------------
    import ..joint_distance_distribution

    include("tests_distance_based/JointDistancesCausalityTest.jl")
    include("tests_distance_based/JointDistanceDistributionTest.jl")
    include("tests_distance_based/JointDistanceDistributionTTest.jl")

    # Cross mapping tests
    # ---------------------
    include("tests_distance_based/CrossMappingTest.jl")
    include("tests_distance_based/ConvergentCrossMappingTest.jl")

    ################################################################
    # Entropy based causality tests
    ################################################################
    include("tests_entropy_based/EntropyBasedCausalityTest.jl")

    # Transfer entropy causality tests
    # ---------------------------------------
    include("tests_entropy_based/TransferEntropyCausalityTest.jl")
    include("tests_entropy_based/VisitationFrequencyTest.jl")
    include("tests_entropy_based/TransferOperatorGridTest.jl")
    include("tests_entropy_based/ApproximateSimplexIntersectionTest.jl")
    include("tests_entropy_based/ExactSimplexIntersectionTest.jl")

    ################################################################
    # Predictive asymmetry causality tests
    ################################################################
    include("tests_predictive_asymmetry/PredictiveAsymmetryTest.jl")

    ################################################################
    # On uncertain data with uncertainties in both index and value
    ################################################################
    include("causality_on_uncertaindata.jl")

    ################################################################
    # Preprocessed data + causality test
    ################################################################
    include("highlevel_tests/PreprocessedDataCausalityTest.jl")
    include("highlevel_tests/BinnedDataCausalityTest.jl")
    include("highlevel_tests/ConstrainedTest.jl")
    include("highlevel_tests/InterpolateBinTest.jl")
    include("highlevel_tests/RandomSequencesTest.jl")

    include("highlevel_tests/causality_BinnedDataCausalityTest.jl")
    include("highlevel_tests/causality_ConstrainedTest.jl")
    include("highlevel_tests/causality_InterpolateBinTest.jl")
    include("highlevel_tests/causality_RandomSequencesTest.jl")

    export causality
end

"""
	CausalityTests

A module defining causality tests.
"""
CausalityTests


text:
using Documenter, BitSAD

makedocs(;
    modules=[BitSAD],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Introduction" => [
            "Getting started" => "getting-started.md",
            "Stochastic Bitstream Walkthrough" => "sbitstream-example.md",
            "Deterministic Bitstream Walkthrough" => "dbitstream-example.md",
        ],
        "Types" => [
            "Abstract Bitstreams" => "bitstream.md",
            "Stochastic Bitstreams" => "sbitstream.md",
            "Deterministic Bitstreams" => "dbitstream.md",
        ],
        "Hardware Generation" => [
            "Getting started" => "hardware-generation.md",
            "Internals" => "hardware-internals.md",
        ],
        "Customizing BitSAD" => [
            "Custom hardware generation" => "custom-hardware.md",
            "Custom `SBitstream` operators" => "custom-soperators.md",
        ],
    ],
    repo="https://github.com/UW-PHARM/BitSAD.jl/blob/{commit}{path}#L{line}",
    sitename="BitSAD.jl",
    authors="Kyle Daruwalla, University of Wisconsin-Madison PHARM Group",
    assets=String[],
)

deploydocs(;
    repo="github.com/UW-PHARM/BitSAD.jl",
)


text:
module SeisImaging
    using Compat
    include("Imaging/Imaging.jl")
end


text:
using StatsBase

export summarize

abstract type Summary end

summarize(nt::NamedTuple) = rmap(summarize, nt)

summarize(tv::TupleVector) = summarize(unwrap(tv))

import Base

###############################################################################
# FullSummary

struct FullSummary{T} <: Summary
    data :: T
end

summarize(x) = FullSummary(x)

Base.show(io, s::FullSummary) = print(io, s.data)

###############################################################################
# RealSummary
export RealSummary
struct RealSummary <: Summary
    μ :: Float64
    σ :: Float64

    function RealSummary(μ,σ)
        μ = isnan(μ) ? 0.0 : μ
        σ = isnan(σ) ? 0.0 : σ
        new(μ,σ)
    end
end

Base.typeinfo_prefix(io::IO, ::AbstractArray{<:Summary}) = ("", false)

function Base.show(io::IO, s::RealSummary)
    io = IOContext(io, :compact => true)
    σ = round(s.σ, sigdigits=2)
    if s.μ == 0 || σ == 0
        μdigits = 2
    else
        μdigits = max(0, ceil(Int, log10(2) * (exponent(s.μ) - exponent(σ))) + 2)
    end
    μ = round(s.μ, sigdigits = μdigits)
    print(io, μ, "±", σ)
end

summarize(x::AbstractVector{<:Real}) = RealSummary(mean_and_std(x)...)

function summarize(x::AbstractArray{<:AbstractArray{T}}) where {T <: Real}
    [RealSummary(μ,σ) for (μ,σ) in zip(mean_and_std(x)...)]
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(abs_tol=1e-15, orderT=8, orderQ=1, max_steps=50_000)),
    "X0"   => Hyperrectangle(low=[0.8, 0.8], high=[1.2, 1.2]),
    "info" => Dict(
        "name" => "jet_engine",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])


text:
#
# Select atoms using vmd selection syntax, with vmd in background
#

# Method to return the selection from the topology file, by calling VMD in the background.

function select( topology :: String, selection :: String; vmd="vmd" )

  local index_list :: String
  local readnext :: Bool = false

  vmd_input = Base.open("./VMDINPUT_TMP.VMD","w")
  Base.write(vmd_input,"mol new \"$topology\" \n")
  Base.write(vmd_input,"set sel [ atomselect top \"$selection\" ] \n")
  Base.write(vmd_input,"puts \"INDEXLIST\" \n")
  Base.write(vmd_input,"set indexes [ \$sel get index ] \n")
  Base.write(vmd_input,"puts \"ENDINDEXLIST\" \n")
  Base.write(vmd_input,"exit \n")
  Base.close(vmd_input)

  vmd_output = read(`$vmd -dispdev text -e ./VMDINPUT_TMP.VMD`, String)

  for line in split(vmd_output,"\n")
    if readnext
      if line == "ENDINDEXLIST"
        error("ERROR: Selection '$selection' does not contain any atom")
      end
      index_list = line
      break
    end
    if line == "INDEXLIST"
      readnext = true
    end
  end
  index_split = split(index_list)
  nsel = length(index_split)
  selection_indexes = Vector{Int64}(undef,nsel)
  for i in 1:nsel
    selection_indexes[i] = parse(Int64,index_split[i]) + 1
  end

  run(`\rm -f ./VMDINPUT_TMP.VMD`)

  return selection_indexes

end

# Method to return the selection if what was provided was the Simulation object.

select(simulation :: Simulation, selection :: String) = select( simulation.psf, selection, vmd=simulation.vmd )

#
# select atoms using VMD selection syntax
#

#function select(atoms :: Vector{PDBEnergy.Atom}, selection :: String)
#
#  sel = selection
#
#  keys = [ "name" ]
#
#  seldata = split(sel)
#  code = "("
#  for str in seldata
#    if str in keys
#      str = replace(sel, key => "atom.$key in [" )
#    else
#
#    end
#    code = "$code $str"
#  end
#
#  code = "[ atom.index for atom in filter( atom -> ( $sel ), atoms ) ]"
#
#  println(code)
#
#
#  #return selection_indexes
#end


text:
using StochasticDiffEq, DiffEqDevTools, Test, LinearAlgebra
Random.seed!(100)
dts = (1/2) .^ (10:-1:2) #14->7 good plot

const σ_const = 0.87
const μ = 1.01

u0 = rand(2)
A = [-2.0 1.0;1.0 -2.0]
B = Diagonal([σ_const for i in 1:2])

function f_commute(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end

function f_commute_analytic(u0,p,t,W)
 tmp = (A+1.01I-(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

function σ(du,u,p,t)
  du[1,1] = σ_const*u[1]
  du[1,2] = σ_const*u[1]
  du[2,1] = σ_const*u[2]
  du[2,2] = σ_const*u[2]
end

ff_commute = SDEFunction(f_commute,σ,analytic=f_commute_analytic)

prob = SDEProblem(ff_commute,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,2))

sol = solve(prob,RKMilCommute(),dt=1/2^(8))
sol = solve(prob,EM(),dt=1/2^(10))

dts = (1/2) .^ (10:-1:3) #14->7 good plot
sim2 = test_convergence(dts,prob,EM(),numMonte=Int(1e2))
sim2 = test_convergence(dts,prob,RKMilCommute(),numMonte=Int(1e2))

sim2.𝒪est[:final] - 1 < 0.2


text:
struct EquationSolveError
    eq
    var
    rhs
end

function Base.showerror(io::IO, ese::EquationSolveError)
    print(io, "EquationSolveError: While solving\n\n\t")
    print(io, ese.eq)
    print(io, "\nfor ")
    printstyled(io, var, bold=true)
    print(io, ", obtained RHS\n\n\tt")
    println(io, rhs)
end

function masked_cumsum!(A::Vector)
    acc = zero(eltype(A))
    for i in eachindex(A)
        iszero(A[i]) && continue
        A[i] = (acc += A[i])
    end
end

function contract_variables(graph::BipartiteGraph, var_eq_matching::Matching, eliminated_variables)
    var_rename = ones(Int64, ndsts(graph))
    eq_rename = ones(Int64, nsrcs(graph))
    for v in eliminated_variables
        eq_rename[var_eq_matching[v]] = 0
        var_rename[v] = 0
    end
    masked_cumsum!(var_rename)
    masked_cumsum!(eq_rename)

    dig = DiCMOBiGraph{true}(graph, var_eq_matching)

    # Update bipartite graph
    var_deps = map(1:ndsts(graph)) do v
        [var_rename[v′] for v′ in neighborhood(dig, v, Inf; dir=:in) if var_rename[v′] != 0]
    end

    nelim = length(eliminated_variables)
    newgraph = BipartiteGraph(nsrcs(graph) - nelim, ndsts(graph) - nelim)
    for e in 𝑠vertices(graph)
        ne = eq_rename[e]
        ne == 0 && continue
        for v in 𝑠neighbors(graph, e)
            newvar = var_rename[v]
            if newvar != 0
                add_edge!(newgraph, ne, newvar)
            else
                for nv in var_deps[v]
                    add_edge!(newgraph, ne, nv)
                end
            end
        end
    end

    return newgraph
end

"""
    algebraic_variables_scc(sys)

Find strongly connected components of algebraic variables in a system.
"""
function algebraic_variables_scc(state::TearingState)
    graph = state.structure.graph
    # skip over differential equations
    algvars = BitSet(findall(v->isalgvar(state.structure, v), 1:ndsts(graph)))
    algeqs = BitSet(findall(map(1:nsrcs(graph)) do eq
        all(v->!isdervar(state.structure, v), 𝑠neighbors(graph, eq))
    end))
    var_eq_matching = complete(maximal_matching(graph, e->e in algeqs, v->v in algvars))
    var_sccs = find_var_sccs(complete(graph), var_eq_matching)

    return var_eq_matching, var_sccs
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using Klara

p = BasicContUnvParameter(:p, logtarget=p::Real -> -abs2(p), diffopts=DiffOptions(mode=:forward))

model = likelihood_model([p], isindexed=false)

mcsampler = NUTS(0.4, maxndoublings=7)

tuner = DualAveragingMCTuner(0.651, 5000)

mcrange = BasicMCRange(nsteps=5000, burnin=1000)

v0 = Dict(:p=>3.11)

outopts = Dict(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept, :ndoublings, :a, :na])

job = BasicMCJob(model, mcsampler, mcrange, v0, tuner=tuner, outopts=outopts)

@time run(job)

chain = output(job)

mean(chain)

diags = diagnostics(chain)

mean(diags[:a]./diags[:na])


text:
# The file contains test examples definitions

#=
    The root-n system
    (not to be confused with cyclic system)!!
=#
function rootn(n; ground=QQ)
    R, xs = PolynomialRing(ground, ["x$i" for i in 1:n])
    ans = [
        sum(map(prod, Combinatorics.combinations(xs, i)))
        for i in 1:n
    ]
    ans[end] -= (-1)^(n - 1)
    ans
end

#=
    The cyclic-n system
    (not to be confused with cyclic system)!!
=#
function cyclicn(n; ground=QQ)
    R, z = PolynomialRing(ground, ["z$i" for i in 1:n])
    [(sum(prod(z[(k-1) % n + 1] for k in j:j+m) for j in 1:n) for m=0:(n-2))...,prod(z)-1]
end

#=
    The katsura-n system
=#
function katsuran(n; ground=QQ)
    R, x = PolynomialRing(ground, ["x$i" for i in 0:n])

    return [
        (sum(x[abs(l)+1]*x[abs(m-l)+1] for l=-n:n if abs(m-l)<=n) -
        x[m+1] for m=0:n-1)...,
        x[1] + 2sum(x[i+1] for i=1:n) - 1
    ]
end

function henrion5(;ground=QQ)

    R, (f1,f2,f3,f4,f5,t) = PolynomialRing(ground, ["f1","f2","f3","f4","f5","t"])
    fs = [
        2*f1*f2*f3*f4*f5-9823275,
        ground(21)//5*f1*f2*f4*f5+ground(16)//5*f1*f3*f4*f5+ground(9)//5*f2*f3*f4*f5+ground(24)//5*f1*f2*f3*f5+5*f4*f3*f1*f2-4465125,
        ground(14)//5*f4*f5*f1+ground(14)//5*f4*f5*f2+ground(8)//5*f3*f4*f5+ground(18)//5*f1*f2*f5+ground(24)//5*f1*f3*f5+ground(18)//5*f2*f3*f5+4*f3*f1*f2+6*f1*f2*f4+6*f3*f4*f1+4*f2*f3*f4-441486,
        ground(7)//5*f4*f5+ground(12)//5*f5*f1+ground(12)//5*f5*f2+ground(12)//5*f5*f3+3*f1*f2+4*f3*f1+4*f4*f1+3*f2*f3+4*f4*f2+3*f3*f4-15498,
        ground(6)//5*f5+2*f4+2*f3+2*f2+2*f1-215,
        f1+2*f2+3*f3+4*f4+5*f5+6*t
    ]
end

function henrion6(;ground=QQ)
    R, (f1,f2,f3,f4,f5,f6) = PolynomialRing(ground, ["f1","f2","f3","f4","f5","f6"])

    fs = [
        2*f1*f2*f3*f4*f5*f6-1404728325,
        6*f5*f4*f3*f1*f2+ground(11)//6*f2*f3*f4*f5*f6+ground(16)//3*f1*f2*f3*f5*f6+ground(9)//2*f1*f2*f4*f5*f6+ground(10)//3*f1*f3*f4*f5*f6+ground(35)//6*f1*f2*f3*f4*f6-648336150,
        5*f4*f3*f1*f2+5*f2*f3*f4*f5+ground(5)//3*f3*f4*f5*f6+8*f1*f2*f3*f5+9*f1*f2*f4*f5+8*f1*f3*f4*f5+4*f1*f2*f5*f6+ground(16)//3*f1*f3*f5*f6+3*f1*f4*f5*f6+4*f2*f3*f5*f6+3*f2*f4*f5*f6+ground(14)//3*f1*f2*f3*f6+7*f1*f2*f4*f6+7*f1*f3*f4*f6+ground(14)//3*f2*f3*f4*f6-67597623,
        6*f1*f2*f5+8*f1*f3*f5+6*f2*f3*f5+ground(8)//3*f5*f6*f3+ground(8)//3*f5*f6*f2+ground(8)//3*f5*f6*f1+ground(7)//2*f1*f2*f6+ground(14)//3*f1*f3*f6+ground(14)//3*f1*f4*f6+ground(7)//2*f2*f3*f6+ground(14)//3*f2*f4*f6+ground(7)//2*f3*f4*f6+6*f4*f5*f1+ground(3)//2*f4*f5*f6+4*f3*f1*f2+4*f2*f3*f4+6*f3*f4*f1+4*f3*f4*f5+6*f1*f2*f4+6*f4*f5*f2-2657700,
        ground(4)//3*f5*f6+ground(7)//3*f6*f1+ground(7)//3*f6*f2+ground(7)//3*f6*f3+ground(7)//3*f6*f4+3*f1*f2+4*f3*f1+4*f4*f1+4*f5*f1+3*f2*f3+4*f4*f2+4*f5*f2+3*f3*f4+4*f5*f3+3*f4*f5-46243,
        ground(7)//6*f6+2*f5+2*f4+2*f3+2*f2+2*f1-358
    ]
end

function henrion7(;ground=QQ)
    R, (f1,f2,f3,f4,f5,f6,f7) = PolynomialRing(ground, ["f1","f2","f3","f4","f5","f6","f7"])

    fs = [
        2*f1*f2*f3*f4*f5*f6*f7-273922023375,
        ground(45)//7*f1*f2*f3*f4*f6*f7+ground(40)//7*f1*f2*f3*f5*f6*f7+ground(33)//7*f1*f2*f4*f5*f6*f7+ground(24)//7*f1*f3*f4*f5*f6*f7+ground(48)//7*f1*f2*f3*f4*f5*f7+7*f6*f5*f4*f3*f1*f2+ground(13)//7*f2*f3*f4*f5*f6*f7-127830277575,
        6*f5*f4*f3*f1*f2+6*f2*f3*f4*f5*f6+ground(12)//7*f3*f4*f5*f6*f7+10*f1*f2*f3*f4*f6+12*f1*f2*f3*f5*f6+12*f1*f2*f4*f5*f6+10*f1*f3*f4*f5*f6+ground(36)//7*f1*f2*f3*f6*f7+ground(54)//7*f1*f2*f4*f6*f7+ground(30)//7*f1*f2*f5*f6*f7+ground(54)//7*f1*f3*f4*f6*f7+ground(40)//7*f1*f3*f5*f6*f7+ground(22)//7*f1*f4*f5*f6*f7+ground(36)//7*f2*f3*f4*f6*f7+ground(30)//7*f2*f3*f5*f6*f7+ground(22)//7*f2*f4*f5*f6*f7+ground(40)//7*f1*f2*f3*f4*f7+ground(64)//7*f1*f2*f3*f5*f7+ground(72)//7*f1*f2*f4*f5*f7+ground(64)//7*f1*f3*f4*f5*f7+ground(40)//7*f2*f3*f4*f5*f7-13829872635,
        -585849123+5*f4*f3*f1*f2+5*f2*f3*f4*f5+5*f3*f4*f5*f6+8*f1*f2*f3*f5+9*f1*f2*f4*f5+8*f1*f3*f4*f5+ground(11)//7*f4*f5*f6*f7+8*f1*f2*f3*f6+12*f1*f2*f4*f6+9*f1*f2*f5*f6+12*f1*f3*f4*f6+12*f1*f3*f5*f6+8*f1*f4*f5*f6+8*f2*f3*f4*f6+9*f2*f3*f5*f6+8*f2*f4*f5*f6+ground(27)//7*f1*f2*f6*f7+ground(36)//7*f1*f3*f6*f7+ground(36)//7*f1*f4*f6*f7+ground(20)//7*f1*f5*f6*f7+ground(27)//7*f2*f3*f6*f7+ground(36)//7*f2*f4*f6*f7+ground(20)//7*f2*f5*f6*f7+ground(27)//7*f3*f4*f6*f7+ground(20)//7*f3*f5*f6*f7+ground(32)//7*f1*f2*f3*f7+ground(48)//7*f1*f2*f4*f7+ground(48)//7*f1*f2*f5*f7+ground(48)//7*f1*f3*f4*f7+ground(64)//7*f1*f3*f5*f7+ground(48)//7*f1*f4*f5*f7+ground(32)//7*f2*f3*f4*f7+ground(48)//7*f2*f3*f5*f7+ground(48)//7*f2*f4*f5*f7+ground(32)//7*f3*f4*f5*f7,
        -11675085+6*f1*f2*f6+8*f1*f3*f6+8*f1*f4*f6+6*f2*f3*f6+8*f2*f4*f6+6*f3*f4*f6+ground(18)//7*f6*f7*f4+ground(18)//7*f6*f7*f3+ground(18)//7*f6*f7*f2+ground(18)//7*f6*f7*f1+ground(24)//7*f1*f2*f7+ground(32)//7*f1*f3*f7+ground(32)//7*f1*f4*f7+ground(32)//7*f1*f5*f7+ground(24)//7*f2*f3*f7+ground(32)//7*f2*f4*f7+ground(32)//7*f2*f5*f7+ground(24)//7*f3*f4*f7+ground(32)//7*f3*f5*f7+ground(24)//7*f4*f5*f7+
        ground(10)//7*f5*f6*f7+6*f1*f2*f5+8*f1*f3*f5+6*f2*f3*f5+6*f5*f6*f3+6*f5*f6*f2+6*f5*f6*f1+6*f1*f2*f4+6*f4*f5*f2+6*f4*f5*f1+4*f4*f5*f6+6*f3*f4*f1+4*f3*f4*f5+4*f3*f1*f2+4*f2*f3*f4,
        ground(9)//7*f6*f7+ground(16)//7*f7*f1+ground(16)//7*f7*f2+ground(16)//7*f7*f3+ground(16)//7*f7*f4+ground(16)//7*f7*f5+3*f1*f2+4*f3*f1+4*f4*f1+4*f5*f1+4*f6*f1+3*f2*f3+4*f4*f2+4*f5*f2+4*f6*f2+3*f3*f4+4*f5*f3+4*f6*f3+3*f4*f5+4*f6*f4+3*f5*f6-116053,
        ground(8)//7*f7+2*f6+2*f5+2*f4+2*f3+2*f2+2*f1-553
    ]
end

function eco5(;ground=GF(2^31-1))
    R, (x1, x2, x3, x4, x5) = PolynomialRing(ground, ["x$i" for i in 1:5])

    fs = [
    (x1 + x1*x2 + x2*x3 + x3*x4)*x5 - 1,
     (x2 + x1*x3 + x2*x4)*x5 - 2,
             (x3 + x1*x4)*x5 - 3,
                       x4*x5 - 4,
           x1 + x2 + x3 + x4 + 1
    ]
end

function eco7(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7) = PolynomialRing(ground, ["x$i" for i in 1:7])

    fs = [
        (x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6)*x7 - 1,
     (x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6)*x7 - 2,
     (x3 + x1*x4 + x2*x5 + x3*x6)*x7 - 3,
     (x4 + x1*x5 + x2*x6)*x7 - 4,
     (x5 + x1*x6)*x7 - 5,
     x6*x7 - 6,
     x1 + x2 + x3 + x4 + x5 + x6 + 1
    ]
    fs
end

function eco10(;ground=QQ)
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x0*x1*x9+x1*x2*x9+x2*x3*x9+x3*x4*x9+x4*x5*x9+x5*x6*x9+x6*x7*x9+x7*x8*x9+x0*x9-1,
        x0*x2*x9+x1*x3*x9+x2*x4*x9+x3*x5*x9+x4*x6*x9+x5*x7*x9+x6*x8*x9+x1*x9-2,
        x0*x3*x9+x1*x4*x9+x2*x5*x9+x3*x6*x9+x4*x7*x9+x5*x8*x9+x2*x9-3,
        x0*x4*x9+x1*x5*x9+x2*x6*x9+x3*x7*x9+x4*x8*x9+x3*x9-4,
        x0*x5*x9+x1*x6*x9+x2*x7*x9+x3*x8*x9+x4*x9-5,
        x0*x6*x9+x1*x7*x9+x2*x8*x9+x5*x9-6,
        x0*x7*x9+x1*x8*x9+x6*x9-7,
        x0*x8*x9+x7*x9-8,
        x8*x9-9,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+1
    ]
    fs
end

function eco11(;ground=QQ)
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x0*x1*x10+x1*x2*x10+x2*x3*x10+x3*x4*x10+x4*x5*x10+x5*x6*x10+x6*x7*x10+x7*x8*x10+x8*x9*x10+x0*x10-1,
        x0*x2*x10+x1*x3*x10+x2*x4*x10+x3*x5*x10+x4*x6*x10+x5*x7*x10+x6*x8*x10+x7*x9*x10+x1*x10-2,
        x0*x3*x10+x1*x4*x10+x2*x5*x10+x3*x6*x10+x4*x7*x10+x5*x8*x10+x6*x9*x10+x2*x10-3,
        x0*x4*x10+x1*x5*x10+x2*x6*x10+x3*x7*x10+x4*x8*x10+x5*x9*x10+x3*x10-4,
        x0*x5*x10+x1*x6*x10+x2*x7*x10+x3*x8*x10+x4*x9*x10+x4*x10-5,
        x0*x6*x10+x1*x7*x10+x2*x8*x10+x3*x9*x10+x5*x10-6,
        x0*x7*x10+x1*x8*x10+x2*x9*x10+x6*x10-7,
        x0*x8*x10+x1*x9*x10+x7*x10-8,
        x0*x9*x10+x8*x10-9,
        x9*x10-10,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+1
    ]
    fs
end

function eco12(;ground=QQ)
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11) = PolynomialRing(ground, ["x$i" for i in 1:12])

    fs = [
        x0*x1*x11+x1*x2*x11+x2*x3*x11+x3*x4*x11+x4*x5*x11+x5*x6*x11+x6*x7*x11+x7*x8*x11+x8*x9*x11+x9*x10*x11+x0*x11-1,
        x0*x2*x11+x1*x3*x11+x2*x4*x11+x3*x5*x11+x4*x6*x11+x5*x7*x11+x6*x8*x11+x7*x9*x11+x8*x10*x11+x1*x11-2,
        x0*x3*x11+x1*x4*x11+x2*x5*x11+x3*x6*x11+x4*x7*x11+x5*x8*x11+x6*x9*x11+x7*x10*x11+x2*x11-3,
        x0*x4*x11+x1*x5*x11+x2*x6*x11+x3*x7*x11+x4*x8*x11+x5*x9*x11+x6*x10*x11+x3*x11-4,
        x0*x5*x11+x1*x6*x11+x2*x7*x11+x3*x8*x11+x4*x9*x11+x5*x10*x11+x4*x11-5,
        x0*x6*x11+x1*x7*x11+x2*x8*x11+x3*x9*x11+x4*x10*x11+x5*x11-6,
        x0*x7*x11+x1*x8*x11+x2*x9*x11+x3*x10*x11+x6*x11-7,
        x0*x8*x11+x1*x9*x11+x2*x10*x11+x7*x11-8,
        x0*x9*x11+x1*x10*x11+x8*x11-9,
        x0*x10*x11+x9*x11-10,
        x10*x11-11,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+1
    ]
    fs
end

function eco13(;ground=QQ)
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12) = PolynomialRing(ground, ["x$i" for i in 1:13])

    fs = [
        x0*x1*x12+x1*x2*x12+x2*x3*x12+x3*x4*x12+x4*x5*x12+x5*x6*x12+x6*x7*x12+x7*x8*x12+x8*x9*x12+x9*x10*x12+x10*x11*x12+x0*x12-1,
        x0*x2*x12+x1*x3*x12+x2*x4*x12+x3*x5*x12+x4*x6*x12+x5*x7*x12+x6*x8*x12+x7*x9*x12+x8*x10*x12+x9*x11*x12+x1*x12-2,
        x0*x3*x12+x1*x4*x12+x2*x5*x12+x3*x6*x12+x4*x7*x12+x5*x8*x12+x6*x9*x12+x7*x10*x12+x8*x11*x12+x2*x12-3,
        x0*x4*x12+x1*x5*x12+x2*x6*x12+x3*x7*x12+x4*x8*x12+x5*x9*x12+x6*x10*x12+x7*x11*x12+x3*x12-4,
        x0*x5*x12+x1*x6*x12+x2*x7*x12+x3*x8*x12+x4*x9*x12+x5*x10*x12+x6*x11*x12+x4*x12-5,
        x0*x6*x12+x1*x7*x12+x2*x8*x12+x3*x9*x12+x4*x10*x12+x5*x11*x12+x5*x12-6,
        x0*x7*x12+x1*x8*x12+x2*x9*x12+x3*x10*x12+x4*x11*x12+x6*x12-7,
        x0*x8*x12+x1*x9*x12+x2*x10*x12+x3*x11*x12+x7*x12-8,
        x0*x9*x12+x1*x10*x12+x2*x11*x12+x8*x12-9,
        x0*x10*x12+x1*x11*x12+x9*x12-10,
        x0*x11*x12+x10*x12-11,
        x11*x12-12,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+1
    ]
    fs
end

#=
    The noon-4 system
=#
function noonn(n; ground=QQ)
    without(x, k) = x[1:end .!= k]

    R, xs = PolynomialRing(ground, ["x$i" for i in 1:n])
    fs = zeros(R, n)
    for i in 1:n
        other = without(xs, i)
        fs[i] = xs[i] * (10*sum(other .^ 2) - 11) + 10
    end
    fs
end

function ku10(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) = PolynomialRing(ground, ["x$i" for i in 1:10])

    fs = [
        5*x1*x2+ 5*x1+ 3*x2+ 55,
        7*x2*x3+ 9*x2+ 9*x3+ 19,
        3*x3*x4+ 6*x3+ 5*x4-4,
        6*x4*x5+ 6*x4+ 7*x5+ 118,
        x5*x6+ 3*x5+ 9*x6+ 27,
        6*x6*x7+ 7*x6+x7+ 72,
        9*x7*x8+ 7*x7+x8+ 35,
        4*x8*x9+ 4*x8+ 6*x9+ 16,
        8*x9*x10+ 4*x9+ 3*x10-51,
        3*x1*x10-6*x1+x10+ 5
    ]
    fs
end

function kinema(;ground=QQ)
    R, (z1, z2, z3, z4, z5, z6, z7, z8, z9) = PolynomialRing(ground, ["z$i" for i in 1:9])

    fs = [
    z1^2 + z2^2 + z3^2 - 12*z1 - 68;
    z4^2 + z5^2 + z6^2 - 12*z5 - 68;
    z7^2 + z8^2 + z9^2 - 24*z8 - 12*z9 + 100;
    z1*z4 + z2*z5 + z3*z6 - 6*z1 - 6*z5 - 52;
    z1*z7 + z2*z8 + z3*z9 - 6*z1 - 12*z8 - 6*z9 + 64;
    z4*z7 + z5*z8 + z6*z9 - 6*z5 - 12*z8 - 6*z9 + 32;
    2*z2 + 2*z3 - z4 - z5 - 2*z6 - z7 - z9 + 18;
    z1 + z2 + 2*z3 + 2*z4 + 2*z6 - 2*z7 + z8 - z9 - 38;
    z1 + z3 - 2*z4 + z5 - z6 + 2*z7 - 2*z8 + 8;
    ]
end

function sparse5(; ground=QQ)
    R, (x1, x2, x3, x4, x5) = PolynomialRing(ground, ["x$i" for i in 1:5])

    fs = [
        x1^2*x2^2*x3^2*x4^2*x5^2 + 3*x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + 3*x2^2 + x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + 3*x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + x3^2 + 3*x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + x3^2 + x4^2 + 3*x5^2 + x1*x2*x3*x4*x5 + 5
    ]
end

function s9_1(; ground=QQ)
    R, (a, b, c, d, e, f, g, h) = PolynomialRing(ground, ["x$i" for i in 1:8])

    fs = [-e*g - 2*d*h,
        9*e + 4*b,
        -4*c*h - 2*e*f - 3*d*g,
        -7*c + 9*a - 8*f,
        -4*d*f - 5*c*g - 6*h - 3*e,
        -5*d - 6*c*f - 7*g + 9*b,
        9*d + 6*a - 5*b,
        9*c - 7*a + 8]
end

function ojika4(;ground=QQ)
    R, (x1, x2, x3) = PolynomialRing(ground, ["x$i" for i in 1:3])

    fs = [x1+x3*x1^3+x1*x3*x2^2-x1*x3,
    10*x2-2*x2*x3*x1^2-x3*x2^3-x2*x3,
    -6*x3^2*x1^4-3*x1^2*x2^2*x3^2-x3^2*x1^2+28*x3*x1^2 - 3*x3^2*x2^4+2*x3^2*x2^2+7*x3*x2^2+x3^2-11*x3+10
    ]
end

function ojika3_d1R2(;ground=QQ)
    R, (x1, x2, x3) = PolynomialRing(ground, ["x$i" for i in 1:3])

    fs = [x1^3*x3+x1*x3*x2^2-x1*x3+x1,
-2*x1^2*x3*x2-x3*x2^3-x3*x2+ 10*x2,
-6*x1^4*x3^2-3*x1^2*x3^2*x2^2-3*x3^2*x2^4-x1^2*x3^2+ 2*x3^2*x2^2+ 28*x1^2*x3+ 7*x3*x2^2+x3^2-11*x3+ 10]
end

function ojika4_d1R2_d2R5(;ground=QQ)
    R, (x1, x2, x3) = PolynomialRing(ground, ["x$i" for i in 1:3])

    fs = [
    x1^3*x3+x1*x3*x2^2-x1*x3+x1,
-2*x1^2*x3*x2-x3*x2^3-x3*x2+ 10*x2,
-6*x1^4*x3^2-3*x1^2*x3^2*x2^2-3*x3^2*x2^4-x1^2*x3^2+ 2*x3^2*x2^2+ 28*x1^2*x3+ 7*x3*x2^2+x3^2-11*x3+ 10
    ]
end

function generate_set(nvariables, exps, nterms, npolys, csz, rng,
                            ground, ordering)
    R, _ = PolynomialRing(
        ground,
        ["x$i" for i in 1:nvariables],
        ordering=ordering
    )

    return filter!(!iszero, [
        map_coefficients(
            c -> ground(data(c) % csz),
            rand(rng, R, exps, nterms)
        )
        for _ in 1:rand(rng, npolys)
    ])
end

function generate_set(nvariables, exps, nterms, npolys, csz, rng,
                        ground::AbstractAlgebra.Rationals, ordering)

    semiground = GF(2^31 - 1)
    R, _ = PolynomialRing(
        semiground,
        ["x$i" for i in 1:nvariables],
        ordering=ordering
    )

    csz = BigInt(csz)
    zzbase = BigInt(9223372036854775837)
    randzz() = rand(-zzbase:zzbase)

    return filter!(!iszero, [
        map_coefficients(
            c -> ground(mod(randzz(), csz), mod(randzz(), csz) + 1),
            rand(rng, R, exps, nterms)
        )
        for _ in 1:rand(rng, npolys)
    ])
end


text:
using Test
using BigCleaner: Pipeline, add_step!

@testset "Pipeline is working as expected" begin
    p = Pipeline("a", push!)
    @test p isa Pipeline
    @test p.table == "a"

    add_step!(p, pop!)
    @test p.steps == [push!, pop!]

end


text:
import MacroTools

function check_observations(choices::ChoiceMap, observations::ChoiceMap)
    for (key, value) in get_values_shallow(observations)
        !has_value(choices, key) && error("Check failed: observed choice at $key not found")
        choices[key] != value && error("Check failed: value of observed choice at $key changed")
    end
    for (key, submap) in get_submaps_shallow(observations)
        check_observations(get_submap(choices, key), submap)
    end
end

is_custom_primitive_kernel(::Function) = false
check_is_kernel(::Function) = false

"""
    @pkern function k(trace, ..; check=false, observations=EmptyChoiceMap())
        ..
        return trace
    end

Declare a Julia function as a primitive stationary kernel.

The first argument of the function should be a trace, and the return value of the function should be a trace.
There should be keyword arguments check and observations.
"""
macro pkern(ex)
    MacroTools.@capture(ex, function f_(args__) body_ end) || error("expected a function")
    escaped_args = map(esc, args)
    quote
        function $(esc(f))($(escaped_args...))
            $(esc(body))
        end
        Gen.is_custom_primitive_kernel(::typeof($(esc(f)))) = true
        Gen.check_is_kernel(::typeof($(esc(f)))) = true
    end
end

function expand_kern_ex(ex)
    ex = MacroTools.postwalk(MacroTools.unblock, ex)
    MacroTools.@capture(ex, function f_(T_, args__) body_ end) || error("expected kernel syntax: function my_kernel(trace, .) .. end")
    trace = T
    toplevel_args = args

    body = MacroTools.postwalk(body) do x
        if MacroTools.@capture(x, for idx_ in range_ body_ end)

            # for loops
            quote
                loop_range = $(esc(range))
                for $(esc(idx)) in loop_range
                    $body
                end
                check && (loop_range != $(esc(range))) && error("Check failed in loop")
            end

        elseif MacroTools.@capture(x, if cond_ body_ end)

            # if .. end
            quote
                cond = $(esc(cond))
                if cond
                    $body
                end
                check && (cond != $(esc(cond))) && error("Check failed in if-end")
            end

        elseif MacroTools.@capture(x, let var_ = rhs_; body_ end)

            # let
            quote
                rhs = $(esc(rhs))
                let $(esc(var)) = rhs
                    $body
                end
                check && (rhs != $(esc(rhs))) && error("Check failed in let")
            end

        elseif MacroTools.@capture(x, let idx_ ~ dist_(args__); body_ end)

            # mixture
            quote
                dist = $(esc(dist))
                args = ($(map(esc, args)...),)
                let $(esc(idx)) = dist($(map(esc, args)...))
                    $body
                end
                check && (dist != $(esc(dist))) && error("Check failed in mixture (distribution)")
                check && (args != ($(map(esc, args)...),)) && error("Check failed in mixture (arguments)")
            end

        elseif MacroTools.@capture(x, T_ ~ k_(T_, args__))

            # applying a kernel
            quote
                check && Gen.check_is_kernel($(esc(k)))
                $(esc(T)) = $(esc(k))(
                    $(esc(T)), $(map(esc, args)...),
                    check=check, observations=observations)[1]
            end

        else

            # leave it as is
            x
        end
    end

    ex = quote
        function $(esc(f))(
                $(esc(trace))::Trace, $(toplevel_args...);
                check=false, observations=EmptyChoiceMap())
            $body
            check && check_observations(get_choices($(esc(trace))), observations)
            metadata = nothing
            ($(esc(trace)), metadata)
        end
        Gen.check_is_kernel(::typeof($(esc(f)))) = true
    end

    ex = MacroTools.postwalk(MacroTools.unblock, ex)

    ex, f
end

"""
    k2 = reversal(k1)

Return the reversal kernel for a given kernel.
"""
function reversal(f)
    check_is_kernel(f)
    error("Reversal for kernel $f is not defined")
end

"""
    @rkern k1 : k2

Declare that two primitive stationary kernels are reversals of one another.

The two kernels must have the same argument type signatures.
"""
macro rkern(ex)
    MacroTools.@capture(ex, k_ : l_) || error("expected a pair of functions")
    quote
        Gen.is_custom_primitive_kernel($(esc(k))) || error("first function is not a custom primitive kernel")
        Gen.is_custom_primitive_kernel($(esc(l))) || error("second function is not a custom primitive kernel")
        Gen.reversal(::typeof($(esc(k)))) = $(esc(l))
        Gen.reversal(::typeof($(esc(l)))) = $(esc(k))
    end
end

function reversal_ex(ex)

    ex = MacroTools.postwalk(MacroTools.unblock, ex)
    MacroTools.@capture(ex, function f_(args__) body_ end) || error("expected a function")
    toplevel_args = args

    # modify the body
    body = MacroTools.postwalk(body) do x
        x = if MacroTools.@capture(x, for idx_ in range_ body_ end)

            # for loops - reverse the order of loop indices
            quote
                for $idx in reverse($range)
                    $body
                end
            end

        elseif MacroTools.@capture(x, T_ ~ k_(T_, args__))

            # applying a kernel - apply the reverse kernel
            quote
                $T ~ (Gen.reversal($k))($T, $(args...))
            end

        elseif isa(x, Expr) && x.head == :block

            # a sequence of things -- reverse the order
            Expr(:block, reverse(x.args)...)

        else

            # leave it as is
            x
        end
    end

    # change the name
    rev = gensym("reverse_kernel")
    ex = quote
        function $rev($(toplevel_args...))
            $body
        end
    end

    ex = MacroTools.postwalk(MacroTools.unblock, ex)
end

"""
    @kern function k(trace, ..)
        ..
    end

Construct a composite MCMC kernel.

The resulting object is a Julia function that is annotated as a composite MCMC kernel, and can be called as a Julia function or applied within other composite kernels.
"""
macro kern(ex)
    kern_ex, kern = expand_kern_ex(ex)
    rev_kern_ex, rev_kern = expand_kern_ex(reversal_ex(ex))
    quote
        # define forward kerel
        $kern_ex

        # define reversal kernel
        $rev_kern_ex

        # bind the reversals for both
        Gen.reversal(::typeof($(esc(kern)))) = $(esc(rev_kern))
        Gen.reversal(::typeof($(esc(rev_kern)))) = $(esc(kern))
    end
end

export @pkern, @rkern, @kern, reversal


text:
#################### Slice Simplex Sampler ####################

#################### Types and Constructors ####################

mutable struct SliceSimplexTune <: SamplerTune
  logf::Union{Function, Missing}
  scale::Float64

  SliceSimplexTune() = new()

  SliceSimplexTune(x::Vector, logf::Union{Function, Missing}; scale::Real=1.0) =
    new(logf, scale)
end

SliceSimplexTune(x::Vector; args...) =
  SliceSimplexTune(x, missing; args...)

const SliceSimplexVariate = SamplerVariate{SliceSimplexTune}

function validate(v::SliceSimplexVariate)
  0 < v.tune.scale <= 1 || throw(ArgumentError("scale is not in (0, 1]"))
  validatesimplex(v)
end

#################### Sampler Constructor ####################

function SliceSimplex(params::ElementOrVector{Symbol}; args...)
  params = asvec(params)
  samplerfx = function(model::Model, block::Integer)
    s = model.samplers[block]
    local node, x
    for key in params
      node = model[key]
      x = unlist(node)

      sim = function(inds::AbstractRange, logf::Function)
        v = SamplerVariate(x[inds], s, model.iter; args...)
        sample!(v, logf)
      end

      logf = function(d::MultivariateDistribution, v::AbstractVector,
                      inds::AbstractRange)
        x[inds] = v
        relist!(model, x, key)
        logpdf(d, v) + logpdf(model, node.targets)
      end

      SliceSimplex_sub!(node.distr, sim, logf)
    end
    nothing
  end
  Sampler(params, samplerfx, SliceSimplexTune())
end

function SliceSimplex_sub!(d::MultivariateDistribution, sim::Function,
                           logf::Function)
  inds = 1:length(d)
  sim(inds, v -> logf(d, v, inds))
end

function SliceSimplex_sub!(D::Array{MultivariateDistribution}, sim::Function,
                           logf::Function)
  inds = 0:0
  for i in 1:length(D)
    d = D[i]
    inds = last(inds) .+ (1:length(d))
    sim(inds, v -> logf(d, v, inds))
  end
end

function SliceSimplex_sub!(d, sim::Function, logf::Function)
  throw(ArgumentError("unsupported distribution structure $(typeof(d))"))
end

#################### Sampling Functions ####################

sample!(v::SliceSimplexVariate) = sample!(v, v.tune.logf)

function sample!(v::SliceSimplexVariate, logf::Function)
  p0 = logf(v.value) + log(rand())
  d = Dirichlet(fill!(similar(v), 1))

  vertices = makefirstsimplex(v, v.tune.scale)
  vb = vertices \ v
  xb = rand(d)
  x = vertices * xb
  while any(x .< 0.0) || any(x .> 1.0) || logf(x) < p0
    vertices = shrinksimplex(vb, xb, v, x, vertices)
    vb = vertices \ v
    xb = rand(d)
    x = vertices * xb
  end
  v[:] = x

  v
end

function makefirstsimplex(x::AbstractVector{Float64}, scale::Real)
  vertices = Matrix{Float64}(I, length(x), length(x))
  vertices[:, 2:end] += (1.0 - scale) * (vertices[:, 1] .- vertices[:, 2:end])
  vertices .+ x .- vertices * rand(Dirichlet(fill!(similar(x), 1)))
end

function shrinksimplex(bx::AbstractVector{Float64}, bc::AbstractVector{Float64},
                       cx::AbstractVector{Float64}, cc::AbstractVector{Float64},
                       vertices::AbstractMatrix{Float64})
  for i in findall(bc .< bx)
    inds = [1:(i - 1); (i + 1):size(vertices, 2)]
    vertices[:, inds] += bc[i] * (vertices[:, i] .- vertices[:, inds])
    bc = vertices \ cc
  end
  vertices
end


text:
using Statistics
using JLD2
using Printf
using GLMakie

using Oceananigans.Utils: prettytime

file = jldopen("full_cubed_sphere_gravity_waves.jld2")

iterations = parse.(Int, keys(file["timeseries/t"]))

# Azimuthal spherical coordinates:
# Convert to λ ∈ [0°, 360°]
# Convert to φ ∈ [0°, 180°] (0° at north pole)
geographic2x(λ, φ, r=1) = r * cosd(λ + 180) * sind(90 - φ)
geographic2y(λ, φ, r=1) = r * sind(λ + 180) * sind(90 - φ)
geographic2z(λ, φ, r=1) = r * cosd(90 - φ)

Nf = file["grid/faces"] |> keys |> length
Nx = file["grid/faces/1/Nx"]
Ny = file["grid/faces/1/Ny"]
Nz = file["grid/faces/1/Nz"]

#####
##### Plot η
#####

size_2d = (Nx, Ny * Nf)

flatten_cubed_sphere(field, size) = reshape(cat(field..., dims=4), size)

λ = flatten_cubed_sphere((file["grid/faces/$n/λᶜᶜᵃ"][2:33, 2:33] for n in 1:6), size_2d)
φ = flatten_cubed_sphere((file["grid/faces/$n/φᶜᶜᵃ"][2:33, 2:33] for n in 1:6), size_2d)

x = geographic2x.(λ, φ)
y = geographic2y.(λ, φ)
z = geographic2z.(λ, φ)

iter = Node(0)

plot_title = @lift @sprintf("Surface gravity waves on a cubed sphere: time = %s", prettytime(file["timeseries/t/" * string($iter)]))

η = @lift flatten_cubed_sphere(file["timeseries/η/" * string($iter)], size_2d)

fig = Figure(resolution = (1920, 1080))
ax = fig[1, 1] = LScene(fig) # make plot area wider
wireframe!(ax, Sphere(Point3f0(0), 0.99f0), show_axis=false)
sf = surface!(ax, x, y, z, color=η, colormap=:balance, colorrange=(-0.01, 0.01))

rotate_cam!(ax.scene, (3π/4, 0, 0))
zoom!(ax.scene, (0, 0, 0), 8, false)

cb1 = fig[1, 2] = Colorbar(fig, sf, label="sea surface height η′ (m)", width=30)

supertitle = fig[0, :] = Label(fig, plot_title, textsize=50)

record(fig, "surface_gravity_waves_on_a_cubed_sphere.mp4", iterations, framerate=15) do i
    @info "Animating iteration $i/$(iterations[end])..."
    iter[] = i
end

#####
##### Plot u
#####

# size_2d = (Nx+1, Ny * Nf)

# λ = flatten_cubed_sphere((file["grid/faces/$n/λᶠᶠᵃ"][2:34, 2:33] for n in 1:6), size_2d)
# φ = flatten_cubed_sphere((file["grid/faces/$n/φᶠᶠᵃ"][2:34, 2:33] for n in 1:6), size_2d)

# x = geographic2x.(λ, φ)
# y = geographic2y.(λ, φ)
# z = geographic2z.(λ, φ)

# iter = Node(0)

# plot_title = @lift @sprintf("Surface gravity waves on a cubed sphere: time = %s", prettytime(file["timeseries/t/" * string($iter)]))

# u = @lift flatten_cubed_sphere(file["timeseries/u/" * string($iter)])

# fig = Figure(resolution = (1920, 1080))
# ax = fig[1, 1] = LScene(fig) # make plot area wider
# wireframe!(ax, Sphere(Point3f0(0), 0.99f0), show_axis=false)
# sf = surface!(ax, x, y, z, color=u, colormap=:balance, colorrange=(-1e-4, 1e-4))

# rotate_cam!(ax.scene, (3π/4, 0, 0))
# zoom!(ax.scene, (0, 0, 0), 8, false)

# cb1 = fig[1, 2] = Colorbar(fig, sf, label="u-velocity (m/s)", width=30)

# supertitle = fig[0, :] = Label(fig, plot_title, textsize=50)

# record(fig, "surface_gravity_waves_on_a_cubed_sphere_u.mp4", iterations[1:20], framerate=15) do i
#     @info "Animating iteration $i/$(iterations[end])..."
#     iter[] = i
# end

#####
##### Plot v
#####

# size_2d = (Nx, (Ny+1) * Nf)

# λ = flatten_cubed_sphere((file["grid/faces/$n/λᶠᶠᵃ"][2:33, 2:34] for n in 1:6), size_2d)
# φ = flatten_cubed_sphere((file["grid/faces/$n/φᶠᶠᵃ"][2:33, 2:34] for n in 1:6), size_2d)

# x = geographic2x.(λ, φ)
# y = geographic2y.(λ, φ)
# z = geographic2z.(λ, φ)

# iter = Node(0)

# plot_title = @lift @sprintf("Surface gravity waves on a cubed sphere: time = %s", prettytime(file["timeseries/t/" * string($iter)]))

# v = @lift flatten_cubed_sphere(file["timeseries/v/" * string($iter)])

# fig = Figure(resolution = (1920, 1080))
# ax = fig[1, 1] = LScene(fig) # make plot area wider
# wireframe!(ax, Sphere(Point3f0(0), 0.99f0), show_axis=false)
# sf = surface!(ax, x, y, z, color=v, colormap=:balance, colorrange=(-1e-4, 1e-4))

# rotate_cam!(ax.scene, (3π/4, 0, 0))
# zoom!(ax.scene, (0, 0, 0), 8, false)

# cb1 = fig[1, 2] = Colorbar(fig, sf, label="v-velocity (m/s)", width=30)

# supertitle = fig[0, :] = Label(fig, plot_title, textsize=50)

# record(fig, "surface_gravity_waves_on_a_cubed_sphere_v.mp4", iterations[1:20], framerate=15) do i
#     @info "Animating iteration $i/$(iterations[end])..."
#     iter[] = i
# end

close(file)


text:
@testset "Basic tests" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    MOI.set(optimizer, MOI.Silent(), true)

    x = MOI.add_variables(optimizer,2)
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z = MOI.VariableIndex(4)
    cz = MOI.ConstraintIndex{MOI.SingleVariable, POI.Parameter}(4)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    @test_throws ErrorException("Cannot constrain a parameter") MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.EqualTo(0.0))

    @test_throws ErrorException("Variable not in the model") MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.GreaterThan(0.0))

    cons1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], y]), 0.0)

    c1 = MOI.add_constraint(optimizer, cons1, MOI.EqualTo(2.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], y]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 2

    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 2

    @test_throws ErrorException("Variable not in the model") MOI.get(optimizer, MOI.VariablePrimal(), z)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))

    @test_throws ErrorException("Parameter not in the model") MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 2

    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 1

    new_obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), new_obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 1

    @test MOI.get(optimizer, MOI.SolverName()) == "ParametricOptimizer with GLPK attached"

    @test MOI.supports(optimizer, MOI.VariableName(), MOI.VariableIndex)
    @test MOI.supports(optimizer, MOI.ConstraintName(), MOI.ConstraintIndex)
    @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    @test MOI.get(optimizer, MOI.VariableName(), x[1]) == ""
    @test MOI.get(optimizer, MOI.ConstraintName(), c1) == ""

end

@testset "Quadratic objective parameter x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [1.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(1))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(1))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], y, z))

    objective_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, x),
                            quad_terms,
                            0.0
                        )

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 1.0, atol = ATOL)
    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 0

    @test_throws ErrorException("Cannot calculate the dual of a multiplicative parameter") MOI.get(optimizer, MOI.ConstraintDual(), cy)
    @test_throws ErrorException("Cannot calculate the dual of a multiplicative parameter") MOI.get(optimizer, MOI.ConstraintDual(), cz)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 2.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(3.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 6.0, atol = ATOL)

    MOI.set(optimizer, POI.ParameterValue(), y, 5)
    MOI.set(optimizer, POI.ParameterValue(), z, 5.0)
    @test_throws ErrorException MOI.set(optimizer, POI.ParameterValue(), MOI.VariableIndex(10872368175), 5.0)
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 25.0, atol = ATOL)

end

@testset "Quadratic objective parameter in affine part" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(1))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(1))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], x[2], x[2]))

    objective_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, [y,z]),
                            quad_terms,
                            0.0
                        )

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 0, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 2.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), 1.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 5.0, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 2.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), 1.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(3.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(5.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(5.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 15.0, atol = ATOL)

end

@testset "Quadratic constraint parameter x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], y, y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], y, z))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], z, z))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, x),
                            quad_terms,
                            0.0
                        )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 50.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 20.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 10.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 42.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 36.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(-1.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 45.0, atol = ATOL)

end

@testset "Quadratic constraint no parameters" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [0.0, 0.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(1.0))
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.LessThan(5.0))
    end

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, x),
                            [MOI.ScalarQuadraticTerm(A[1,2], x[1], x[2])],
                            0.0
                        )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(9.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 11.8, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 5.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 1.8, atol = ATOL)

end

@testset "Quadratic constraint parameter x variable" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], y, y))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, x),
                            quad_terms,
                            0.0
                        )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 30.25, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 0.5, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 29.25, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 22.0, atol = ATOL)
end

@testset "Quadratic constraint variable x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], y, y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], y, x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], x[1], x[1]))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, x),
                            quad_terms,
                            0.0
                        )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 30.25, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 0.5, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 29.25, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 22.0, atol = ATOL)
end

@testset "Quadratic constraint variable x variable + parameter in affine part" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], x[2], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, [x[1], y]),
                            quad_terms,
                            0.0
                        )

    cons_index = MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 9.0664, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4.3665, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 1/3, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), MOI.get(optimizer, MOI.ConstraintDual(), cons_index), atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 8.6609, atol = ATOL)
end

@testset "Quadratic constraint variable x variable + parameter in affine part - variation to assess duals" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 2.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], x[2], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, [x[1], y]),
                            quad_terms,
                            0.0
                        )

    cons_index = MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 9.0664, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4.3665, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 1/3, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 2*MOI.get(optimizer, MOI.ConstraintDual(), cons_index), atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 8.2376, atol = ATOL)
end

@testset "Special cases of getters" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,1], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1,2], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2,2], x[2], y))

    constraint_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(a, [x[1], y]),
                            quad_terms,
                            0.0
                        )

    cons_index = MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), [MOI.ScalarQuadraticTerm(A[2,2], x[2], y)], 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(optimizer, MOI.ObjectiveFunctionType()) == MOI.ScalarQuadraticFunction{Float64}
    @test MOI.get(optimizer, MOI.NumberOfVariables()) == 3

end


text:
# Autogenerated wrapper script for oneTBB_jll for x86_64-apple-darwin
export libtbb, libtbbmalloc, libtbbmalloc_proxy

JLLWrappers.@generate_wrapper_header("oneTBB")
JLLWrappers.@declare_library_product(libtbb, "@rpath/libtbb.12.dylib")
JLLWrappers.@declare_library_product(libtbbmalloc, "@rpath/libtbbmalloc.2.dylib")
JLLWrappers.@declare_library_product(libtbbmalloc_proxy, "@rpath/libtbbmalloc_proxy.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtbb,
        "lib/libtbb.12.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc,
        "lib/libtbbmalloc.2.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc_proxy,
        "lib/libtbbmalloc_proxy.2.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@everywhere using DistributedArrays
include("../utils.jl")

# Get the cell id for a given x. Note cells are labeled 1,2,3,...
@everywhere function getCellID2d(
    p::Array{Float64,1},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    bdryConstraints::Bool = false,
  )
  #clamp
  p1::Float64 = min(nCols-1e-12, max(0.0, p[1]))
  p2::Float64 = min(nRows-1e-12, max(0.0, p[2]))

  cellId::Int64 = round(Int, floor(Int64, p1/incX)+ nRows/incY*floor(Int64, p2/incY))
  cellId *= 4

  #out of bounds case
  if !bdryConstraints
    #left
    if p[1] ≤ 0.0
      if (p[2]≤0.0 && p[2]/incY<p[1]/incX)
        cellId += 1
      elseif (p[2]≥nRows && (p[2]-nRows)/incY>-p[1]/incX)
        cellId += 3
      else
        cellId += 4
      end
      return cellId
    end

    #right
    if p[1] ≥ nCols
      if (p[2]≤0.0 && -p[2]/incY>(p[1]-nCols)/incX)
        cellId += 1
      elseif (p[2]≥nRows && (p[2]-nRows)/incY>(p[1]-nCols)/incX)
        cellId += 3
      else
        cellId += 2
      end
      return cellId
    end

    #up
    if p[2]≤0.0
      cellId += 1
      return cellId
    end

    #down
    if p[2]≥nRows
      cellId += 3
      return cellId
    end
  end

  #inbounds
  x::Float64 = mod(p1, incX)/incX
  y::Float64 = mod(p2, incY)/incY
  if (x<y)
    if (1-x<y)
      cellId += 3 #down
    else
      cellId += 4 #left
    end
  elseif (1-x<y)
    cellId += 2 #right
  else
    cellId += 1 #top
  end
  cellId
end

# does Ap+B
@everywhere function getVelocity2d(
    p::Array{Float64,1},
    As::Array{Float64,3},
    cellId::Int64,
    v::Array{Float64,1} = Array(Float64,2),
  )
  v[1] = As[cellId,1,1]*p[1] + As[cellId,1,2]*p[2] + As[cellId,1,3]
  v[2] = As[cellId,2,1]*p[1] + As[cellId,2,2]*p[2] + As[cellId,2,3]
  v
end

@everywhere function calculateNumericalPhi(
    p::Array{Float64,1},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    nSteps::Int64,
    bdryConstraints::Bool = false,
    midpoint::Array{Float64,1} = Array(Float64,2),
    velocity::Array{Float64,1} = Array(Float64,2),
  )

  yn::Array{Float64,1} = p
  Δt::Float64 = t/nSteps
  cellId::Int64 = 1
  for i::Int64 = 1:nSteps
    cellId = getCellID2d(yn, nCols, nRows, incX, incY, bdryConstraints)
    getVelocity2d(yn, As, cellId, velocity)
    midpoint[1] = yn[1]+Δt/2*velocity[1]
    midpoint[2] = yn[2]+Δt/2*velocity[2]
    yn[1] += Δt*velocity[1]
    yn[2] += Δt*velocity[2]
  end
  yn
end

# Result is computed in pts_fwd
@everywhere function solverODE(
    pts_src::Array{Float64,2},
    pts_fwd::Array{Float64,2},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    nSteps::Int64,
    bdryConstraints::Bool = false
  )

  nPts::Int64 = size(pts_src,1)
  for i::Int64=1:nPts
    pts_fwd[i,:] = calculateNumericalPhi(squeeze(pts_src[i,:],1), t, As, nCols, nRows, incX, incY, nSteps, bdryConstraints)
  end
  pts_fwd
end

# Result is computed in pts_fwd
@everywhere function solverODEParallel(
    pts_src::Array{Float64,2},
    pts_fwd::DArray{Float64,2},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    nSteps::Int64,
    bdryConstraints::Bool = false
  )

  nPts::Int64 = size(pts_src,1)
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  # spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(pts_fwd)[1];
      localArray = localpart(pts_fwd);
      tic();
      for i::Int64 = 1:length(indicesOnThread);
        localArray[i,:] = calculateNumericalPhi(squeeze(pts_src[indicesOnThread[i],:],1), t, As, nCols, nRows, incX, incY, nSteps, bdryConstraints);
      end;
      toc();
      ))
  end

  # wait for the workers to finish
  pmap(wait, refs)
  pts_fwd
end

# expm seems faster than this, Assume 2d matrix is A[i,1:2,1:2], Stores result in expT[i,1:2,1:2]
@everywhere function expm2x2(A::Array{Float64,3}, expT::Array{Float64,3}, i::Int64)
  a::Float64=A[i,1,1]
  b::Float64=A[i,1,2]
  c::Float64=A[i,2,1]
  d::Float64=A[i,2,2]

  Δtmp::Float64 = (a-d)^2 + 4*b*c
  avgADexp::Float64 = exp((a+d)/2)
  Δ::Float64
  sinhΔoverΔ::Float64

  if Δtmp == 0
    expT[i,1,1] = (1 + (a-d)/2) * avgADexp
    expT[i,1,2] = b * avgADexp
    expT[i,2,1] = c * avgADexp
    expT[i,2,2] = (1 - (a-d)/2) * avgADexp

  elseif Δtmp > 0
    Δ = sqrt(Δtmp) / 2
    coshΔ::Float64 = cosh(Δ)
    sinhΔ::Float64 = sinh(Δ)
    sinhΔoverΔ = sinhΔ / Δ

    expT[i,1,1] = (coshΔ + (a-d)/2 * sinhΔoverΔ) * avgADexp
    expT[i,1,2] = b * sinhΔoverΔ  * avgADexp
    expT[i,2,1] = c * sinhΔoverΔ  * avgADexp
    expT[i,2,2] = (coshΔ - (a-d)/2 * sinhΔoverΔ) * avgADexp
  else
    Δ = sqrt(-Δtmp) / 2
    cosΔ::Float64 = cos(Δ)
    sinΔ::Float64 = sin(Δ)
    sinΔoverΔ = sinΔ / Δ

    expT[i,1,1] = (cosΔ + (a-d)/2 * sinΔoverΔ) * avgADexp
    expT[i,1,2] = b * sinΔoverΔ * avgADexp
    expT[i,2,1] = c * sinΔoverΔ * avgADexp
    expT[i,2,2] = (cosΔ - (a-d)/2 * sinΔoverΔ) * avgADexp
  end
end

# Stores result in expT[i,:,:], Assume 2d matrix is A[i,:,:]
@everywhere function expm3x3(A::Array{Float64,3}, expT::Array{Float64,3}, i::Int64)
  det2x2::Float64 = A[i,1,1]*A[i,2,2] - A[i,1,2]*A[i,2,1]
  if det2x2 != 0.0
    expm2x2(A, expT, i)

    a = A[i,2,2]/det2x2*A[i,1,3] - A[i,1,2]/det2x2*A[i,2,3]
    b = -A[i,2,1]/det2x2*A[i,1,3] + A[i,1,1]/det2x2*A[i,2,3]
    expT[i,1,3] = (expT[i,1,1]-1)*a + expT[i,1,2]*b
    expT[i,2,3] = expT[i,2,1]*a + (expT[i,2,2]-1)*b

    #last row 0,0,1
    expT[i,3,1] = expT[i,3,2] = 0
    expT[i,3,3] = 1
  else
    Ã::Array{Float64,2} = Array(Float64,3,3)
    Ã[1:2,:] = squeeze(A[i,:,:],1)
    Ã[3,1] = Ã[3,2] = Ã[3,3] = 0
    expT[i,:,:] = expm(Ã)
  end
end

@everywhere function expm3x3(A::Array{Float64,2})
  det2x2::Float64 = A[1,1]*A[2,2] - A[1,2]*A[2,1]
  expT::Array{Float64,2} = Array(Float64,3,3)
  if det2x2 != 0.0
    expT[1:2,1:2] = expm2x2(A[1:2,1:2])

    a = A[2,2]/det2x2*A[1,3] - A[1,2]/det2x2*A[2,3]
    b = -A[2,1]/det2x2*A[1,3] + A[1,1]/det2x2*A[2,3]
    expT[1,3] = (expT[1,1]-1)*a + expT[1,2]*b
    expT[2,3] = expT[2,1]*a + (expT[2,2]-1)*b

    #last row 0,0,1
    expT[3,1] = expT[3,2] = 0
    expT[3,3] = 1
  else
    Ã::Array{Float64,2} = Array(Float64,3,3)
    Ã[1:2,:] = A
    Ã[3,1] = Ã[3,2] = Ã[3,3] = 0
    expT[:,:] = expm(Ã)
  end
  expT
end

@everywhere function solver2dPhi(
    p::Array{Float64,1},
    t::Float64,
    As::Array{Float64,3},
    expT::Array{Float64,3},
    expTHalf::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    bdryConstraints::Bool = false,
    midpoint::Array{Float64,1} = Array(Float64,2),
    velocity::Array{Float64,1} = Array(Float64,2),
    pTemp::Array{Float64,1} = Array(Float64,2),
  )
  pPrev::Array{Float64,1} = p
  cellId::Int64 = getCellID2d(pPrev, nCols, nRows, incX, incY, bdryConstraints)
  for j::Int64 = 1:Nsteps
    getVelocity2d(pPrev, expT, cellId, pTemp) # expT[cellId] * p
    if cellId == getCellID2d(pTemp, nCols, nRows, incX, incY, bdryConstraints)
      pPrev[1] = pTemp[1]
      pPrev[2] = pTemp[2]
    else
      getVelocity2d(pPrev, expTHalf, cellId, pTemp)
      if cellId == getCellID2d(pTemp, nCols, nRows, incX, incY, bdryConstraints)
        pPrev[1] = pTemp[1]
        pPrev[2] = pTemp[2]
        pPrev = calculateNumericalPhi(pPrev, Δt/2, As, nCols, nRows, incX, incY, round(Int64, nSteps/2), bdryConstraints, midpoint, velocity)
        cellId = getCellID2d(pPrev, nCols, nRows, incX, incY, bdryConstraints)

      else
        pPrev = calculateNumericalPhi(pPrev, Δt, As, nCols, nRows, incX, incY, nSteps, bdryConstraints, midpoint, velocity)
        cellId = getCellID2d(pPrev, nCols, nRows, incX, incY, bdryConstraints)
      end
    end
  end
  pPrev
end

# Stores result in pts_fwd
@everywhere function solver2d(
    pts_src::Array{Float64,2},
    pts_fwd::Array{Float64,2},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    bdryConstraints::Bool = false,
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
  )
  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs, expT, i)
  end

  expTHalf = Array(Float64, size(As,1), 3, 3)
  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs/2, expTHalf, i)
  end

  nPts::Int64 = size(pts_src,1)
  midpoint::Array{Float64,1} = Array(Float64,2)
  velocity::Array{Float64,1} = Array(Float64,2)
  currPoint::Array{Float64,1} = Array(Float64,2)
  pTemp::Array{Float64,1} = Array(Float64,2)
  for i::Int64=1:nPts
    currPoint[1] = pts_src[i,1]
    currPoint[2] = pts_src[i,2]
    pts_fwd[i,:] = solver2dPhi(currPoint, t, As, expT, expTHalf, nCols, nRows, incX, incY, Nsteps, Δt, nSteps, bdryConstraints, midpoint, velocity, pTemp)
  end
  pts_fwd
end

# Stores result in pts_fwd (note: Distributed Array)
@everywhere function solver2dParallel(
    pts_src::Array{Float64,2},
    pts_fwd::DArray{Float64,2},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    bdryConstraints::Bool = false,
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
  )

  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs, expT, i)
  end

  expTHalf = Array(Float64, size(As,1), 3, 3)
  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs/2, expTHalf, i)
  end

  nPts::Int64 = size(pts_src,1)
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  #spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(pts_fwd)[1];
      localArray = localpart(pts_fwd);
      midpoint::Array{Float64,1} = Array(Float64,2);
      velocity::Array{Float64,1} = Array(Float64,2);
      currPoint::Array{Float64,1} = Array(Float64,2);
      pTemp::Array{Float64,1} = Array(Float64,2);
      for i::Int64 = 1:length(indicesOnThread);
        currPoint[1] = pts_src[indicesOnThread[i],1]
        currPoint[2] = pts_src[indicesOnThread[i],2]
        localArray[i,:] = solver2dPhi(currPoint, t, As, expT, expTHalf, nCols, nRows, incX, incY, Nsteps, Δt, nSteps, bdryConstraints, midpoint, velocity, pTemp);
      end;
      ))
  end

  # wait for the workers to finish
  map(wait, refs)
  pts_fwd
end

# This function retuns the range indexes assigned to this worker
@everywhere function myrange(q::SharedArray)
    idx = indexpids(q)
    if idx == 0
        # This worker is not assigned a piece
        return 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in linspace(0,size(q,1),nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

# Stores result in pts_fwd (note: Shared Array)
@everywhere function solver2dParallelShared(
    pts_src::Array{Float64,2},
    pts_fwd::SharedArray{Float64,2},
    t::Float64,
    As::Array{Float64,3},
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    bdryConstraints::Bool = false,
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
  )

  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs, expT, i)
  end

  expTHalf = Array(Float64, size(As,1), 3, 3)
  for i::Int64 = 1:size(As, 1)
    expm3x3(ΔtAs/2, expTHalf, i)
  end

  nPts::Int64 = size(pts_src,1)
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  # spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      midpoint::Array{Float64,1} = Array(Float64,2);
      velocity::Array{Float64,1} = Array(Float64,2);
      currPoint::Array{Float64,1} = Array(Float64,2);
      pTemp::Array{Float64,1} = Array(Float64,2);
      for i::Int64 = myrange(pts_fwd);
        currPoint[1] = pts_src[i,1];
        currPoint[2] = pts_src[i,2];
        pts_fwd[i,:] = solver2dPhi(currPoint, t, As, expT, expTHalf, nCols, nRows, incX, incY, Nsteps, Δt, nSteps, bdryConstraints, midpoint, velocity, pTemp);
      end;
      ))
  end

  # wait for the workers to finish
  map(wait, refs)
  pts_fwd
end


text:
module RocketMaxOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: max()" begin

    println("Testing: operator max()")

    run_proxyshowcheck("Max", max())

    run_testset([
        (
            source      = from(1:42) |> max(),
            values      = @ts([ 42, c ]),
            source_type = Union{Int, Nothing}
        ),
        (
            source      = from(1:42) |> max(from = 100),
            values      = @ts([ 100, c ]),
            source_type = Union{Int, Nothing}
        ),
        (
            source      = from(1:42) |> max() |> some(),
            values      = @ts([ 42, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> max(),
            values      = @ts([ nothing, c ]),
            source_type = Union{Int, Nothing}
        ),
        (
            source      = faulted(String, "e") |> max(),
            values      = @ts(e("e")),
            source_type = Union{Nothing, String}
        ),
        (
            source = never() |> max(),
            values = @ts()
        )
    ])

end

end