text:
include("constants.jl")
include("utils.jl")

"""
Convert a byte array into an integer array. The number of bytes forming an integer
is defined by num
in_bytes: the input bytes
num: the number of bytes per int
"""
function ints(in_bytes::Array{UInt8,1}, num)
    out_arr = Int[]
    if num == 1
        out_arr = htonofarray(reinterpret(Int8,in_bytes))
    elseif num == 2
        out_arr = htonofarray(reinterpret(Int16,in_bytes))
    elseif num == 4
        out_arr = htonofarray(reinterpret(Int32,in_bytes))
    else
        throw(ArgumentError("Cannot convert bytes to integer array. Number of bytes parsed is wrong"))
    end
    return out_arr
end

"""
Convert floating points to integers using a multiplier.
in_floats: the input floats
multiplier: the multiplier to be used for conversion. Corresponds to the precisison.
"""
function ints(in_floats::Array{Float32,1}, multiplier)
    return Int32[Int32(round(x * multiplier)) for x in in_floats]
end

"""
Convert an array of chars to an array of ints.
"""
function ints(in_chars::Array{Char,1})
    return Int32[Int32(x) for x in in_chars]
end

"""
Convert an integer array into a byte arrays.
"""
function bytes(in_ints::Array{T,1}) where {T<:Integer}
    out_array=UInt8[]
    return append!(out_array,reinterpret(UInt8,ntohofarray(in_ints)))
end

"""
Convert integers to floats by division.
in_ints: the integer array
divider: the divider
"""
function floats(in_ints::Array{T,1}, divider) where {T<:Integer}
    return Float32[x/divider for x in in_ints]
end

"""
Convert integers to chars.
"""
function chars(in_ints::Array{T,1}) where {T<:Integer}
    return Char[Char(x) for x in in_ints]
end

"""
Convert a list of bytes to a list of strings. Each string is of length mmtf.CHAIN_LEN
"""
function decodechainlist(in_bytes::Array{UInt8,1})
    tot_strings = div.(length(in_bytes),CHAIN_LEN)
    out_strings = String[]
    for i=1:tot_strings
        out_s = in_bytes[(i*CHAIN_LEN)-(CHAIN_LEN-1) : i * CHAIN_LEN]
        stripped_s = strip(String(out_s),NULL_BYTE)
        push!(out_strings,String(stripped_s))
    end
    return out_strings
end

"""
Convert a list of strings to a list of byte arrays.
"""
function encodechainlist(in_strings::Array{String,1})
    out_bytes=UInt8[]
    for in_s in in_strings
        append!(out_bytes, transcode(UInt8,in_s))
        for i=1:(CHAIN_LEN - length(in_s))
            append!(out_bytes, transcode(UInt8,"$(NULL_BYTE)"))
        end
    end
    return out_bytes
end

"""
Unpack an array of integers using recursive indexing.
"""
function recursiveindexdecode(int_array::Array{Int16,1})
    out_arr = Int32[]
    encoded_ind = 1
    while encoded_ind <= length(int_array)
        decoded_val = 0
        while int_array[encoded_ind]==MAX_SHORT || int_array[encoded_ind]==MIN_SHORT
            decoded_val += int_array[encoded_ind]
            encoded_ind+=1
            if int_array[encoded_ind]==0
                break
            end
        end
        decoded_val += int_array[encoded_ind]
        encoded_ind+=1
        push!(out_arr,Int32(decoded_val))
    end
    return out_arr
end

"""
Pack an integer array using recursive indexing.
"""
function recursiveindexencode(int_array::Array{Int32,1})
    out_arr = Int16[]
    for curr in int_array
        if curr >= 0
            while curr >= MAX_SHORT
                push!(out_arr,MAX_SHORT)
                curr -=  MAX_SHORT
            end
        else
            while curr <= MIN_SHORT
                push!(out_arr,MIN_SHORT)
                curr += Int(abs(MIN_SHORT))
            end
        end
        push!(out_arr, Int16(curr))
    end
    return out_arr
end


text:
module OpenAIGym

using PyCall
import PyCall: hasproperty
using Reexport
@reexport using MinimalRLCore
import MinimalRLCore: AbstractEnvironment
# @reexport using Reinforce
# import Reinforce:
#     MouseAction, MouseActionSet,
#     KeyboardAction, KeyboardActionSet

export
    GymEnv,
    render,
    close

struct DiscreteSet{T<:AbstractArray} <: AbstractSet{T}
    items::T
end
randtype(s::DiscreteSet) = eltype(s.items)
Base.rand(s::DiscreteSet, dims::Integer...) = rand(s.items, dims...)
Base.in(x, s::DiscreteSet) = x in s.items
Base.length(s::DiscreteSet) = length(s.items)
Base.getindex(s::DiscreteSet, i::Int) = s.items[i]
Base.:(==)(s1::DiscreteSet, s2::DiscreteSet) = s1.items == s2.items
Base.isempty(s::DiscreteSet) = Base.isempty(s.items)

# --------------------------------------------------------------

abstract type AbstractGymEnv <: AbstractEnvironment end

"A simple wrapper around the OpenAI gym environments to add to the Reinforce framework"
mutable struct GymEnv{T} <: AbstractGymEnv
    name::Symbol
    ver::Symbol
    pyenv::PyObject   # the python "env" object
    pystep::PyObject  # the python env.step function
    pyreset::PyObject # the python env.reset function
    pystate::PyObject # the state array object referenced by the PyArray state.o
    pystepres::PyObject # used to make stepping the env slightly more efficient
    info::PyObject    # store it as a PyObject for speed, since often unused
    state::T
    reward::Float64
    total_reward::Float64
    actions::AbstractSet
    done::Bool
    function GymEnv{T}(name, ver, pyenv, pystate, state, seed) where T
        env = new{T}(name, ver, pyenv, pyenv."step", pyenv."reset",
                     pystate, PyNULL(), PyNULL(), state)
        env.pyenv.seed(seed)
        MinimalRLCore.reset!(env)
        env
    end
end

use_pyarray_state(envname::Symbol) = !(envname ∈ (:Blackjack,))

function GymEnv(name::Symbol, ver::Symbol = :v0;
                stateT = ifelse(use_pyarray_state(name), PyArray, PyAny),
                seed=0)
    if PyCall.ispynull(pysoccer) && name ∈ (:Soccer, :SoccerEmptyGoal)
        copy!(pysoccer, pyimport("gym_soccer"))
    end

    GymEnv(name, ver, pygym.make("$name-$ver"), stateT, seed)
end

GymEnv(name::AbstractString; kwargs...) =
    GymEnv(Symbol.(split(name, '-', limit = 2))...; kwargs...)

function GymEnv(name::Symbol, ver::Symbol, pyenv, stateT, seed)
    pystate = pycall(pyenv."reset", PyObject)
    state = convert(stateT, pystate)
    T = typeof(state)
    GymEnv{T}(name, ver, pyenv, pystate, state, seed)
end

function Base.show(io::IO, env::GymEnv)
  println(io, "GymEnv $(env.name)-$(env.ver)")
  if hasproperty(env.pyenv, :class_name)
    println(io, "  $(env.pyenv.class_name())")
  end
  println(io, "  r  = $(env.reward)")
  print(  io, "  ∑r = $(env.total_reward)")
end

# --------------------------------------------------------------

"""
    close(env::AbstractGymEnv)
"""
Base.close(env::AbstractGymEnv) =
	!ispynull(env.pyenv) && env.pyenv.close()

# --------------------------------------------------------------

"""
    render(env::AbstractGymEnv; mode = :human)

# Arguments

- `mode`: `:human`, `:rgb_array`, `:ansi`
"""
render(env::AbstractGymEnv, args...; kwargs...) =
    pycall(env.pyenv.render, PyAny; kwargs...)

# --------------------------------------------------------------

function actionset(A::PyObject)
    if hasproperty(A, :n)
        # choose from n actions
        DiscreteSet(0:A.n-1)
    elseif hasproperty(A, :spaces)
        # a tuple of action sets
        sets = [actionset(a) for a in A.spaces]
        TupleSet(sets...)
    elseif hasproperty(A, :high)
        # continuous interval
        IntervalSet{Vector{Float64}}(A.low, A.high)
        # if A[:shape] == (1,)  # for now we only support 1-length vectors
        #     IntervalSet{Float64}(A[:low][1], A[:high][1])
        # else
        #     # @show A[:shape]
        #     lo,hi = A[:low], A[:high]
        #     # error("Unsupported shape for IntervalSet: $(A[:shape])")
        #     [IntervalSet{Float64}(lo[i], hi[i]) for i=1:length(lo)]
        # end
    elseif hasproperty(A, :actions)
        # Hardcoded
        TupleSet(DiscreteSet(A.actions))
    else
        @show A
        @show keys(A)
        error("Unknown actionset type: $A")
    end
end

function MinimalRLCore.get_actions(env::AbstractGymEnv, s′)
    actionset(env.pyenv.action_space)
end

pyaction(a::Vector) = Any[pyaction(ai) for ai=a]
pyaction(a) = a

"""
`reset!(env::GymEnv)` reset the environment
"""
function MinimalRLCore.reset!(env::GymEnv, args...)
    pycall!(env.pystate, env.pyreset, PyObject)
    convert_state!(env)
    env.reward = 0.0
    env.total_reward = 0.0
    env.actions = MinimalRLCore.get_actions(env, nothing)
    env.done = false
    return env.state
end

"""
    step!(env::GymEnv, a)

take a step in the enviroment
"""
function MinimalRLCore.environment_step!(env::GymEnv, a, args...)
    pyact = pyaction(a)
    pycall!(env.pystepres, env.pystep, PyObject, pyact)

    env.pystate, env.reward, env.done, env.info =
        convert(Tuple{PyObject, Float64, Bool, PyObject}, env.pystepres)

    convert_state!(env)

    env.total_reward += env.reward
    return (env.reward, env.state)
end

convert_state!(env::GymEnv{T}) where T =
    env.state = convert(T, env.pystate)

convert_state!(env::GymEnv{<:PyArray}) =
    env.state = PyArray(env.pystate)

# Reinforce.finished(env::GymEnv)     = env.done
# Reinforce.finished(env::GymEnv, s′) = env.done
MinimalRLCore.is_terminal(env::GymEnv) = env.done
MinimalRLCore.get_reward(env::GymEnv) = env.reward
MinimalRLCore.get_state(env::GymEnv) = env.state

# --------------------------------------------------------------

const pygym    = PyNULL()
const pysoccer = PyNULL()

function __init__()
    # the copy! puts the gym module into `pygym`, handling python ref-counting
    copy!(pygym, pyimport("gym"))
end

end # module

ext: jl


text:
using Documenter, MerCounting, Pkg

makedocs(
    modules = [MerCounting],
    format = Documenter.HTML(),
    sitename = "MerCounting.jl",
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ),
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Installation" => "man/installation.md",
            "Basic types and counting" => "man/basic_counting.md"
        ],
        "API" => [
            "MerCount" => "api/MerCount.md"
        ]
    ],

)

deploydocs(
    repo = "github.com/BioJulia/MerCounting.jl.git",
    push_preview = true,
    deps = nothing,
    make = nothing
)


text:
"""
    AbstractRasterSeries <: DimensionalData.AbstractDimensionalArray

Abstract supertype for high-level `DimensionalArray` that hold `RasterStacks`, `Rasters`,
or the paths they can be loaded from. `RasterSeries` are indexed with dimensions
as with a `AbstractRaster`. This is useful when you have multiple files containing
rasters or stacks of rasters spread over dimensions like time and elevation.

As much as possible, implementations should facilitate loading entire
directories and detecting the dimensions from metadata.

This allows syntax like for a series of stacks of arrays:

```julia
RasterSeries[Time(Near(DateTime(2001, 1))][:temp][Y(Between(70, 150)), X(Between(-20,20))] |> plot`
```

[`RasterSeries`](@ref) is the concrete implementation.
"""
abstract type AbstractRasterSeries{T,N,D,A} <: AbstractDimArray{T,N,D,A} end

# Interface methods ####################################################

DD.metadata(A::AbstractRasterSeries) = NoMetadata()
DD.name(A::AbstractRasterSeries) = NoName()
DD.label(A::AbstractRasterSeries) = ""

"""
    modify(f, series::AbstractRasterSeries)

Apply function `f` to the data of the child object.
If the child is an `AbstractRasterStack` the function will
be passed on to its child `AbstractRaster`s.

`f` must return an idenically sized array.

This method triggers a complete rebuild of all objects,
and disk based objects will be transferred to memory.

This is useful for swapping out array backend for an
entire series to `CuArray` from CUDA.jl to copy data to a GPU,
and potentially other types like `DAarray` from Distributed.jl.
"""
DD.modify(f, A::AbstractRasterSeries) = map(child -> modify(f, child), values(A))

"""
    RasterSeries <: AbstractRasterSeries

    RasterSeries(arrays::AbstractArray{<:AbstractRaster}, dims; kw...)
    RasterSeries(stacks::AbstractArray{<:AbstractRasterStack}, dims; kw...)
    RasterSeries(filepaths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
    RasterSeries(dirpath:::AbstractString, dims; ext, child, duplicate_first, kw...)

Concrete implementation of [`AbstractRasterSeries`](@ref).
`RasterSeries` are an array of `Raster` or `RasterStack`, along some dimension(s).

# Arguments

- `dims`: series dimension/s.

# Keywords

- `refdims`: existing reference dimension/s.
- `child`: constructor of child objects for use with filenames are passed in,
    can be `Raster` or `RasterStack`. Defaults to `Raster`.
- `duplicate_first::Bool`: wether to duplicate the dimensions and metadata of the
    first file with all other files. This can save load time with a large
    series where dimensions are essentially identical. `true` by default to improve
    load times. If you need exact metadata, set to `false`.
- `ext`: filename extension such as ".tiff" to find when only a directory path is passed in.
- `kw`: keywords passed to the child constructor [`Raster`](@ref) or [`RasterStack`](@ref)
    if only file names are passed in.
"""
struct RasterSeries{T,N,D,R,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
    data::A
    dims::D
    refdims::R
end
function RasterSeries(data::AbstractArray{<:Union{AbstractRasterStack,AbstractRaster}}, dims;
    refdims=()
)
    RasterSeries(data, DD.format(dims, data), refdims)
end
function RasterSeries(filenames::NamedTuple{K}, dims; kw...) where K
    RasterSeries(map((fns...) -> NamedTuple{K}(fns), values(filenames)...), dims; kw...)
end
function RasterSeries(filenames::AbstractArray{<:Union{AbstractString,NamedTuple}}, dims;
    refdims=(), duplicate_first=true, child=nothing, resize=nothing, kw...
)
    childtype = if isnothing(child)
        eltype(filenames) <: NamedTuple ? RasterStack : Raster
    else
        child
    end
    data = if duplicate_first
        # We assume all dims, metadata and missingvals are the same over the series
        # We just load the first object, and swap in the filenames of the others.
        data1 = if childtype <: AbstractRaster
            childtype(first(filenames); kw...)
        else
            childtype(first(filenames); resize, kw...)
        end
        map(filenames) do fn
            swap_filename(data1, fn)
        end
    else
        # Load everything separately
        if childtype <: AbstractRaster
            [childtype(fn; kw...) for fn in filenames]
        else
            [childtype(fn; resize, kw...) for fn in filenames]
        end
    end
    return RasterSeries(data, DD.format(dims, data); refdims)
end
function RasterSeries(dirpath::AbstractString, dims; ext=nothing, kw...)
    filepaths = filter_ext(dirpath, ext)
    RasterSeries(filepaths, dims; kw...)
end

@inline function DD.rebuild(
    A::RasterSeries, data, dims::Tuple, refdims=(), name=nothing, metadata=nothing,
)
    RasterSeries(data, dims, refdims)
end
@inline function DD.rebuild(
    A::RasterSeries;
    data=parent(A), dims=dims(A), refdims=refdims(A), name=nothing, metadata=nothing,
)
    RasterSeries(data, dims, refdims)
end

@deprecate series(args...; kw...) RasterSeries(args...; kw...)

# Swap in the filename/s of an object for another filename, wherever it is.
# This is used to use already loaded metadata of one file with another
# file that is similar or identical besides tha actual raster data.
swap_filename(x, filename) = rebuild(x, data=swap_filename(parent(x), filename))
swap_filename(x::NamedTuple, filenames::NamedTuple) = map(swap_filename, x, filenames)
swap_filename(x::FileStack, filename::AbstractString) = @set x.filename = filename
swap_filename(x::FileArray, filename::AbstractString) = @set x.filename = filename
function swap_filename(x::AbstractArray, filename::AbstractString)
    # The `FileArray` is wrapped, so use Flatten.jl to update it wherever it is
    ignore = Union{Dict,Set,Base.MultiplicativeInverses.SignedMultiplicativeInverse}
    Flatten.modify(x, FileArray, ignore) do fa
        @set fa.filename = filename
    end
end


text:
# -----------------------------------------------------------------------------
##### CYLINDER/SHELL #####-----------------------------------------------------
function shellregioncylinder(dataobject::PartDataType;
                            radius::Array{<:Real,1}=[0.,0.],
                            height::Real=0.,
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            direction::Symbol=:z,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    radius_in  = radius[1]
    radius_out = radius[2]
    if radius_in == 0. || radius_out == 0. || height == 0. || in(0., center)
        error("[Mera]: given radius, height or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift, height_shift = prep_cylindrical_shellranges(dataobject.info, center, radius_in, radius_out, height, range_unit, verbose)

    if inverse == false
        sub_data = filter(p-> sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2)
                                    >= ( radius_in_shift*boxlen )  &&

                              sqrt( (p.x -  cx_shift*boxlen)^2 +
                              (p.y -  cy_shift*boxlen )^2)
                              <= ( radius_out_shift*boxlen )  &&

                            abs(p.z - cz_shift*boxlen) <= ( height_shift*boxlen),
                                dataobject.data)
    elseif inverse == true
        sub_data = filter(p-> sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2)
                                    < ( radius_in_shift*boxlen ) ||

                              sqrt( (p.x -  cx_shift*boxlen)^2 +
                              (p.y -  cy_shift*boxlen )^2)
                              > ( radius_out_shift*boxlen )  ||

                            abs(p.z - cz_shift*boxlen) > ( height_shift*boxlen),
                                dataobject.data)
        ranges = dataobject.ranges
    end

    printtablememory(sub_data, verbose)

    partdata = PartDataType()
    partdata.data = sub_data
    partdata.info = dataobject.info
    partdata.lmin = dataobject.lmin
    partdata.lmax = dataobject.lmax
    partdata.boxlen = dataobject.boxlen
    partdata.ranges = ranges
    partdata.selected_partvars = dataobject.selected_partvars
    partdata.used_descriptors  = dataobject.used_descriptors
    partdata.scale = dataobject.scale
    return partdata

end

# -----------------------------------------------------------------------------
##### SPHERE/SHELL #####-------------------------------------------------------
function shellregionsphere(dataobject::PartDataType;
                            radius::Array{<:Real,1}=[0.,0.],
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    radius_in  = radius[1]
    radius_out = radius[2]
    if radius_in == 0. || radius_out == 0. || in(0., center)
        error("[Mera]: given inner and outer radius or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift = prep_spherical_shellranges(dataobject.info, center, radius_in, radius_out, range_unit, verbose)

    if inverse == false
        sub_data = filter(p-> sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2 +
                                    (p.z - cz_shift*boxlen)^2 )
                                    >= ( radius_in_shift*boxlen ) &&

                                    sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2 +
                                    (p.z - cz_shift*boxlen)^2 )
                                    <= ( radius_out_shift*boxlen ),
                                dataobject.data)
    elseif inverse == true
        sub_data = filter(p-> sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2 +
                                    (p.z - cz_shift*boxlen)^2 )
                                    < ( radius_in_shift*boxlen ) ||

                                    sqrt( (p.x -  cx_shift*boxlen)^2 +
                                    (p.y -  cy_shift*boxlen )^2 +
                                    (p.z - cz_shift*boxlen)^2 )
                                    > ( radius_out_shift*boxlen ),
                                dataobject.data)
        ranges = dataobject.ranges
    end

    printtablememory(sub_data, verbose)

    partdata = PartDataType()
    partdata.data = sub_data
    partdata.info = dataobject.info
    partdata.lmin = dataobject.lmin
    partdata.lmax = dataobject.lmax
    partdata.boxlen = dataobject.boxlen
    partdata.ranges = ranges
    partdata.selected_partvars = dataobject.selected_partvars
    partdata.used_descriptors  = dataobject.used_descriptors
    partdata.scale = dataobject.scale
    return partdata
end


text:
# Test variable initializers
@testset "Variable Initializers" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1])
    @infinite_variable(m, x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    data = DiscreteMeasureData(par, [0.5, 0.5], [0, 1])
    @constraint(m, c1, x + z - 2 <= 0)
    @constraint(m, c2, measure(x + y, data) - w == 0)
    @constraint(m, c3, x0 + y0 == 5)
    tm = optimizer_model(m)
    # test _initialize_hold_variables
    @testset "_initialize_hold_variables" begin
        @hold_variable(m)
        @test isa(IOTO._initialize_hold_variables(tm, m), Nothing)
        @test length(transcription_data(tm).hold_to_var) == 2
        @test transcription_variable(tm, z) isa VariableRef
        @test transcription_variable(tm, w) isa VariableRef
        @test name(transcription_variable(tm, z)) == name(z)
        @test name(transcription_variable(tm, w)) == name(w)
        @test has_lower_bound(transcription_variable(tm, z))
        @test has_upper_bound(transcription_variable(tm, z))
        @test is_binary(transcription_variable(tm, z))
        @test is_fixed(transcription_variable(tm, w))
        @test is_integer(transcription_variable(tm, w))
        @test start_value(transcription_variable(tm, w)) == 1.
    end
    # test _list_supports
    @testset "_list_supports" begin
        @test IOTO._list_supports((par, pars)) == [[0, 1], supports(pars)]
    end
    # test _make_supports
    @testset "_make_supports" begin
        expected = sort([(0, 0), (0, 1), (1, 0), (1, 1)])
        @test sort(IOTO._make_supports((par, par))) == expected
        @test IOTO._make_supports((par, pars)) isa Vector
    end
    # test _initialize_infinite_variables
    @testset "_initialize_infinite_variables" begin
        @infinite_variable(m, parameter_refs = (par,))
        @test isa(IOTO._initialize_infinite_variables(tm, m), Nothing)
        @test length(transcription_data(tm).infinite_to_vars) == 2
        @test transcription_variable(tm, x) isa Vector{VariableRef}
        @test transcription_variable(tm, y) isa Vector{VariableRef}
        @test name(transcription_variable(tm, x)[1]) == "x(support: 1)"
        @test name(transcription_variable(tm, y)[3]) == "y(support: 3)"
        @test has_lower_bound(transcription_variable(tm, x)[1])
        @test is_binary(transcription_variable(tm, y)[2])
        @test is_fixed(transcription_variable(tm, y)[4])
        @test is_integer(transcription_variable(tm, x)[2])
        @test start_value(transcription_variable(tm, y)[1]) == 0.
        @test supports(x) == [(0,), (1,)]
        @test length(supports(y)) == 4
    end
    # test _update_point_mapping
    @testset "_update_point_mapping" begin
        # test first point variable
        @test isa(IOTO._update_point_mapping(tm, x0, x, parameter_values(x0)),
                  Nothing)
        @test transcription_data(tm).point_to_var[x0] == transcription_variable(tm, x)[1]
        # test second point variable
        @test isa(IOTO._update_point_mapping(tm, y0, y, parameter_values(y0)),
                  Nothing)
        @test transcription_data(tm).point_to_var[y0] == transcription_variable(tm, y)[1]
        # test error
        @test_throws ErrorException IOTO._update_point_mapping(tm, x0, x, (0.5,))
    end
    # test _update_point_info
    @testset "_update_point_info" begin
        vref = transcription_variable(tm, x0)
        # test lower bound update
        fix(vref, 0, force = true)
        set_lower_bound(x0, 0)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test lower_bound(vref) == 0
        # test upper bound update
        fix(vref, 0, force = true)
        set_upper_bound(x0, 0)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test upper_bound(vref) == 0
        # test fix update
        fix(x0, 0, force = true)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test fix_value(vref) == 0
        # test binary update
        unset_integer(x0)
        set_binary(x0)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test is_binary(vref)
        # test integer update
        unset_binary(x0)
        set_integer(x0)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test is_integer(vref)
        # test start value
        set_start_value(x0, 0.5)
        @test isa(IOTO._update_point_info(tm, x0), Nothing)
        @test start_value(vref) == 0.5
        # Undo changes
        unfix(vref)
        unfix(x0)
        set_lower_bound(vref, 0)
        set_lower_bound(x0, 0)
        set_start_value(vref, NaN)
        set_start_value(x0, NaN)
        delete!(transcription_data(tm).point_to_var, x0)
        delete!(transcription_data(tm).point_to_var, y0)
    end
    # test _map_point_variables
    @testset "_map_point_variables" begin
        @test isa(IOTO._map_point_variables(tm, m), Nothing)
        @test length(transcription_data(tm).infinite_to_vars) == 2
        @test transcription_variable(tm, x0) isa VariableRef
        @test transcription_variable(tm, y0) isa VariableRef
        @test name(transcription_variable(tm, x0)) == "x(support: 1)"
        @test name(transcription_variable(tm, y0)) == "y(support: 1)"
        @test has_lower_bound(transcription_variable(tm, x0))
        @test is_integer(transcription_variable(tm, x0))
        @test has_lower_bound(transcription_variable(tm, y0))
        @test has_upper_bound(transcription_variable(tm, y0))
        @test is_integer(transcription_variable(tm, y0))
        @test start_value(transcription_variable(tm, y0)) == 0.
    end
end

# Test parameter search methods
@testset "Parameter Search Methods" begin
    # initialize model
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= x <= 1)
    @infinite_parameter(m, 0 <= y <= 1)
    @infinite_parameter(m, 0 <= z[1:3] <= 1, container = SparseAxisArray)
    supp1 = (1, convert(JuMPC.SparseAxisArray, [2, 3, 4]))
    tup1 = (y, z)
    supp2 = (1, 2, convert(JuMPC.SparseAxisArray, [3, 4, 5]))
    tup2 = (x, y, z)
    supp3 = (1, 2, JuMPC.SparseAxisArray(Dict((1,) => 3, (3,) => 5)))
    tup3 = (x, y, JuMPC.SparseAxisArray(Dict((1,) => z[1], (3,) => z[3])))
    # test _parameter_tuple_index
    @testset "_parameter_tuple_index" begin
        @test IOTO._parameter_tuple_index(x, tup1) == -1
        @test IOTO._parameter_tuple_index(x, tup2) == 1
        @test IOTO._parameter_tuple_index(y, tup2) == 2
        @test IOTO._parameter_tuple_index(z[2], tup2) == 3
        @test IOTO._parameter_tuple_index(z[2], tup3) === 3
        @test IOTO._parameter_tuple_index(z[3], tup3) == 3
    end
    # test _parameter_value
    @testset "_parameter_value" begin
        @test IOTO._parameter_value(x, supp1, tup1) === NaN
        @test IOTO._parameter_value(x, supp2, tup2) == 1
        @test IOTO._parameter_value(y, supp2, tup2) == 2
        @test IOTO._parameter_value(z[2], supp2, tup2) == 4
        @test IOTO._parameter_value(z[2], supp3, tup3) === NaN
        @test IOTO._parameter_value(z[3], supp3, tup3) == 5
    end
end

# Test "_map_to_variable"
@testset "_map_to_variable" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    index = m.next_var_index + 1
    m.reduced_info[index] = ReducedInfiniteInfo(y, Dict(1 => 1))
    yr = ReducedInfiniteVariableRef(m, index)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @variable(tm, d)
    @variable(tm, e)
    @variable(tm, f)
    @variable(tm, g)
    @variable(tm, h)
    supp1 = convert(JuMPC.SparseAxisArray, [0, 0])
    supp2 = convert(JuMPC.SparseAxisArray, [1, 1])
    # test for FiniteVariableRefs
    @testset "FiniteVariableRef" begin
        # test cannot find
        @test_throws ErrorException IOTO._map_to_variable(x0, (), (), tm)
        @test_throws ErrorException IOTO._map_to_variable(y0, (), (), tm)
        @test_throws ErrorException IOTO._map_to_variable(z, (), (), tm)
        @test_throws ErrorException IOTO._map_to_variable(w, (), (), tm)
        # prepare transcrition model with mapped variables
        tm.ext[:TransData].point_to_var[x0] = a
        tm.ext[:TransData].point_to_var[y0] = b
        tm.ext[:TransData].hold_to_var[z] = c
        tm.ext[:TransData].hold_to_var[w] = d
        # test the function
        @test IOTO._map_to_variable(x0, (), (), tm) == a
        @test IOTO._map_to_variable(y0, (), (), tm) == b
        @test IOTO._map_to_variable(z, (), (), tm) == c
        @test IOTO._map_to_variable(w, (), (), tm) == d
    end
    # test for InfiniteVariableRefs
    @testset "InfiniteVariableRef" begin
        # prepare transcrition model with mapped variables
        tm.ext[:TransData].infinite_to_vars[x] = [a, e]
        tm.ext[:TransData].infinite_to_vars[y] = [b, f, g, h]
        tm.ext[:TransData].infvar_to_supports[x] = [(0,), (1,)]
        tm.ext[:TransData].infvar_to_supports[y] = [(0, supp1), (0, supp2),
                                                    (1, supp1), (1, supp2)]
        # test cannot find
        @test_throws ErrorException IOTO._map_to_variable(x, (0.5, supp1),
                                                          (par, pars), tm)
        @test_throws ErrorException IOTO._map_to_variable(y, (0.5, supp1),
                                                          (par, pars), tm)
        # test the function
        @test IOTO._map_to_variable(x, (-0., supp1), (par, pars), tm) == a
        @test IOTO._map_to_variable(x, (1., supp1), (par, pars), tm) == e
        @test IOTO._map_to_variable(y, (-0., supp1), (par, pars), tm) == b
        @test IOTO._map_to_variable(y, (0, supp2), (par, pars), tm) == f
        @test IOTO._map_to_variable(y, (1., supp1), (par, pars), tm) == g
        @test IOTO._map_to_variable(y, (1, supp2), (par, pars), tm) == h
    end
    # test for ReducedInfiniteVariableRefs
    @testset "ReducedInfiniteVariableRef" begin
        # test cannot find
        supp3 = convert(JuMPC.SparseAxisArray, [0.5, 0.5])
        @test_throws ErrorException IOTO._map_to_variable(yr, (0.5, supp3),
                                                          (par, pars), tm)
        # test the function
        @test IOTO._map_to_variable(yr, (0., supp1), (par, pars), tm) == g
        @test IOTO._map_to_variable(yr, (supp1,), (pars,), tm) == g
        @test IOTO._map_to_variable(yr, (0., supp2), (par, pars), tm) == h
        @test IOTO._map_to_variable(yr, (supp2,), (pars,), tm) == h
    end
    # test for ParameterRefs
    @testset "ParameterRef" begin
        # test cannot find
        @test_throws ErrorException IOTO._map_to_variable(par, (supp1,),
                                                          (pars,), tm)
        # test the function
        @test IOTO._map_to_variable(par, (0., supp1), (par, pars), tm) == 0
        @test IOTO._map_to_variable(pars[2], (0., supp1), (par, pars), tm) == 0
    end
end

# Test _support(s)_in_bounds
@testset "Support Bound Checking" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    supp1 = convert(JuMPC.SparseAxisArray, [0, 0])
    supp2 = convert(JuMPC.SparseAxisArray, [1, 1])
    # test _support_in_bounds
    @testset "_support_in_bounds" begin
        # test in bounds
        @test IOTO._support_in_bounds((0, supp1), (par, pars), Dict())
        # test out of bounds
        bounds = Dict(par => IntervalSet(0.5, 1))
        @test !IOTO._support_in_bounds((0, supp1), (par, pars), bounds)
        bounds = Dict(par => IntervalSet(0, 0.5))
        @test !IOTO._support_in_bounds((.7, supp1), (par, pars), bounds)
        # test parameter bounds not relavent
        @test IOTO._support_in_bounds((supp1,), (pars,), bounds)
    end
    @testset "_supports_in_bounds" begin
        supps = [(0, supp1), (0, supp2), (1, supp1), (1, supp2)]
        bounds = Dict(par => IntervalSet(0.5, 1))
        @test IOTO._supports_in_bounds(supps, (par, pars), bounds) == [false,
                                                              false, true, true]
    end
end

# Test _truncate_supports
@testset "_make_transcription_function" begin
    # initialize the models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    m.next_var_index += 1
    m.reduced_info[m.next_var_index] = ReducedInfiniteInfo(y, Dict(1 => 1))
    yr = ReducedInfiniteVariableRef(m, m.next_var_index)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    meas3 = measure(y^2 + z, data1)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @variable(tm, d)
    @variable(tm, e)
    @variable(tm, f)
    @variable(tm, g)
    @variable(tm, h)
    supp1 = convert(JuMPC.SparseAxisArray, [0, 0])
    supp2 = convert(JuMPC.SparseAxisArray, [1, 1])
    # transcribe the variables
    tm.ext[:TransData].point_to_var[x0] = a
    tm.ext[:TransData].point_to_var[y0] = b
    tm.ext[:TransData].hold_to_var[z] = c
    tm.ext[:TransData].hold_to_var[w] = d
    tm.ext[:TransData].infinite_to_vars[x] = [a, e]
    tm.ext[:TransData].infinite_to_vars[y] = [b, f, g, h]
    tm.ext[:TransData].infvar_to_supports[x] = [(0,), (1,)]
    tm.ext[:TransData].infvar_to_supports[y] = [(0, supp1), (0, supp2),
                                                (1, supp1), (1, supp2)]
    # test FiniteVariableRefs
    @testset "FiniteVariableRef" begin
        # test point variables
        @test IOTO._make_transcription_function(x0, tm) == (a, ())
        @test IOTO._make_transcription_function(y0, tm) == (b, ())
        # test hold variables
        @test IOTO._make_transcription_function(z, tm) == (c, ())
        @test IOTO._make_transcription_function(w, tm) == (d, ())
    end
    # test InfiniteVariableRefs
    @testset "InfiniteVariableRef" begin
        # test without bounds
        @test IOTO._make_transcription_function(x, tm) == ([a, e], (par,),
                                                           [(0,), (1,)])
        @test IOTO._make_transcription_function(y, tm) == ([b, f, g, h],
                                                           (par, pars),
                                                           [(0, supp1), (0, supp2),
                                                            (1, supp1), (1, supp2)])
        # test with bounds
        bounds = Dict(par => IntervalSet(0.5, 1))
        @test IOTO._make_transcription_function(x, tm, bounds) == ([e], (par,),
                                                                   [(1,)])
        @test IOTO._make_transcription_function(y, tm, bounds) == ([g, h],
                                                                   (par, pars),
                                                                   [(1, supp1),
                                                                   (1, supp2)])
    end
    # test ParameterRefs
    @testset "ParameterRef" begin
        # test without bounds
        @test IOTO._make_transcription_function(par, tm) == ([0, 1], (par,), [0, 1])
        @test IOTO._make_transcription_function(pars[1], tm) == ([0, 1], (pars[1],),
                                                                 [0, 1])
        # test with bounds
        bounds = Dict(par => IntervalSet(0.5, 1))
        @test IOTO._make_transcription_function(par, tm, bounds) == ([1], (par,),
                                                                     [1])
        @test IOTO._make_transcription_function(pars[1], tm, bounds) == ([0, 1],
                                                             (pars[1],), [0, 1])
    end
    # test Finite GenericAffExprs
    @testset "Finite AffExpr" begin
        # prepare expressions
        aff1 = x0 + 2z - 3
        aff2 = 3z - w
        # test expressions
        @test IOTO._make_transcription_function(aff1, tm)[1] == a + 2c - 3
        @test IOTO._make_transcription_function(aff2, tm)[1] == 3c - d
    end
    # test Finite GenericQuadExprs
    @testset "Finite QuadExpr" begin
        # prepare expressions
        quad1 = x0^2 - y0 * w + 2z - 3
        quad2 = 2z * w + 3z - w
        # test expressions
        @test IOTO._make_transcription_function(quad1, tm)[1] == a^2 - b * d + 2c - 3
        @test IOTO._make_transcription_function(quad2, tm)[1] == 2c * d + 3c - d
    end
    # test General GenericAffExprs
    @testset "General AffExpr" begin
        # test only having 1 variable that is finite
        @test IOTO._make_transcription_function(z - 2, tm)[1] == c - 2
        @test IOTO._make_transcription_function(2w - w, tm)[1] == d + 0
        # test having only 1 variable that is infinite
        @test IOTO._make_transcription_function(x - 2, tm) == ([a - 2, e - 2],
                                                           (par,), [(0,), (1,)])
        @test IOTO._make_transcription_function(par - 2, tm) == ([-2., -1.],
                                                                 (par,), [0, 1])
        # test without bounds
        expr = meas1 + 2x - par + z
        expected = ([3a + e - 2d + c, a + 3e - 2d + c - 1], (par,), [(0,), (1,)])
        @test IOTO._make_transcription_function(expr, tm) == expected
        # test with bounds
        bounds = Dict(par => IntervalSet(0.5, 1))
        expr = meas2 + 2x - pars[2] + z
        expected = ([b + g + 2e + c, f + h + 2e + c - 1], (pars, par),
                    [(supp1, 1), (supp2, 1)])
        @test IOTO._make_transcription_function(expr, tm, bounds) == expected
    end
    # test General GenericQuadExprs
    @testset "General QuadExpr" begin
        # test only having 1 variable that is finite
        expr = zero(GenericQuadExpr{Float64, GeneralVariableRef})
        expr.aff = z - 2
        @test IOTO._make_transcription_function(expr, tm)[1] == c - 2
        expr.aff = 2w - w
        @test IOTO._make_transcription_function(expr, tm)[1] == d + 0
        # test having only 1 variable that is infinite
        expr.aff = x - 2
        @test IOTO._make_transcription_function(expr, tm) == ([a - 2, e - 2],
                                                           (par,), [(0,), (1,)])
        expr.aff = par - 2
        @test IOTO._make_transcription_function(expr, tm) == ([-2., -1.], (par,),
                                                              [0, 1])
        # test without bounds
        expr = x0^2 + 2 * z * y0 - par * x + w * par + 2 * par^2 - 2par + meas1 - 3
        expected = ([a^2 + 2 * c * b + a + e - 2d - 3,
                     a^2 + 2 * c * b + a + 2e - d - 3], (par,), [(0,), (1,)])
        expected[1][2].aff.terms[e] = 0. # hack to get 0e
        @test IOTO._make_transcription_function(expr, tm) == expected
        # test without bounds
        bounds = Dict(par => IntervalSet(0.5, 1))
        expr = x0^2 + 2 * z * y0 - pars[1] * x + w * par + 2 * par^2 - 2par + meas2 - 3
        expected = ([a^2 + 2 * c * b + b + g + d - 3,
                     a^2 + 2 * c * b + f + h - e + d - 3], (par, pars),
                     [(1, supp1), (1, supp2)])
        @test IOTO._make_transcription_function(expr, tm, bounds) == expected
        # test AffExpr dispatch to QuadExpr
        bounds = Dict(par => IntervalSet(0.5, 1))
        expr = 2meas3 - x + 2
        expected = ([2 * (b^2 + g^2 + 2c) - e + 2, 2 * (f^2 + h^2 + 2c) - e + 2],
                    (par, pars), [(1, supp1), (1, supp2)])
        @test IOTO._make_transcription_function(expr, tm, bounds) == expected
    end
    # test MeasureFiniteVariableRef Expressions
    @testset "MeasureFiniteVariableRef Exprs" begin
        # test finite AffExpr
        expr = meas1 + z - 2w
        expected = (a + e - 4d + c, ())
        @test IOTO._make_transcription_function(expr, tm) == expected
        # test infinite AffExpr
        expr = meas1 + x - 2w
        expected = ([2a + e - 4d, a + 2e - 4d], (par,), [(0,), (1,)])
        @test IOTO._make_transcription_function(expr, tm) == expected
        # test finite QuadExpr
        expr = meas1 + z^2 - 2w
        expected = (c^2 + a + e - 4d, ())
        @test IOTO._make_transcription_function(expr, tm) == expected
        # test infinite QuadExpr
        expr = meas1 + x^2 - 2w
        expected = ([a^2 + a + e - 4d, e^2 + a + e - 4d], (par,), [(0,), (1,)])
        @test IOTO._make_transcription_function(expr, tm) == expected
    end
    # test MeasrureRefs
    @testset "MeasureRef" begin
        @test IOTO._make_transcription_function(meas1, tm)[1] == a + e - 2d
        @test IOTO._make_transcription_function(meas2, tm)[1] == [b + g, f + h]
    end
    # test numeric AffExpr
    @testset "Numeric AffExpr" begin
        expr = zero(AffExpr) + 42
        @test IOTO._make_transcription_function(expr, tm) == (expr, ())
    end
    # test the fallback
    @testset "Fallback" begin
        expr = zero(QuadExpr)
        @test_throws ErrorException IOTO._make_transcription_function(expr, tm)
    end
end

# Test _set_objective
@testset "_set_objective" begin
    # initialize the models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @variable(tm, d)
    @variable(tm, e)
    @variable(tm, f)
    @variable(tm, g)
    @variable(tm, h)
    supp1 = convert(JuMPC.SparseAxisArray, [0, 0])
    supp2 = convert(JuMPC.SparseAxisArray, [1, 1])
    # transcribe the variables
    tm.ext[:TransData].hold_to_var[z] = c
    tm.ext[:TransData].hold_to_var[w] = d
    tm.ext[:TransData].infinite_to_vars[x] = [a, e]
    tm.ext[:TransData].infinite_to_vars[y] = [b, f, g, h]
    tm.ext[:TransData].infvar_to_supports[x] = [(0,), (1,)]
    tm.ext[:TransData].infvar_to_supports[y] = [(0, supp1), (0, supp2),
                                                (1, supp1), (1, supp2)]
    # test empty objective
    @test isa(IOTO._set_objective(tm, m), Nothing)
    @test objective_function(tm) == zero(AffExpr)
    # test finite objective
    @objective(m, Min, z - w)
    @test isa(IOTO._set_objective(tm, m), Nothing)
    @test objective_function(tm) == c - d
    @test objective_sense(tm) == MOI.MIN_SENSE
    # test measure objective
    @objective(m, Max, z - 2meas1)
    @test isa(IOTO._set_objective(tm, m), Nothing)
    @test objective_function(tm) == c - 2a - 2e + 4d
    @test objective_sense(tm) == MOI.MAX_SENSE
    # test error with infinite objective
    @objective(m, Max, z - 2meas2)
    @test_throws ErrorException IOTO._set_objective(tm, m)
end

# Test the helper functions for constraint transcription
@testset "Constraint Helpers" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, x(par) >= 0, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - z, data1)
    tm = optimizer_model(m)
    @variable(tm, xt[1:2])
    @variable(tm, zt)
    # Setup up the constraints
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @constraint(tm, c1t[i = 1:2], xt[i] + supports(par)[i] - zt == 0)
    @constraint(tm, c2t, zt >= -3)
    @constraint(tm, c3t, xt[1] + xt[2] - 3zt == 0)
    # Test _set_mapping (InfiniteConstraintRefs)
    @testset "_set_mapping (Infinite)" begin
        @test isa(IOTO._set_mapping(c1, c1t), Nothing)
        @test tm.ext[:TransData].infinite_to_constrs[c1] == [c1t[1], c1t[2]]
    end
    # Test _set_mapping (MeasureConstraintRefs that are infinite)
    @testset "_set_mapping (Infinite Measures)" begin
        @test isa(IOTO._set_mapping(c3, [c3t]), Nothing)
        @test tm.ext[:TransData].measure_to_constrs[c3] == [c3t]
    end
    # Test _set_mapping (MeasureConstraintRefs that are finite)
    @testset "_set_mapping (Finite Measures)" begin
        @test isa(IOTO._set_mapping(c3, c3t), Nothing)
        @test tm.ext[:TransData].measure_to_constrs[c3] == [c3t]
    end
    # Test _set_mapping (FiniteConstraintRefs)
    @testset "_set_mapping (Finite)" begin
        @test isa(IOTO._set_mapping(c2, c2t), Nothing)
        @test tm.ext[:TransData].finite_to_constr[c2] == c2t
    end
    # Test _set_supports (InfiniteConstraintRefs)
    @testset "_set_supports (Infinite)" begin
        @test isa(IOTO._set_supports(tm, c1, [(0,), (1,)]), Nothing)
        @test tm.ext[:TransData].infconstr_to_supports[c1] == [(0,), (1,)]
    end
    # Test _set_supports (MeasureConstraintRefs)
    @testset "_set_supports (Measure)" begin
        @test isa(IOTO._set_supports(tm, c3, [(0,), (1,)]), Nothing)
        @test tm.ext[:TransData].measconstr_to_supports[c3] == [(0,), (1,)]
    end
    # Test _set_parameter_refs (InfiniteConstraintRefs)
    @testset "_set_parameter_refs (Infinite)" begin
        @test isa(IOTO._set_parameter_refs(tm, c1, (par,)), Nothing)
        @test tm.ext[:TransData].infconstr_to_params[c1] == (par,)
    end
    # Test _set_parameter_refs (MeasureConstraintRefs)
    @testset "_set_parameter_refs (Measure)" begin
        @test isa(IOTO._set_parameter_refs(tm, c3, (par,)), Nothing)
        @test tm.ext[:TransData].measconstr_to_params[c3] == (par,)
    end
    # Test _root_name (constraints)
    @testset "_root_name (Constraints)" begin
        # test simple name return
        @test InfiniteOpt._root_name(c1) == "c1"
        # test anonymous name
        @test InfiniteOpt._root_name(@constraint(m, z == 1)) == "noname"
        # test array of constraints
        @test InfiniteOpt._root_name(@constraint(m, a[1:2], z == 1)[1]) == "a"
    end
end

# Test _set_constraints
@testset "_set_constraints" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @variable(tm, d)
    @variable(tm, e)
    @variable(tm, f)
    @variable(tm, g)
    @variable(tm, h)
    supp1 = convert(JuMPC.SparseAxisArray, [0, 0])
    supp2 = convert(JuMPC.SparseAxisArray, [1, 1])
    # transcribe the variables
    tm.ext[:TransData].point_to_var[x0] = a
    tm.ext[:TransData].point_to_var[y0] = b
    tm.ext[:TransData].hold_to_var[z] = c
    tm.ext[:TransData].hold_to_var[w] = d
    tm.ext[:TransData].infinite_to_vars[x] = [a, e]
    tm.ext[:TransData].infinite_to_vars[y] = [b, f, g, h]
    tm.ext[:TransData].infvar_to_supports[x] = [(0,), (1,)]
    tm.ext[:TransData].infvar_to_supports[y] = [(0, supp1), (0, supp2),
                                                (1, supp1), (1, supp2)]
    # Setup up the constraints
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @BDconstraint(m, c4(par in [0.5, 1]), meas2 - 2y0 + x <= 1)
    @constraint(m, c5, meas2 == 0)
    # test the main function
    @test isa(IOTO._set_constraints(tm, m), Nothing)
    @test length(list_of_constraint_types(tm)) == 3
    @test constraint_object(tm.ext[:TransData].infinite_to_constrs[c1][1]).func == a - c
    @test constraint_object(tm.ext[:TransData].infinite_to_constrs[c1][2]).func == e - c
    @test constraint_object(tm.ext[:TransData].finite_to_constr[c2]).func == c + a
    @test constraint_object(tm.ext[:TransData].infinite_to_constrs[c4][1]).func == -b + g + e
    @test name(tm.ext[:TransData].finite_to_constr[c2]) == "c2"
    @test name(tm.ext[:TransData].infinite_to_constrs[c1][1]) == "c1(Support: 1)"
    @test tm.ext[:TransData].infconstr_to_params[c1] == (par, )
    @test tm.ext[:TransData].infconstr_to_params[c4] == (pars, par)
    @test tm.ext[:TransData].infconstr_to_supports[c1] == [(0,), (1,)]
    @test tm.ext[:TransData].measconstr_to_params[c5] == (pars,)
end

# Test variable info mappers
@testset "Variable Info Mappers" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    tm = optimizer_model(m)
    IOTO._initialize_hold_variables(tm, m)
    IOTO._initialize_infinite_variables(tm, m)
    IOTO._map_point_variables(tm, m)
    xt = transcription_variable(tm, x)
    yt = transcription_variable(tm, y)
    zt = transcription_variable(tm, z)
    wt = transcription_variable(tm, w)
    # test _map_info_constraints (FiniteVariableRef)
    @testset "_map_info_constraints (Finite)" begin
        # test z
        @test isa(IOTO._map_info_constraints(z, tm), Nothing)
        @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(z)] == LowerBoundRef(zt)
        @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(z)] == UpperBoundRef(zt)
        @test tm.ext[:TransData].finite_to_constr[BinaryRef(z)] == BinaryRef(zt)
        # test w
        @test isa(IOTO._map_info_constraints(w, tm), Nothing)
        @test tm.ext[:TransData].finite_to_constr[FixRef(w)] == FixRef(wt)
        @test tm.ext[:TransData].finite_to_constr[IntegerRef(w)] == IntegerRef(wt)
    end
    # test _map_info_constraints (InfiniteVariableRef)
    @testset "_map_info_constraints (Infinite)" begin
        # test x
        @test isa(IOTO._map_info_constraints(x, tm), Nothing)
        @test tm.ext[:TransData].infinite_to_constrs[LowerBoundRef(x)] == LowerBoundRef.(xt)
        @test tm.ext[:TransData].infconstr_to_params[LowerBoundRef(x)] == parameter_refs(x)
        @test tm.ext[:TransData].infconstr_to_supports[LowerBoundRef(x)] == [(0,), (1,)]
        @test tm.ext[:TransData].infinite_to_constrs[UpperBoundRef(x)] == UpperBoundRef.(xt)
        @test tm.ext[:TransData].infconstr_to_params[UpperBoundRef(x)] == parameter_refs(x)
        @test tm.ext[:TransData].infconstr_to_supports[UpperBoundRef(x)] == [(0,), (1,)]
        @test tm.ext[:TransData].infinite_to_constrs[IntegerRef(x)] == IntegerRef.(xt)
        @test tm.ext[:TransData].infconstr_to_params[IntegerRef(x)] == parameter_refs(x)
        @test tm.ext[:TransData].infconstr_to_supports[IntegerRef(x)] == [(0,), (1,)]
        # test y
        @test isa(IOTO._map_info_constraints(y, tm), Nothing)
        @test tm.ext[:TransData].infinite_to_constrs[FixRef(y)] == FixRef.(yt)[2:end]
        @test tm.ext[:TransData].infconstr_to_params[FixRef(y)] == parameter_refs(y)
        @test tm.ext[:TransData].infconstr_to_supports[FixRef(y)] == supports(y)[2:end]
        @test tm.ext[:TransData].infinite_to_constrs[BinaryRef(y)] == [BinaryRef(yt[i]) for i = 2:4]
        @test tm.ext[:TransData].infconstr_to_params[BinaryRef(y)] == parameter_refs(y)
        @test tm.ext[:TransData].infconstr_to_supports[BinaryRef(y)] == supports(y)[2:end]
    end
end

# Test _map_variable_info_constraints
@testset "_map_variable_info_constraints" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    tm = optimizer_model(m)
    IOTO._initialize_hold_variables(tm, m)
    IOTO._initialize_infinite_variables(tm, m)
    IOTO._map_point_variables(tm, m)
    xt = transcription_variable(tm, x)
    yt = transcription_variable(tm, y)
    zt = transcription_variable(tm, z)
    wt = transcription_variable(tm, w)
    # test running the function
    @test isa(IOTO._map_variable_info_constraints(tm, m), Nothing)
    # test z
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(z)] == LowerBoundRef(zt)
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(z)] == UpperBoundRef(zt)
    @test tm.ext[:TransData].finite_to_constr[BinaryRef(z)] == BinaryRef(zt)
    # test w
    @test tm.ext[:TransData].finite_to_constr[FixRef(w)] == FixRef(wt)
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(w)] == IntegerRef(wt)
    # test x0
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(x0)] == LowerBoundRef(xt[1])
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(x0)] == UpperBoundRef(xt[1])
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(x0)] == IntegerRef(xt[1])
    # test y0
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(y0)] == LowerBoundRef(yt[1])
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(y0)] == UpperBoundRef(yt[1])
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(y0)] == IntegerRef(yt[1])
    # test x
    @test tm.ext[:TransData].infinite_to_constrs[LowerBoundRef(x)] == LowerBoundRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[LowerBoundRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[LowerBoundRef(x)] == [(0,), (1,)]
    @test tm.ext[:TransData].infinite_to_constrs[UpperBoundRef(x)] == UpperBoundRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[UpperBoundRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[UpperBoundRef(x)] == [(0,), (1,)]
    @test tm.ext[:TransData].infinite_to_constrs[IntegerRef(x)] == IntegerRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[IntegerRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[IntegerRef(x)] == [(0,), (1,)]
    # test y
    @test tm.ext[:TransData].infinite_to_constrs[FixRef(y)] == FixRef.(yt)[2:end]
    @test tm.ext[:TransData].infconstr_to_params[FixRef(y)] == parameter_refs(y)
    @test tm.ext[:TransData].infconstr_to_supports[FixRef(y)] == supports(y)[2:end]
    @test tm.ext[:TransData].infinite_to_constrs[BinaryRef(y)] == [BinaryRef(yt[i]) for i = 2:4]
    @test tm.ext[:TransData].infconstr_to_params[BinaryRef(y)] == parameter_refs(y)
    @test tm.ext[:TransData].infconstr_to_supports[BinaryRef(y)] == supports(y)[2:end]
end

# Test TranscriptionModel constructor
@testset "TranscriptionModel (Full Build)" begin
    # initialize model
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1],
                        container = SparseAxisArray)
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    @finite_parameter(m, fin, 0)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    @set_parameter_bounds(z, (pars == 0, par == 0))
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @BDconstraint(m, c4(par in [0.5, 1]), meas2 - 2y0 + x + fin <= 1)
    @constraint(m, c5, meas2 == 0)
    @objective(m, Min, x0 + meas1)
    # test basic usage
    @test isa(set_optimizer_model(m, TranscriptionModel(m)), Nothing)
    tm = optimizer_model(m)
    # test hold variables
    @test length(transcription_data(tm).hold_to_var) == 2
    @test transcription_variable(tm, z) isa VariableRef
    @test transcription_variable(tm, w) isa VariableRef
    @test name(transcription_variable(tm, z)) == name(z)
    @test name(transcription_variable(tm, w)) == name(w)
    @test has_lower_bound(transcription_variable(tm, z))
    @test has_upper_bound(transcription_variable(tm, z))
    @test is_binary(transcription_variable(tm, z))
    @test is_fixed(transcription_variable(tm, w))
    @test is_integer(transcription_variable(tm, w))
    @test start_value(transcription_variable(tm, w)) == 1.
    # test infinite variables
    @test length(transcription_data(tm).infinite_to_vars) == 2
    @test transcription_variable(tm, x) isa Vector{VariableRef}
    @test transcription_variable(tm, y) isa Vector{VariableRef}
    @test name(transcription_variable(tm, x)[1]) == "x(support: 1)"
    @test name(transcription_variable(tm, y)[3]) == "y(support: 3)"
    @test has_lower_bound(transcription_variable(tm, x)[1])
    @test is_binary(transcription_variable(tm, y)[2])
    @test is_fixed(transcription_variable(tm, y)[4])
    @test is_integer(transcription_variable(tm, x)[2])
    @test start_value(transcription_variable(tm, y)[1]) == 0.
    @test supports(x) == [(0,), (1,)]
    @test length(supports(y)) == 4
    # test point variables
    @test length(transcription_data(tm).infinite_to_vars) == 2
    @test transcription_variable(tm, x0) isa VariableRef
    @test transcription_variable(tm, y0) isa VariableRef
    @test name(transcription_variable(tm, x0)) == "x(support: 1)"
    @test name(transcription_variable(tm, y0)) == "y(support: 1)"
    @test has_lower_bound(transcription_variable(tm, x0))
    @test is_integer(transcription_variable(tm, x0))
    @test has_lower_bound(transcription_variable(tm, y0))
    @test has_upper_bound(transcription_variable(tm, y0))
    @test is_integer(transcription_variable(tm, y0))
    @test start_value(transcription_variable(tm, y0)) == 0.
    # test objective
    xt = transcription_variable(tm, x)
    wt = transcription_variable(tm, w)
    @test objective_function(tm) == 2xt[1] + xt[2] - 2wt
    @test objective_sense(tm) == MOI.MIN_SENSE
    # test constraints
    zt = transcription_variable(tm, z)
    yt = transcription_variable(tm, y)
    @test constraint_object(tm.ext[:TransData].infinite_to_constrs[c1][1]).func == xt[1] - zt
    @test constraint_object(tm.ext[:TransData].finite_to_constr[c2]).func == zt + xt[1]
    @test constraint_object(tm.ext[:TransData].infinite_to_constrs[c4][1]).func == -yt[1] + yt[2] + xt[2]
    @test name(tm.ext[:TransData].finite_to_constr[c2]) == "c2"
    @test name(tm.ext[:TransData].infinite_to_constrs[c1][1]) == "c1(Support: 1)"
    @test tm.ext[:TransData].infconstr_to_params[c1] == (par, )
    @test tm.ext[:TransData].infconstr_to_params[c4] == (pars, par, fin)
    @test tm.ext[:TransData].infconstr_to_supports[c1] == [(0,)]
    @test tm.ext[:TransData].measconstr_to_params[c5] == (pars,)
    # test info constraints
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(z)] == LowerBoundRef(zt)
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(z)] == UpperBoundRef(zt)
    @test tm.ext[:TransData].finite_to_constr[BinaryRef(z)] == BinaryRef(zt)
    @test tm.ext[:TransData].finite_to_constr[FixRef(w)] == FixRef(wt)
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(w)] == IntegerRef(wt)
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(x0)] == LowerBoundRef(xt[1])
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(x0)] == UpperBoundRef(xt[1])
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(x0)] == IntegerRef(xt[1])
    @test tm.ext[:TransData].finite_to_constr[LowerBoundRef(y0)] == LowerBoundRef(yt[1])
    @test tm.ext[:TransData].finite_to_constr[UpperBoundRef(y0)] == UpperBoundRef(yt[1])
    @test tm.ext[:TransData].finite_to_constr[IntegerRef(y0)] == IntegerRef(yt[1])
    @test tm.ext[:TransData].infinite_to_constrs[LowerBoundRef(x)] == LowerBoundRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[LowerBoundRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[LowerBoundRef(x)] == [(0,), (1,)]
    @test tm.ext[:TransData].infinite_to_constrs[UpperBoundRef(x)] == UpperBoundRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[UpperBoundRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[UpperBoundRef(x)] == [(0,), (1,)]
    @test tm.ext[:TransData].infinite_to_constrs[IntegerRef(x)] == IntegerRef.(xt)
    @test tm.ext[:TransData].infconstr_to_params[IntegerRef(x)] == parameter_refs(x)
    @test tm.ext[:TransData].infconstr_to_supports[IntegerRef(x)] == [(0,), (1,)]
    @test tm.ext[:TransData].infinite_to_constrs[FixRef(y)] == FixRef.(yt)[2:end]
    @test tm.ext[:TransData].infconstr_to_params[FixRef(y)] == parameter_refs(y)
    @test tm.ext[:TransData].infconstr_to_supports[FixRef(y)] == supports(y)[2:end]
    @test tm.ext[:TransData].infinite_to_constrs[BinaryRef(y)] == [BinaryRef(yt[i]) for i = 2:4]
    @test tm.ext[:TransData].infconstr_to_params[BinaryRef(y)] == parameter_refs(y)
    @test tm.ext[:TransData].infconstr_to_supports[BinaryRef(y)] == supports(y)[2:end]
    # test with solver
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    @test isa(set_optimizer_model(m, TranscriptionModel(m, mockoptimizer)), Nothing)
end


text:
####### SIMULATION OF MARGINAL ###########################################

### +++ INITIALIZATION +++

### import packages
using Random, LinearAlgebra, RDatasets, MixedModels, Distributions, Optim, CSV, Dates

### define functions
vcov2 = function(x, xbar, σ, τ, n, γ)
    m = length(xbar)
    a = γ / τ^2 * m # 1-1 entry
    d = (sum(x.^2) - n * γ * sum(xbar.^2) ) / σ^2 # 2-2 entry
    b = γ / τ^2 / n * sum(x) # 2-1 entry
    [d -b; -b a] / (a * d - b^2) # return inverse (X^t V^-1 X)^-1
end

beta = function(vcov1, x, xbar, y, ybar, σ, τ, n, γ)
    c1 = γ / τ^2 * sum(ybar) # = 1^t V^-1 y
    c2 = (x' * y - n * γ * sum(ybar .* xbar))/σ^2 # = x^t V^-1 y
    vcov1 * [c1 c2]'  # (X^tV^-1X)^-1 X^tV^-1y
end

### fix parameters
p = 2
nsim = 10000
Δ = [0:0.1:2;]
Random.seed!(12345)

### +++ PER SPECIFICATION SETTING +++

### to compile
τ1 = 8^0.5
σ1 = 4^0.5
m = 10
n = 1

### variable parameters
N = m * n

β0 = [-4.9, 0.03] # rand(Uniform(0, 1), p)
x = rand(Uniform(0, 10), N)
xbar = vec(sum(reshape(x, n, :), dims=1) / n)

### +++ KNOWN DELTA +++

###     MARGINAL COVARIANCE MATRIX

### helpful parameters
γ0 = τ1^2 / (τ1^2 + σ1^2 / n)
vcov0 = vcov2(x, xbar, σ1, τ1, n, γ0) # = (X^tV^-1X)^-1
d0 = (1 - γ0) .* hcat(fill(1, m), xbar) # d' = l' - b'X
dd0 = d0' * d0

### assemble covariance matrix parts
g1 = γ0 * σ1^2 / n

### assemble inverse of covariance matrix with Woodbury
Σ0diag = g1
Σ0i = (I - d0 * inv(Σ0diag * inv(vcov0) + dd0) * d0' ) / Σ0diag

###     QUANTILES
qm = cquantile(Chisq(m), 0.05) # quantile of alpha = 0.05

### simulation
Coverage = zeros(2, length(Δ))
for i in 1:nsim

    ### +++ EMPIRICAL +++

    ### marginal evaluation
    v = τ1 * randn(m)
    u = β0[1] .+ β0[2] * x + repeat(v, inner = n)
    μ0 = β0[1] .+ β0[2] .* xbar + v

    ### variable data
    y = u + σ1 * randn(N) # conditional
    ybar = vec(sum(reshape(y, n, :), dims=1) / n)
    data = DataFrame(x = x, v = repeat(categorical(v), inner = n), y = y)

    ### fit model
    model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)

    ### make sure the variance parameters are positive
    τ = maximum([VarCorr(model).σρ.v.σ[1], 0.001])
    σ = maximum([VarCorr(model).s, 0.001])

    ###     MARGINAL

    ### helpful parameters
    γ = τ^2 / (τ^2 + σ^2 / n)
    vcov1 = vcov2(x, xbar, σ, τ, n, γ) # = (X^tV^-1X)^-1
    d = (1-γ) .* hcat(fill(1, m), xbar) # d' = l' - b'X
    dd = d' * d
    a = n * τ^2 + σ^2

    ### assemble diagonal matrix parts
    g1 = γ * σ^2 / n
    g3 = 2 * σ^4 / τ^2 * γ / N/ (n - 1)

    ### assemble inverse of covariance matrix with Woodbury
    Σdiag = g1 + 2 * g3
    #Σ = Σdiag * I + d * vcov1 * d'
    Σi = (I - d * inv(Σdiag * inv(vcov1) + dd) * d' ) / Σdiag

    ###     COMPUTE MEAN
    β1 = beta(vcov0, x, xbar, y, ybar, σ1, τ1, n, γ0)
    μ1 = (1 - γ0) * (β1[1] .+ β1[2] .* xbar) + γ0 * ybar

    β2 = beta(vcov1, x, xbar, y, ybar, σ, τ, n, γ) # = fixef(model)
    μ2 = (1 - γ) * (β2[1] .+ β2[2] .* xbar) + γ * ybar

    ###     COMPUTE STATISTICS
    t1 = (μ1 - μ0) .+ Δ'
    t2 = (μ2 - μ0) .+ Δ'
    Coverage[1,:] += (sum(t1' .* (t1' * Σ0i), dims = 2) .< qm) / nsim
    Coverage[2,:] += (sum(t2' .* (t2' * Σi), dims = 2) .< qm) / nsim
end


text:
@testset "Metadata" begin
    sample_text1 = "This is a string"
    sample_text2 = "This is also a string"

    sd1 = StringDocument(sample_text1)
    sd2 = StringDocument(sample_text2)
    crps = Corpus([sd1, sd2])

    # Single document metadata getters
    @test isequal(StringAnalysis.name(sd1), "")
    @test isequal(StringAnalysis.language(sd1), Languages.English())
    @test isequal(StringAnalysis.author(sd1), "")
    @test isequal(StringAnalysis.timestamp(sd1), "")
    @test isequal(StringAnalysis.id(sd1), "")
    @test isequal(StringAnalysis.publisher(sd1), "")
    @test isequal(StringAnalysis.published_year(sd1), "")
    @test isequal(StringAnalysis.edition_year(sd1), "")
    @test isequal(StringAnalysis.documenttype(sd1), "")
    @test isequal(StringAnalysis.note(sd1), "")

    # Single document metadata setters
    StringAnalysis.name!(sd1, "Document")
    StringAnalysis.language!(sd1, Languages.German())
    StringAnalysis.author!(sd1, "Author")
    StringAnalysis.timestamp!(sd1, "Time")
    StringAnalysis.id!(sd1, "ID")
    StringAnalysis.publisher!(sd1, "Publisher")
    StringAnalysis.published_year!(sd1, "Publishing Year")
    StringAnalysis.edition_year!(sd1, "Edition Year")
    StringAnalysis.documenttype!(sd1, "Type")
    StringAnalysis.note!(sd1, "A Note")

    @test isequal(StringAnalysis.name(sd1), "Document")
    @test isequal(StringAnalysis.language(sd1), Languages.German())
    @test isequal(StringAnalysis.author(sd1), "Author")
    @test isequal(StringAnalysis.timestamp(sd1), "Time")
    @test isequal(StringAnalysis.id(sd1), "ID")
    @test isequal(StringAnalysis.publisher(sd1), "Publisher")
    @test isequal(StringAnalysis.published_year(sd1), "Publishing Year")
    @test isequal(StringAnalysis.edition_year(sd1), "Edition Year")
    @test isequal(StringAnalysis.documenttype(sd1), "Type")
    @test isequal(StringAnalysis.note(sd1), "A Note")

    # Metadata getters for an entire corpus
    @test isequal(StringAnalysis.names(crps), ["Document", ""])
    @test isequal(StringAnalysis.languages(crps), [Languages.German(), Languages.English()])
    @test isequal(StringAnalysis.authors(crps), ["Author", ""])
    @test isequal(StringAnalysis.timestamps(crps), ["Time", ""])
    @test isequal(StringAnalysis.ids(crps), ["ID", ""])
    @test isequal(StringAnalysis.publishers(crps), ["Publisher", ""])
    @test isequal(StringAnalysis.published_years(crps), ["Publishing Year", ""])
    @test isequal(StringAnalysis.edition_years(crps), ["Edition Year", ""])
    @test isequal(StringAnalysis.documenttypes(crps), ["Type", ""])
    @test isequal(StringAnalysis.notes(crps), ["A Note", ""])

    # Metadata setters for an entire corpus
    StringAnalysis.names!(crps, "Document")
    StringAnalysis.languages!(crps, Languages.Spanish())
    StringAnalysis.authors!(crps, "Author")
    StringAnalysis.timestamps!(crps, "Time")
    StringAnalysis.ids!(crps, "ID")
    StringAnalysis.publishers!(crps, "Publisher")
    StringAnalysis.published_years!(crps, "Publishing Year")
    StringAnalysis.edition_years!(crps, "Edition Year")
    StringAnalysis.documenttypes!(crps, "Type")
    StringAnalysis.notes!(crps, "A Note")
    @test isequal(StringAnalysis.names(crps), ["Document", "Document"])
    @test isequal(StringAnalysis.languages(crps), [Languages.Spanish(), Languages.Spanish()])
    @test isequal(StringAnalysis.authors(crps), ["Author", "Author"])
    @test isequal(StringAnalysis.timestamps(crps), ["Time", "Time"])
    @test isequal(StringAnalysis.ids(crps), ["ID", "ID"])
    @test isequal(StringAnalysis.publishers(crps), ["Publisher", "Publisher"])
    @test isequal(StringAnalysis.published_years(crps), ["Publishing Year", "Publishing Year"])
    @test isequal(StringAnalysis.edition_years(crps), ["Edition Year", "Edition Year"])
    @test isequal(StringAnalysis.documenttypes(crps), ["Type", "Type"])
    @test isequal(StringAnalysis.notes(crps), ["A Note", "A Note"])

    StringAnalysis.names!(crps, ["Unnamed Document", "Unnamed Document"])
    StringAnalysis.languages!(crps, [Languages.English(), Languages.English()])
    StringAnalysis.authors!(crps, ["Unknown Author", "Unknown Author"])
    StringAnalysis.timestamps!(crps, ["Unknown Time", "Unknown Time"])
    StringAnalysis.ids!(crps, ["Unknown ID", "Unknown ID"])
    StringAnalysis.publishers!(crps, ["Unknown Publisher", "Unknown Publisher"])
    StringAnalysis.published_years!(crps, ["Unknown Publishing Year", "Unknown Publishing Year"])
    StringAnalysis.edition_years!(crps, ["Unknown Edition Year", "Unknown Edition Year"])
    StringAnalysis.documenttypes!(crps, ["Unknown Type", "Unknown Type"])
    StringAnalysis.notes!(crps, ["A Note", "Another Note"])
    @test isequal(StringAnalysis.names(crps), ["Unnamed Document", "Unnamed Document"])
    @test isequal(StringAnalysis.languages(crps), [Languages.English(), Languages.English()])
    @test isequal(StringAnalysis.authors(crps), ["Unknown Author", "Unknown Author"])
    @test isequal(StringAnalysis.timestamps(crps), ["Unknown Time", "Unknown Time"])
    @test isequal(StringAnalysis.ids(crps), ["Unknown ID", "Unknown ID"])
    @test isequal(StringAnalysis.publishers(crps), ["Unknown Publisher", "Unknown Publisher"])
    @test isequal(StringAnalysis.published_years(crps), ["Unknown Publishing Year", "Unknown Publishing Year"])
    @test isequal(StringAnalysis.edition_years(crps), ["Unknown Edition Year", "Unknown Edition Year"])
    @test isequal(StringAnalysis.documenttypes(crps), ["Unknown Type", "Unknown Type"])
    @test isequal(StringAnalysis.notes(crps), ["A Note", "Another Note"])
end


text:
using Random
using Distributions
using StaticArrays

mutable struct CartPole <: AbstractEnvironment{SizedVector{4,Float64},Int}
    max_steps::Int
    steps::Int
    state::Vector{Float64}
    done::Bool
    rng::MersenneTwister

    function CartPole(max_steps::Int; seed=Union{Nothing,Int} = nothing)
        max_steps < 1 && throw(DomainError(max_steps, "max_steps may not be lesser then 1"))
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
        rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

        new(max_steps, 0, Vector{Float64}(undef, 4), false, rng)
    end
end

CartPole(; steps, seed) = CartPole(steps; seed)
CartPoleV1(;seed=nothing) = CartPole(500; seed)

    function reset!(env::CartPole)
    rand!(env.rng, Uniform(-0.05, 0.05), env.state)
    env.done = false
    env.steps = 0

    return 0.0, SizedVector{4}(env.state), env.done
end

observation_type(::CartPole) = SizedVector{4,Float64}
observation_extrema(::CartPole) = (SVector(-4.8, -4.0, -0.418, -3.0), SVector(4.8, 4.0, 0.418, 3.0))
observation_length(::CartPole) = 4
action_type(::CartPole) = Int
action_set(::CartPole) = SA[-1,1]

# TODO: Check if env has been reseted
function step!(env::CartPole, action::Int)
    env.done && error("Cannot step after episode's end. Reset environment")
    (action != -1 && action != 1) && throw(DomainError(action, "Invalid action"))

    gravity = 9.8
    masscart = 1.0
    masspole = 0.1
    total_mass = masspole + masscart
    length = 0.5  # actually half the pole's length
    polemass_length = masspole * length
    force_mag = 10.0
    tau = 0.02  # seconds between state updates

  # Angle at which to fail the episode
    theta_threshold_radians = 12 * 2 * pi / 360
    x_threshold = 2.4

    x, x_dot, theta, theta_dot = env.state
    force = (action == 1) ? force_mag : -force_mag
    costheta = cos(theta)
    sintheta = sin(theta)

    temp = (force + polemass_length * theta_dot^2 * sintheta) / total_mass
    thetaacc = (gravity * sintheta - costheta * temp) / (length * (4.0 / 3.0 - masspole * costheta^2 / total_mass))
    xacc = temp - polemass_length * thetaacc * costheta / total_mass

    env.state[1] += tau * x_dot
    env.state[2] += tau * xacc
    env.state[3] += tau * theta_dot
    env.state[4] += tau * thetaacc
    env.steps += 1
    env.done = x < -x_threshold ||
               x > x_threshold ||
               theta < -theta_threshold_radians ||
               theta > theta_threshold_radians ||
               env.steps ≥ env.max_steps

    # If no copy is made, all transitions in a buffer are going to end up
    # pointing to the same state.
    # The copying might be unecessary if I were to use StructArrays in the
    # buffer (???)
    # TODO: Think of a more elegant and/or efficient solution for the problem
    #       described above.
    return 1.0, SizedVector{4}(env.state[:]), env.done
end


text:
"""
    Student's t

* Parametrization
parametrized in \\nu

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
TDist

function score!(score_til::Matrix{T}, y::T, ::Type{TDist}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = 0.5 * (((y^2)*(param[t, 1] + 1)/(param[t, 1] * y^2 + param[t, 1]^2)) -
                      1/param[t, 1] -
                      log((y^2)/(param[t, 1]) + 1) -
                      digamma(param[t, 1] / 2) + digamma((param[t, 1] + 1) / 2))
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T
    aux.fisher[1, 1] = 0.5 * (0.5 * trigamma(0.5 * param[t, 1]) - 0.5 * trigamma( 0.5 * (param[t, 1] + 1.0) ) -
                             (param[t, 1] + 5.0) / (param[t, 1] * (param[t, 1] + 3.0) * (param[t, 1] + 1.0)))
    return
end

function log_likelihood(::Type{TDist}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = zero(T)
    for t in 1:n
        loglik -= 0.5 * log(param[t, 1]) + logbeta(0.5, param[t, 1]/2) + ((param[t, 1] + 1)/2) *
                  log(1 + (y[t]^2)/param[t, 1])
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{TDist}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
    return
end

# utils
function update_dist(::Type{TDist}, param::Matrix{T}, t::Int) where T
    return TDist(param[t, 1])
end

function params_sdm(d::TDist)
    return Distributions.params(d)
end

function num_params(::Type{TDist})
    return 1
end


text:
module AbstractMergers

export AbstractMerger

"""
    AbstractMerger

Abstract base type for mergers.

# Interface

Any subtype should support to combine a forward stream with other streams that
can be accessed through a state dictionary via their `Splitter` name.
"""
abstract type AbstractMerger end

end # module AbstractMergers


text:
module Ditherings

using Images
export Quantise
export ZeroOne
export ZeroOne_PerChannel
export FloydSteinbergDither4Sample
export FloydSteinbergDither12Sample

function Quantise(pixel)
    shift = 4
    scale = 255.0
    r = Int64.(round(scale * (red(pixel)  )))>>shift
    g = Int64.(round(scale * (green(pixel))))>>shift
    b = Int64.(round(scale * (blue(pixel) )))>>shift
    r = Float64.((r<<shift)/scale)
    g = Float64.((g<<shift)/scale)
    b = Float64.((b<<shift)/scale)
    RGB(r,g,b)
end

function ZeroOne(pixel)
    r = red(pixel)
    g = green(pixel)
    b = blue(pixel)

    if(r + b + g >= 1.5)
        return RGB(1.0, 1.0, 1.0)
    else
        return RGB(0.0, 0.0, 0.0)
    end
end

function ZeroOne_PerChannel(pixel)
    r = red(pixel)
    g = green(pixel)
    b = blue(pixel)

    r= r>0.5 ? 1.0 : 0.0
    g= g>0.5 ? 1.0 : 0.0
    b= b>0.5 ? 1.0 : 0.0

    return RGB(r, g, b)

end

function FloydSteinbergDither4Sample(_img, PaletteFunction, weights=[7/16,3/16,5/16,1/16])
    width, height = size(_img);
    inputtype = typeof(_img)
    newimg = RGB{Float64}.(_img)

    for y in 2:height-1
        for x in 2:width-1
            oldpix = newimg[x,y]
            newpix = PaletteFunction(oldpix)
            newimg[x,y] = newpix
            quant_error = oldpix - newpix
            newimg[x+1,y  ] = newimg[x+1,y  ] + (quant_error * weights[1])
            newimg[x-1,y+1] = newimg[x-1,y+1] + (quant_error * weights[2])
            newimg[x  ,y+1] = newimg[x  ,y+1] + (quant_error * weights[3])
            newimg[x+1,y+1] = newimg[x+1,y+1] + (quant_error * weights[4])
        end
    end

    newimg
end

function FloydSteinbergDither12Sample(_img, PaletteFunction, weights = [0.243228, 0.0810761, 0.0417918, 0.0573652, 0.243228, 0.0573652, 0.0417918, 0.0347469, 0.0417918, 0.0810761, 0.0417918, 0.0347469])
    width, height = size(_img);
    inputtype = typeof(_img)
    newimg = RGB{Float64}.(_img)

    for y in 3:height-2
        for x in 3:width-2
            oldpix = newimg[x,y]
            newpix = PaletteFunction(oldpix)
            newimg[x,y] = newpix
            quant_error = oldpix - newpix
            newimg[x+1,y  ] = newimg[x+1,y  ] + (quant_error * weights[1] )
            newimg[x+2,y  ] = newimg[x+2,y  ] + (quant_error * weights[2] )

            newimg[x-2,y+1] = newimg[x-2,y+1] + (quant_error * weights[3] )
            newimg[x-1,y+1] = newimg[x-1,y+1] + (quant_error * weights[4] )
            newimg[x  ,y+1] = newimg[x  ,y+1] + (quant_error * weights[5] )
            newimg[x+1,y+1] = newimg[x+1,y+1] + (quant_error * weights[6] )
            newimg[x+2,y+1] = newimg[x+2,y+1] + (quant_error * weights[7] )

            newimg[x-2,y+2] = newimg[x-2,y+2] + (quant_error * weights[8] )
            newimg[x-1,y+2] = newimg[x-1,y+2] + (quant_error * weights[9] )
            newimg[x  ,y+2] = newimg[x  ,y+2] + (quant_error * weights[10] )
            newimg[x+1,y+2] = newimg[x+1,y+2] + (quant_error * weights[11] )
            newimg[x+2,y+2] = newimg[x+2,y+2] + (quant_error * weights[12] )
        end
    end

    newimg
end

function HorizontalDither3Sample(_img, PaletteFunction, weights=[0.5,0.4,0.1])
    width, height = size(_img);
    inputtype = typeof(_img)
    newimg = RGB{Float64}.(_img)

    for y in 2:height-3
        for x in 2:width-3
            oldpix = newimg[x,y]
            newpix = PaletteFunction(oldpix)
            newimg[x,y] = newpix
            quant_error = oldpix - newpix
            newimg[x+1,y  ] = newimg[x+1,y  ] + (quant_error * weights[1])
            newimg[x+2,y  ] = newimg[x+2,y  ] + (quant_error * weights[2])
            newimg[x+3,y  ] = newimg[x+3,y  ] + (quant_error * weights[3])
        end
    end

    newimg
end

function VerticalDither3Sample(_img, PaletteFunction, weights=[0.5,0.4,0.1])
    width, height = size(_img);
    inputtype = typeof(_img)
    newimg = RGB{Float64}.(_img)

    for y in 2:height-3
        for x in 2:width-3
            oldpix = newimg[x,y]
            newpix = PaletteFunction(oldpix)
            newimg[x,y] = newpix
            quant_error = oldpix - newpix
            newimg[x,y+1  ] = newimg[x,y+1  ] + (quant_error * weights[1])
            newimg[x,y+2  ] = newimg[x,y+2  ] + (quant_error * weights[2])
            newimg[x,y+3  ] = newimg[x,y+3  ] + (quant_error * weights[3])
        end
    end

    newimg
end

end


text:
# -----------------------------------------------------------
#
# Implementing a nodal basis for DG sparse grids,
# Multidimensional case
#
# -----------------------------------------------------------

#=

	We continue the work started in 1D_Nodal_Basis.jl by applying the
	tensor product construction together with the sparse cutoff
	to constract a way to translate between the modal basis used in
	linear wave evolution and the nodal/point basis necessary for
	efficiently handling nonlinearities in the time evolution equations

=#

# Check whether (level1, cell1) overlaps with (level2, cell2)
function relevant_cell_1D(k::Int, level1::Int, cell1::Int, level2::Int, cell2::Int)
    left1  = (cell1-1)/(1<<max(0, level1-1))
    right1 = cell1/(1<<max(0, level1-1))
    left2  = (cell2-1)/(1<<max(0, level2-1))
    right2 = cell2/(1<<max(0, level2-1))
    return (left1 <= left2 && right1 >= right2) || (left1 >= left2 && right1 <= right2)
end

# Multidimensional version of the above
function relevant_cell{D}(k::Int, level1::CartesianIndex{D}, cell1::CartesianIndex{D},
								  level2::CartesianIndex{D}, cell2::CartesianIndex{D})
	bool = true
	for d in 1:D
		bool &= relevant_cell_1D(k, level1[d]-1, cell1[d], level2[d]-1, cell2[d])
	end
	return bool
end

# Takes an input (level, cell, mode) and yielding the
# corresponding index in the vector of coefficients for the nodal/point basis
function get_point_index(k::Int, l::Int, c::Int, m::Int)
	if l == 0
		return m
	end
	return k + (k-1) * (1<<(l-1) - 1) + (k-1) * (c-1) + m
end

function get_modal_index(k::Int, l::Int, c::Int, m::Int)
	if l == 0
		return m
	end
	return k * (1<<(l-1) + (c-1) ) + m
end

function inner_loop_to_points{D, T<:Real}(c2::Int, k::Int, l::CartesianIndex{D},
									 level::CartesianIndex{D},
									 cell::CartesianIndex{D},
									 mode::CartesianIndex{D},
									 I::Array{Int, 1},
									 V::Array{T, 1},
									 base_indices::Array{Int, 1},
									 point_mat_1D::Array{T, 2})
	cells::NTuple{D, Int} = ntuple(i -> 1<<max(0, l[i]-2), D)
	modes::NTuple{D, Int} = ntuple(i->l[i]==1?k:k-1,D)
	cellrange = CartesianRange(cells)
	moderange = CartesianRange(modes)
	for c in cellrange
		!relevant_cell(k, l, c, level, cell) && (c2 += prod(modes); continue)

		for m in moderange
			val = one(T)
			for d in 1:D
				index = get_point_index(k, l[d]-1, c[d], m[d])
				val *= point_mat_1D[index, base_indices[d]]
			end
			(abs(val) < eps(T)) && (c2 += 1; continue)
			push!(I, c2)
			push!(V, val)
			c2 += 1
		end
	end
	return c2
end

function inner_loop_to_modal{D, T<:Real}(c2::Int, k::Int, l::CartesianIndex{D},
									 level::CartesianIndex{D},
									 cell::CartesianIndex{D},
									 mode::CartesianIndex{D},
									 I::Array{Int, 1},
									 V::Array{T, 1},
									 base_indices::Array{Int, 1},
									 point_mat_1D::Array{T, 2})
	cells::NTuple{D, Int} = ntuple(i -> 1<<max(0, l[i]-2), D)
	modes::NTuple{D, Int} = ntuple(i -> k, D)
	cellrange = CartesianRange(cells)
	moderange = CartesianRange(modes)
	for c in cellrange
		!relevant_cell(k, l, c, level, cell) && (c2 += prod(modes); continue)

		for m in moderange
			val = one(T)
			for d in 1:D
				index = get_modal_index(k, l[d]-1, c[d], m[d])
				val *= point_mat_1D[index, base_indices[d]]
			end
			(abs(val) < eps(T)) && (c2 += 1; continue)
			push!(I, c2)
			push!(V, val)
			c2 += 1
		end
	end
	return c2
end

# Evaluates a given multidimensional basis element in the modal (DG Sparse) basis
# on all the points of the corresponding sparse grid
function eval_points{D, T<:Real}(k::Int, n::Int,
								level::CartesianIndex{D},
								cell::CartesianIndex{D},
								mode::CartesianIndex{D},
								modal_indices::Array{Int, 1},
								point_mat_1D::Array{T, 2};
								scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	ls::NTuple{D, Int} = ntuple(i->(n+1),D)
	I  = Int[]
	V  = T[]
	c2 = 1

	for l in CartesianRange(ls)
		cutoff(l) && continue

		c2 = inner_loop_to_points(c2, k, l, level, cell, mode, I, V, modal_indices, point_mat_1D)
	end
	return sparsevec(I, V)
end

# Using the above method, builds the transformation matrix
# From the modal sparse basis to
# the "point basis" of values of a given function at the points on the sparse grid
function modal2points{D}(k::Int, n::Int,
							srefVD::Array{NTuple{3,CartesianIndex{D}},1};
							scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	len = get_size(D, k, n; scheme=scheme)
	modal_ref = D2Vref(1, k, n)
	point_mat_1D = full(modal2points_1D(k, n))
	I = Int[]
	J = Int[]
	V = Real[]
	for c1 in 1:len
		lcm   = srefVD[c1]
		level = lcm[1]
		cell  = lcm[2]
		mode  = lcm[3]
		modal_indices = [modal_ref[(CartesianIndex{1}(level[d],),
								  CartesianIndex{1}(cell[d],),
								  CartesianIndex{1}(mode[d],))] for d in 1:D]
		basis_eval = eval_points(k, n, level, cell, mode,
									modal_indices,
									point_mat_1D;
									scheme=scheme)
		for (c2, val) in zip(findnz(basis_eval)...)
			push!(I, c2)
			push!(J, c1)
			push!(V, val)
		end
	end
	return sparse(I, J, V)
end

# Version of the above function with simpler arguments
function modal2points(D::Int, k::Int, n::Int; scheme="sparse")
	VD = V2Dref(D, k, n; scheme=scheme)
	return modal2points(k, n, VD; scheme=scheme)
end

# Given a specific collocation point, this constructs the corresponding
# combination of nodal basis functions that vanish on all collcation points but that one
function eval_nodal{D, T<:Real}(k::Int, n::Int, level::CartesianIndex{D},
								cell::CartesianIndex{D},
								mode::CartesianIndex{D},
								inv_mat_1D::Array{T, 2};
								scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	ls::NTuple{D, Int} = ntuple(i->(n+1),D)
	point_indices = [get_point_index(k, level[d]-1, cell[d], mode[d]) for d in 1:D]

	I = Int[]
	V = T[]
	c2 = 1
	for l in CartesianRange(ls)
		cutoff(l) && continue

		c2 = inner_loop_to_points(c2, k, l, level, cell, mode, I, V, point_indices, inv_mat_1D)
	end
	return sparsevec(I, V)
end

# Using the above, this gives the transformation to go from the point basis
# to the nodal basis of functions
function points2nodal(D::Int, k::Int, n::Int; scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	ls	   = ntuple(i->(n+1),D)
	hier_ref = D2Vref(1, k, n)
	inv_mat_1D = full(points2nodal_1D(k, n))
	I = Int[]
	J = Int[]
	V = Real[]

	c1 = 1
	for l in CartesianRange(ls)
		cutoff(l) && continue

		cells = ntuple(i -> 1<<max(0, l[i]-2), D)
		modes = ntuple(i->l[i]==1?k:(k-1),D)
		for c in CartesianRange(cells)
			for m in CartesianRange(modes)
				coeffs = eval_nodal(k, n, l, c, m,
									inv_mat_1D;
									scheme=scheme)
				for (c2, val) in zip(findnz(coeffs)...)
					push!(I, c2)
					push!(J, c1)
					push!(V, val)
				end
				c1 += 1
			end
		end
	end
	return dropzeros!(sparse(I, J, V))
end

# Given a nodal basis function, this expresses it in terms of the (sparse) DG Basis
function eval_DG{D, T<:Real}(k::Int, n::Int,
								level::CartesianIndex{D},
								cell::CartesianIndex{D},
								mode::CartesianIndex{D},
								hier_ref::Dict{NTuple{3, CartesianIndex{1}}, Int},
								inv_nodal_mat_1D::Array{T, 2};
								scheme="sparse")
	nodal_indices = [get_point_index(k, level[d]-1, cell[d], mode[d]) for d in 1:D]
	cutoff = get_cutoff(scheme, D, n)
	#len = get_size(D, k, n; scheme=scheme)
	I = Int[]
	V = T[]
	ls::NTuple{D, Int} = ntuple(i->(n+1),D)

	c2 = 1
	for l in CartesianRange(ls)
		cutoff(l) && continue

		c2 = inner_loop_to_modal(c2, k, l, level, cell, mode, I, V, nodal_indices, inv_nodal_mat_1D)
	end
	return sparsevec(I, V)
end

# This gives the transformation matrix between the (sparse) nodal basis
# and the (sparse) DG basis already developed in the previous paper
function nodal2modal{D}(k::Int, n::Int,
							srefVD::Array{NTuple{3, CartesianIndex{D}},1};
							scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	ls     = ntuple(i->(n+1), D)
	hier_ref = D2Vref(1, k, n)
	inv_point_mat_1D = full(nodal2modal_1D(k, n))
	I = Int[]
	J = Int[]
	V = Real[]

	c1 = 1
	for l in CartesianRange(ls)
		cutoff(l) && continue

		cells = ntuple(i -> 1<<max(0, l[i]-2), D)
		modes = ntuple(i->l[i]==1?k:k-1,D)
		for c in CartesianRange(cells)
			for m in CartesianRange(modes)
				coeffs = eval_DG(k, n, l, c, m,
									hier_ref,
									# srefVD,
									inv_point_mat_1D;
									scheme=scheme)
				for (c2, val) in zip(findnz(coeffs)...)
					push!(I, c2)
					push!(J, c1)
					push!(V, val)
				end
				c1 += 1
			end
		end
	end

	return dropzeros!(sparse(I, J, V))
end

# Version of the above function with simpler arguments
function nodal2modal(D::Int, k::Int, n::Int; scheme="sparse")
	VD = V2Dref(D, k, n; scheme=scheme)
	return nodal2modal(k, n, VD; scheme=scheme)
end

function points2modal(D::Int, k::Int, n::Int)
	return nodal2modal(D, k, n)*points2nodal(D, k, n)
end


text:
##
# Based on Roberto Garrappa's codes
function JDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N, METH, METH_fft, Probl)

    # I defined them bc they couldn't sum themselves in a loop (e.g. nxi += r)
    nxi = 0
    nxf = 0
    nyi = 0
    nyf = 0

    nxi = nx0
    nxf = nx0 + L * r - 1
    nyi = ny0
    nyf = ny0 + L * r - 1
    is = 1

    s_nxi = Int64.(zeros(L))
    s_nxf = Int64.(zeros(L))
    s_nyi = Int64.(zeros(L))
    s_nyf = Int64.(zeros(L))

    s_nxi[is] = nxi
    s_nxf[is] = nxf
    s_nyi[is] = nyi
    s_nyf[is] = nyf

    i_triangolo = 0
    stop = 0
    stop = false # We have to define stop (e.g. by zeros) before the loop to keep changes inside the loop!

    while !stop

        stop = ((nxi + r - 1 == nx0 + L * r - 1) || (nxi + r - 1 >= Nr - 1)) # It stops when current triangle ends at the end of the square

        zn = JQuadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, METH_fft, Probl)

        y, fy = JTriangolo(nxi, nxi + r - 1, t, y, fy, zn, N, METH, Probl)
        i_triangolo += 1

        if !stop

            if nxi + r - 1 == nxf   # The triangle ends where the square ends -> you level down

                i_Delta = Int64.(ff[i_triangolo])
                Delta = i_Delta * r

                nxi = s_nxf[is] + 1
                nxf = s_nxf[is]  + Delta
                nyi = s_nxf[is] - Delta + 1
                nyf = s_nxf[is]

                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf

            else # The triangle ends before the square -> a square is made next to it

                nxi += r
                nxf = nxi + r - 1
                nyi = nyf + 1
                nyf += r

                is += 1

                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf

            end

        end

    end

    return y, fy

end

##
# Based on Roberto Garrappa's codes
function JQuadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, METH_fft, Probl)

    coef_end = nxf - nyi + 1
    i_fft = Int64.(log2(coef_end / METH.r))
    funz_beg = nyi + 1
    funz_end = nyf + 1
    Nnxf = min(N, nxf)

    vett_funz = zeros(Probl.problem_size, coef_end)

    vett_funz[:, 1:funz_end-funz_beg+1] = ifelse(nyi == 0,
    [zeros(Probl.problem_size, 1) fy[:, funz_beg + 1:funz_end]], # Evaluation of the lowest square
    fy[:, funz_beg:funz_end])                                    # we have this "else" only in Jacob

    vett_funz_fft = fft(vett_funz, 2)
    zzn = zeros(Probl.problem_size, coef_end)

    for i in 1:Probl.problem_size

        i_β = min(Probl.β_length, i)
        Z = METH_fft.an_fft[i_β, METH_fft.index_fft[1, i_fft]:METH_fft.index_fft[2, i_fft]] .* vett_funz_fft[i, :]
        zzn[i, :] = real(ifft(Z))

    end

    zzn = zzn[:, nxf - nyf + 1:end]
    zn[:, nxi + 1:Nnxf + 1] += zzn[:, 1:Nnxf - nxi + 1]

    return zn

end

##
# Based on Roberto Garrappa's codes
function JTriangolo(nxi, nxf, t, y, fy, zn, N, METH, Probl)

    for n in nxi:min(N, nxf)

        # Evaluation of the predictor
        Φ = zeros(Probl.problem_size, 1)

        for j in nxi:n - 1

            Φ += METH.an[1:Probl.β_length,n - j + 1] .* fy[:, j + 1]

        end

        St = taylor_expansion(t[n + 1], Probl.ic)

        Φ_n = St + METH.hα2 .* (METH.a0[1:Probl.β_length, n + 1] .* fy[:, 1] .+ zn[:, n + 1] .+ Φ)
        yn0 = y[:, n]
        fn0 = f_value(Probl.f_fun(t[n + 1], yn0, Probl.param...), Probl.problem_size)
        Jfn0 = Probl.JF(t[n + 1], yn0, Probl.param...)
        Gn0 = yn0 - METH.hα2 .* fn0 - Φ_n

        stop = 0
        stop = false
        mu_it = 0

        yn1 = zeros(Probl.problem_size, 1)
        fn1 = zeros(Probl.problem_size, 1)

        while !stop

            JGn0 = I(Probl.problem_size) .- Diagonal(METH.hα2) * Jfn0
            yn1 = yn0 - inv(JGn0) * Gn0
            fn1 = f_value(Probl.f_fun(t[n + 1], yn1, Probl.param...), Probl.problem_size)
            Gn1 = yn1 - METH.hα2 .* fn1 - Φ_n

            mu_it += 1

            if METH.StopIt == "Convergence"

                stop = (norm(yn1 - yn0, Inf) < METH.μTol || norm(Gn1, Inf) <  METH.μTol)

                if (mu_it > METH.itmax && !stop)

                    stop = true

                end

            else

                stop = (mu_it == METH.μ)

            end

            yn0 = yn1
            Gn0 = Gn1

            if !stop

                Jfn0 = Probl.JF(t[n + 1], yn0, Probl.param...)

            end

        end

        y[:, n + 1] = yn1
        fy[:, n + 1] = fn1

    end

    return y, fy

end


text:
module ABM4OSN

using LightGraphs
using DataFrames
using Statistics
using Random
using StatsBase

include("config.jl")
include("tweet.jl")
include("agent.jl")
include("network.jl")
include("simulation.jl")

include("add_input.jl")
include("create_agents.jl")
include("drop_input.jl")
include("generate_agent_props.jl")
include("like.jl")
include("log_network.jl")
include("publish_tweet.jl")
include("retweet.jl")
include("set_inactive.jl")
include("update_check_regularity.jl")
include("update_feed.jl")
include("update_opinion.jl")
include("update_perceiv_publ_opinion.jl")

greet() = print("Hello World!")

end # module

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
abstract AbstractGrid{D}

type Uniform{D,S} <: AbstractGrid{D} end
Uniform{D,T<:Integer}(x::NTuple{D,T}) = Uniform{D,Val{x}}
Uniform(x::Int) = Uniform{1,Val{tuple(x)}}

type Irregular{D,N} <: AbstractGrid{D} end
Irregular(D::Integer,N::Integer) = Irregular{D,N}
Irregular(N::Integer) = Irregular{1,N}


text:
cd(@__DIR__); include("setups/simpletree.jl");
gr(dpi=200)

plt = plot(0:(N - 1), 𝛌, c = :black, lw = 1, ylims = [0, 5], legend = false,
        shape = :circle, frame = :box, xlab = latexstring("l"), ms = 4,
        ylab = latexstring("\\lambda_l"), size = (500, 400))
savefig(plt, "../figs/SimpleTree_eigenvalues.png")