text:
const modificationtests = Dict{String, Function}()

"""
    set_function_single_variable(model::MOI.ModelLike, config::TestConfig)

Test that modifying the function of a `SingleVariable`-in-`LessThan` constraint
throws a [`SettingSingleVariableFunctionNotAllowed`](@ref) error.
"""
function set_function_single_variable(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x, y
        maxobjective: 1.0x + 1.0y
        c: x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    # We test this after the creation of every `SingleVariable` constraint
    # to ensure a good coverage of corner cases.
    @test c.value == x.value
    err = MOI.SettingSingleVariableFunctionNotAllowed()
    func = MOI.SingleVariable(y)
    @test_throws err MOI.set(model, MOI.ConstraintFunction(), c, func)
end
modificationtests["set_function_single_variable"] = set_function_single_variable

"""
    solve_set_singlevariable_lessthan(model::MOI.ModelLike, config::TestConfig)

Test set modification SingleVariable-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_set_singlevariable_lessthan(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c: x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    @test c.value == x.value
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -1.0)]
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(2.0)
    test_model_solution(model, config;
        objective_value   = 2.0,
        variable_primal   = [(x, 2.0)],
        constraint_primal = [(c, 2.0)],
        constraint_dual   = [(c, -1.0)]
    )
end
modificationtests["solve_set_singlevariable_lessthan"] = solve_set_singlevariable_lessthan

"""
    solve_transform_singlevariable_lessthan(model::MOI.ModelLike, config::TestConfig)

Test set transformation of a SingleVariable-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_transform_singlevariable_lessthan(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c: x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    @test c.value == x.value
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -1.0)]
    )
    c2 = MOI.transform(model, c, MOI.GreaterThan(2.0))
    @test !MOI.is_valid(model, c)
    @test MOI.is_valid(model, c2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    test_model_solution(model, config;
        objective_value   = 2.0,
        variable_primal   = [(x, 2.0)],
        constraint_primal = [(c2, 2.0)],
        constraint_dual   = [(c2, 1.0)]
    )
end
modificationtests["solve_transform_singlevariable_lessthan"] = solve_transform_singlevariable_lessthan

"""
    solve_set_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)

Test modifying set of ScalarAffineFunction-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_set_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c: 1.0x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -1.0)]
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(2.0)
    test_model_solution(model, config;
        objective_value   = 2.0,
        variable_primal   = [(x, 2.0)],
        constraint_primal = [(c, 2.0)],
        constraint_dual   = [(c, -1.0)]
    )
end
modificationtests["solve_set_scalaraffine_lessthan"] = solve_set_scalaraffine_lessthan

"""
    solve_coef_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)

Test modifying a variable coefficient in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_coef_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c: 1.0x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -1.0)]
    )
    MOI.modify(model, c, MOI.ScalarCoefficientChange(x, 2.0))
    test_model_solution(model, config;
        objective_value   = 0.5,
        variable_primal   = [(x, 0.5)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -0.5)]
    )
end
modificationtests["solve_coef_scalaraffine_lessthan"] = solve_coef_scalaraffine_lessthan

"""
    solve_func_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)

Test setting the function in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_func_scalaraffine_lessthan(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c: 1.0x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -1.0)]
    )
    MOI.set(model, MOI.ConstraintFunction(), c,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 0.0)
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 0.0)
    test_model_solution(model, config;
        objective_value   = 0.5,
        variable_primal   = [(x, 0.5)],
        # constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, -0.5)]
    )
end
modificationtests["solve_func_scalaraffine_lessthan"] = solve_func_scalaraffine_lessthan

"""
    solve_func_vectoraffine_nonneg(model::MOI.ModelLike, config::TestConfig)

Test setting the function in a VectorAffineFunction-in-Nonnegatives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_func_vectoraffine_nonneg(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x, y
        minobjective: 1.0x + 2.0y
        c: [1.0x, 2.0y] in Nonnegatives(2)
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    test_model_solution(model, config;
        objective_value   = 0.0,
        variable_primal   = [(x, 0.0), (y, 0.0)],
        constraint_primal = [(c, [0.0, 0.0])]
    )
    MOI.set(model, MOI.ConstraintFunction(), c,
            MOI.VectorAffineFunction([
                  MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                  MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y))
                ],
                [-1.0, -1.5]
            )
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.VectorAffineFunction([
                  MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                  MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y))
                ],
                [-1.0, -1.5]
            )
    test_model_solution(model, config;
        objective_value   = 2.5,
        variable_primal   = [(x, 1.0), (y, 0.75)],
        constraint_primal = [(c, [0.0, 0.0])]
    )
end
modificationtests["solve_func_vectoraffine_nonneg"] = solve_func_vectoraffine_nonneg

"""
    solve_const_vectoraffine_nonpos(model::MOI.ModelLike, config::TestConfig)

Test modifying the constant term in a VectorAffineFunction-in-Nonpositives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function solve_const_vectoraffine_nonpos(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x, y
        maxobjective: 1.0x + 2.0y
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.add_constraint(model,
            MOI.VectorAffineFunction([
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y))
                ],
                [0.0, 0.0]
            ),
            MOI.Nonpositives(2)
        )
    test_model_solution(model, config;
        objective_value   = 0.0,
        variable_primal   = [(x, 0.0), (y, 0.0)],
        constraint_primal = [(c, [0.0, 0.0])]
    )
    MOI.modify(model, c, MOI.VectorConstantChange([-1.0, -1.5]))
    test_model_solution(model, config;
        objective_value   = 2.5,
        variable_primal   = [(x, 1.0), (y, 0.75)],
        constraint_primal = [(c, [0.0, 0.0])]
    )
end
modificationtests["solve_const_vectoraffine_nonpos"] = solve_const_vectoraffine_nonpos

"""
    solve_multirow_vectoraffine_nonpos(model::MOI.ModelLike, config::TestConfig)

Test modifying the variable coefficients in a
VectorAffineFunction-in-Nonpositives constraint. If `config.solve=true` confirm
that it solves correctly.
"""
function solve_multirow_vectoraffine_nonpos(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.add_constraint(model,
            MOI.VectorAffineFunction([
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x))
                ],
                [-1.0, -1.0]
            ),
            MOI.Nonpositives(2)
        )
    test_model_solution(model, config;
        objective_value   = 0.5,
        variable_primal   = [(x, 0.5)],
        constraint_primal = [(c, [-0.5, 0.0])]
    )
    MOI.modify(model, c, MOI.MultirowChange(x, [(1,4.0), (2,3.0)]))
    test_model_solution(model, config;
        objective_value   = 0.25,
        variable_primal   = [(x, 0.25)],
        constraint_primal = [(c, [0.0, -0.25])]
    )
end
modificationtests["solve_multirow_vectoraffine_nonpos"] = solve_multirow_vectoraffine_nonpos

"""
    solve_const_scalar_objective(model::MOI.ModelLike, config::TestConfig)

Test the constant of a scalaraffine objective. If `config.solve=true` confirm
that it solves correctly.
"""
function solve_const_scalar_objective(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x + 2.0
        c1: 1.0x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    test_model_solution(model, config;
        objective_value   = 3.0,
        variable_primal   = [(x, 1.0)]
    )
    MOI.modify(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarConstantChange(3.0)
    )
    test_model_solution(model, config;
        objective_value   = 4.0,
        variable_primal   = [(x, 1.0)]
    )
end
modificationtests["solve_const_scalar_objective"] = solve_const_scalar_objective

"""
    solve_coef_scalar_objective(model::MOI.ModelLike, config::TestConfig)

Test modifying a variable coefficient in a scalaraffine objective. If
`config.solve=true` confirm that it solves correctly.
"""
function solve_coef_scalar_objective(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    MOIU.loadfromstring!(model,"""
        variables: x
        maxobjective: 1.0x
        c1: 1.0x <= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)]
    )
    MOI.modify(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarCoefficientChange(x, 3.0)
    )
    test_model_solution(model, config;
        objective_value   = 3.0,
        variable_primal   = [(x, 1.0)]
    )
end
modificationtests["solve_coef_scalar_objective"] = solve_coef_scalar_objective

function delete_variable_with_single_variable_obj(model::MOI.ModelLike,
                                                  config::TestConfig)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(model,"""
        variables: x, y
        minobjective: x
        c: x >= 1.0
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.get(model, MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}, "c")
    @test c.value == x.value
    MOI.delete(model, y)
    test_model_solution(model, config;
        objective_value   = 1.0,
        variable_primal   = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual   = [(c, 1.0)]
    )
end
modificationtests["delete_variable_with_single_variable_obj"] = delete_variable_with_single_variable_obj

"""
    delete_variables_in_a_batch(model::MOI.ModelLike, config::TestConfig)

Test deleting many variables in a batch (i.e. using the delete method which
takes a vector of variable references). If `config.solve=true` confirm that it
solves correctly.
"""
function delete_variables_in_a_batch(model::MOI.ModelLike,
                                     config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(model,"""
        variables: x, y, z
        minobjective: 1.0 * x + 2.0 * y + 3.0 * z
        c1: x >= 1.0
        c2: y >= 1.0
        c3: z >= 1.0
    """)
    x, y, z = MOI.get(model, MOI.ListOfVariableIndices())
    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test MOI.is_valid(model, z)
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 6.0 atol=atol rtol=rtol
    end
    MOI.delete(model, [x, z])
    @test !MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test !MOI.is_valid(model, z)
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
    end
end
modificationtests["delete_variables_in_a_batch"] = delete_variables_in_a_batch

@moitestset modification


text:
# This example demonstrates the use of the high-level Matlab-like interface
#
# First of all we create an array with top-of the atmosphere radiation data

using NetCDF

include("toa.jl")

# Define longitudes and latitudes, day and timesteps
lat = collect(-89:89)
lon = collect(0:359)
day = 1
tim = collect(0:23)

# Create radiation array
rad = Float64[g_pot(x2, x1, day, x3) for x1 in lon, x2 in lat, x3 in tim]

# Define some attributes of the variable (optionlal)
varatts = Dict("longname" => "Radiation at the top of the atmosphere", "units" => "W/m^2")
lonatts = Dict("longname" => "Longitude", "units" => "degrees east")
latatts = Dict("longname" => "Latitude", "units" => "degrees north")
timatts = Dict("longname" => "Time", "units" => "hours since 01-01-2000 00:00:00")

# Now we create the file radiation.nc and call the variable rad, define also the dimensions the variables depends on

isfile("radiation.nc") && rm("radiation.nc")
nccreate(
    "radiation.nc",
    "rad",
    "lon",
    lon,
    lonatts,
    "lat",
    lat,
    latatts,
    "time",
    tim,
    timatts,
    atts = varatts,
)

# Now we can write values to the file

ncwrite(rad, "radiation.nc", "rad")

# Reading the whole dataset is done by:

x = ncread("radiation.nc", "rad")
println("Successfully read an array of size ", size(x))

# Additional
# Reading parts of a file, for example reading only time steps 5 and 6 can be done with:

x = ncread("radiation.nc", "rad", [1, 1, 5], [-1, -1, 2])
println("Successfully read an array of size ", size(x))

# here the first array [1,1,5] gives the starting position for reading while the second array [1,1,2] gives the number of blocks to be read along each dimension.
# If the count is set to -1 the whole dimension is read.

# The same can be done for ncwrite

x = ones(Float64, 1, 360, 1)
# ncwrite(x,"radiation.nc","rad",[60,1,4])
# This would set all values at 60deg longitude at the 4th time step to 1.

ext: jl


text:
import GeometryBasics

"""
    struct Mesh{N, T, PT <: Polyhedron{T}} <: GeometryBasics.GeometryPrimitive{N, T}
        polyhedron::PT
        coordinates::Union{Nothing, Vector{GeometryBasics.Point{3, T}}}
        faces::Union{Nothing, Vector{GeometryBasics.TriangleFace{Int}}}
        normals::Union{Nothing, Vector{GeometryBasics.Point{3, T}}}
    end

Mesh wrapper type that inherits from `GeometryPrimitive` to be used for plotting
a polyhedron. Note that `Mesh(p)` is type unstable but one can use `Mesh{3}(p)`
instead if it is known that `p` is defined in a 3-dimensional space.
"""
mutable struct Mesh{N, T, PT <: Polyhedron{T}} <: GeometryBasics.GeometryPrimitive{N, T}
    polyhedron::PT
    coordinates::Union{Nothing, Vector{GeometryBasics.Point{N, T}}}
    faces::Union{Nothing, Vector{GeometryBasics.TriangleFace{Int}}}
    normals::Union{Nothing, Vector{GeometryBasics.Point{N, T}}}
end
function Mesh{N}(polyhedron::Polyhedron{T}) where {N, T}
    return Mesh{N, T, typeof(polyhedron)}(polyhedron, nothing, nothing, nothing)
end
function Mesh(polyhedron::Polyhedron, ::StaticArrays.Size{N}) where N
    return Mesh{N[1]}(polyhedron)
end
function Mesh(polyhedron::Polyhedron, N::Int)
    # This is type unstable but there is no way around that,
    # use polyhedron built from StaticArrays vector to avoid that.
    return Mesh{N}(polyhedron)
end
function Mesh(polyhedron::Polyhedron)
    return Mesh(polyhedron, FullDim(polyhedron))
end

function fulldecompose!(mesh::Mesh)
    if mesh.coordinates === nothing
        mesh.coordinates, mesh.faces, mesh.normals = fulldecompose(mesh)
    end
    return
end

# Creates a scene for the vizualisation to be used to truncate the lines and rays
function scene(vr::VRep, ::Type{T}) where T
    # First compute the smallest rectangle containing the P-representation (i.e. the points).
    (xmin, xmax) = extrema(map((x)->x[1], points(vr)))
    (ymin, ymax) = extrema(map((x)->x[2], points(vr)))
    (zmin, zmax) = extrema(map((x)->x[3], points(vr)))
    width = max(xmax-xmin, ymax-ymin, zmax-zmin)
    if width == zero(T)
        width = 2
    end
    scene = GeometryBasics.HyperRectangle{3, T}([(xmin + xmax) / 2 - width,
                                                 (ymin + ymax) / 2 - width,
                                                 (zmin + zmax) / 2 - width],
                                                2 * width * ones(T, 3))
    # Intersection of rays with the limits of the scene
    (start, r) -> begin
        ray = coord(r)
        λ = nothing
        min_scene = minimum(scene)
        max_scene = maximum(scene)
        for i in 1:3
            r = ray[i]
            if !iszero(r)
                cur = max((min_scene[i] - start[i]) / r, (max_scene[i] - start[i]) / r)
                if λ === nothing || cur < λ
                    λ = cur
                end
            end
        end
        start + λ * ray
    end
end

function _isdup(zray, triangles)
    for tri in triangles
        normal = tri[2]
        if isapproxzero(cross(zray, normal)) && dot(zray, normal) > 0 # If A[j,:] is almost 0, it is always true...
            # parallel and equality or inequality and same sense
            return true
        end
    end
    false
end
_isdup(poly, hidx, triangles) = _isdup(get(poly, hidx).a, triangles)

function fulldecompose(poly_geom::Mesh{3}, ::Type{T}) where T
    poly = poly_geom.polyhedron
    exit_point = scene(poly, T)

    triangles = Tuple{Tuple{Vector{T},Vector{T},Vector{T}}, Vector{T}}[]

    function decomposeplane(hidx)
        h = get(poly, hidx)
        # xray should be the rightmost ray
        xray = nothing
        # xray should be the leftmost ray
        yray = nothing
        zray = h.a
        isapproxzero(zray) && return

        # Checking rays
        counterclockwise(a, b) = dot(cross(a, b), zray)
        face_vert = pointtype(poly)[]
        for x in points(poly)
            if _isapprox(dot(x, zray), h.β)
                push!(face_vert, x)
            end
        end
        hull, lines, rays = _planar_hull(3, face_vert, incidentlines(poly, hidx), incidentrays(poly, hidx), counterclockwise, r -> cross(zray, r))
        if isempty(lines)
            if length(hull) + length(rays) < 3
                return
            end
            @assert length(rays) <= 2
            if !isempty(rays)
                if length(rays) + length(hull) >= 2
                    push!(hull, exit_point(last(hull), last(rays)))
                end
                push!(hull, exit_point(first(hull), first(rays)))
            end
        else
            if length(hull) == 2
                @assert length(lines) == 1 && isempty(rays)
                a, b = hull
                line = first(lines)
                empty!(hull)
                push!(hull, exit_point(a, line))
                push!(hull, exit_point(a, -line))
                push!(hull, exit_point(b, -line))
                push!(hull, exit_point(b, line))
            else
                @assert length(hull) == 1 && isempty(rays)
                center = first(hull)
                empty!(hull)
                a = first(lines)
                b = nothing
                if length(lines) == 2
                    @assert isempty(rays)
                    b = last(lines)
                elseif !isempty(rays)
                    @assert length(lines) == 1
                    @assert length(rays) == 1
                    b = linearize(first(rays))
                end
                push!(hull, exit_point(center, a))
                if b !== nothing
                    push!(hull, exit_point(center, b))
                end
                push!(hull, exit_point(center, -a))
                if b !== nothing && length(lines) == 2 || length(rays) >= 2
                    @assert length(rays) == 2
                    push!(hull, exit_point(center, -b))
                end
            end
        end

        if length(hull) >= 3
            a = pop!(hull)
            b = pop!(hull)
            while !isempty(hull)
                c = pop!(hull)
                push!(triangles, ((a, b, c), zray))
                b = c
            end
        end
    end

    for hidx in eachindex(hyperplanes(poly))
        decomposeplane(hidx)
    end
    # If there is already a triangle, his normal is an hyperplane and it is the only face
    if isempty(triangles)
        for hidx in eachindex(halfspaces(poly))
            if !_isdup(poly, hidx, triangles)
                decomposeplane(hidx)
            end
        end
    end

    ntri = length(triangles)
    pts = Vector{GeometryBasics.Point{3, T}}(undef, 3ntri)
    faces = Vector{GeometryBasics.TriangleFace{Int}}(undef, ntri)
    ns = Vector{GeometryBasics.Point{3, T}}(undef, 3ntri)
    for i in 1:ntri
        tri = pop!(triangles)
        normal = tri[2]
        for j = 1:3
            idx = 3*(i-1)+j
            #ns[idx] = -normal
            ns[idx] = normal
        end
        faces[i] = collect(3*(i-1) .+ (1:3))
        k = 1
        for k = 1:3
            # reverse order of the 3 vertices so that if I compute the
            # normals with the `normals` function, they are in the good
            # sense.
            # I know I don't use the `normals` function but I don't know
            # what is the OpenGL convention so I don't know if it cares
            # about the order of the vertices.
            pts[3*i-k+1] = tri[1][k]
        end
    end
    # If the type of ns is Rational, it also works.
    # The normalized array in in float but then it it recast into Rational
    map!(normalize, ns, ns)
    (pts, faces, ns)
end

fulldecompose(poly::Mesh{N, T}) where {N, T} = fulldecompose(poly, typeof(one(T)/2))

GeometryBasics.coordinates(poly::Mesh) = (fulldecompose!(poly); poly.coordinates)
GeometryBasics.faces(poly::Mesh) = (fulldecompose!(poly); poly.faces)
GeometryBasics.texturecoordinates(poly::Mesh) = nothing
GeometryBasics.normals(poly::Mesh) = (fulldecompose!(poly); poly.normals)


text:
#Demo 6: Fourier Transforms
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
dbvsf = cons(:a, labels = set(yaxis="Amplitude (dB)", xaxis="Frequency (Hz)"))

#==Input data
===============================================================================#
tfund = 1e-9 #Fundamental
osr = 20 #(fundamental) oversampling ratio

#==Computations
===============================================================================#
tvec = timespace(:ts, tfund/osr, :tfund, tfund)
Δt = step(tvec)
sigvec = rand(DataFloat, length(tvec))

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("padding", collect(0:10:100))
]

#Generate data:
sigSpec = fill(DataHR, sweeplist) do padding
	#t-vector (+padding):
	tpad = 0:Δt:(Δt*((length(tvec)-1)+padding))
	#x(t)-vector:
	xtpad = zeros(length(tpad))
	xtpad[1:length(sigvec)] = sigvec
	#Padded time-domain signal:
	pad = DataTime(tpad, xtpad, tperiodic=false)

	Fpad = freqdomain(pad) #Frequency-domain signal (DataFreq)
	return abs2(fspectrum(Fpad)) #Spectrum (DataF1)
end

#Informative:
fmax = maximum(xval(sigSpec))

fstep = DataHR{DataFloat}(sweeplist) #Create empty data
for inds in subscripts(fstep)
	fstep.elem[inds...] = sigSpec.elem[inds...].x[2]
end

@show fmax
@show fstep

#==Generate plot
===============================================================================#
plot_t = push!(cons(:plot, vvst, title="Time Domain"),
	cons(:wfrm, DataF1(collect(tvec), sigvec)),
)
plot_f = push!(cons(:plot, dbvsf, title="Sampled Frequency Spectrum"),
	cons(:wfrm, dB20(sigSpec)),
)

pcoll = push!(cons(:plotcoll, title="Fourier Transform vs Padding"), plot_t, plot_f)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using LoggingExtras
using Test
using Test: collect_test_logs, TestLogger

using Base.CoreLogging
using Base.CoreLogging: BelowMinLevel, Debug, Info, Warn, Error

@testset "Tee" begin
    @testset "contructor" begin
        @testset "mixed types" begin
            @test TeeLogger(TestLogger(), NullLogger()) isa TeeLogger
        end

        @testset "errors if given nonloggers" begin
            @test_throws Exception TeeLogger(stdout, stderr)
        end
    end
    @testset "basic use with compositional levels" begin
        testlogger_info = TestLogger(min_level=Info)
        testlogger_warn = TestLogger(min_level=Warn)

        with_logger(TeeLogger(testlogger_warn, testlogger_info)) do
            @info "info1"
            @warn "warn1"
            @info "info2"
        end
        @test length(testlogger_info.logs) == 3
        @test length(testlogger_warn.logs) == 1
    end
end

@testset "File" begin
    mktempdir() do dir
        filepath = joinpath(dir, "log")
        with_logger(FileLogger(filepath)) do
            @info "first"
            @warn "second"
            @info "third"
        end
        logtext = String(read(filepath))
        @test occursin("first", logtext)
        @test occursin("second", logtext)
        @test occursin("third", logtext)
    end
end

@testset "Active Filter" begin
    testlogger = TestLogger()
    yodawg_filter(logargs) = startswith(logargs.message, "Yo Dawg!")

    filtered_logger = ActiveFilteredLogger(yodawg_filter, testlogger)

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
    end
    @test length(testlogger.logs) == 2
end

@testset "Early Filter" begin
    testlogger = TestLogger()
    filtered_logger = EarlyFilteredLogger(testlogger) do logargs
        logargs.level == Info  # Only exactly Info, nothing more nothing less
    end

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
    end
    @test length(testlogger.logs) == 3
end

@testset "MinLevel Filter" begin
    testlogger = TestLogger()
    filtered_logger = MinLevelLogger(testlogger, Warn)

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
        @error "MISTAKES WERE MADE"
    end
    @test length(testlogger.logs) == 2
end

@testset "Transformer" begin
    testlogger = TestLogger(min_level=Error)
    transformer_logger = TransformerLogger(testlogger) do log_msg
        # We are going to transform all warnings into errors
        if log_msg.level == Warn
            return merge(log_msg, (; level=Error))
        else
            return log_msg
        end
    end

    with_logger(transformer_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
        @error "MISTAKES WERE MADE"
    end
    @test length(testlogger.logs) == 2
end

@testset "Deprecations" begin
    testlogger = TestLogger(min_level=BelowMinLevel)

    demux_logger = DemuxLogger(testlogger)
    @test demux_logger isa TeeLogger
    @test Set(demux_logger.loggers) == Set([testlogger, global_logger()])

    demux_logger = DemuxLogger(testlogger; include_current_global=true)
    @test demux_logger isa TeeLogger
    @test Set(demux_logger.loggers) == Set([testlogger, global_logger()])

    demux_logger = DemuxLogger(testlogger; include_current_global=false)
    @test demux_logger isa TeeLogger
    @test Set(demux_logger.loggers) == Set([testlogger])
end


text:
"""
    showfeasible

Collects a matrix of only feasible outcomes given the parameter space and the
constraints. Use `collect` to output a matrix of all possible matches for
parameter space.

### Examples
```julia
julia> myset = logicalparse("a, b, c in 1:3")
a,b,c in 1:3             feasible outcomes 27 ✓          :3 3 3

julia> myset = logicalparse("a == b; c > b", myset)
a == b                   feasible outcomes 9 ✓           :1 1 3
c > b                    feasible outcomes 3 ✓           :2 2 3

julia> showfeasible(myset)
3×3 Array{Int64,2}:
 1  1  2
 1  1  3
 2  2  3
```
"""
showfeasible(x::LogicalCombo) = x[:,:,:]


text:
"""
    pl = plot_results(params, pcof; [casename = "test", savefiles = false, samplerate = 32])

Create array of plot objects that can be visualized by, e.g., `display(pl[1])`.

# Arguments
- `params::objparams`: Object holding problem definition
- `pcof::Array{Float64,1}`: Parameter vector
- `casename::String`: Default: `"test"`. String used in plot titles and in file names
- `savefiles::Bool`: Default: `false`.Set to `true` to save plots on files with automatically generated filenames
- `samplerate:: Int64`: Default: `32` samples per unit time (ns). Sample rate for generating plots.
"""
function plot_results(params::objparams, pcof::Array{Float64,1}; casename::String = "test", savefiles::Bool = false, samplerate:: Int64 = 32)
    # Set default font sizes
    fnt = Plots.font("Helvetica", 12)
    lfnt = Plots.font("Helvetica", 10)
    Plots.default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=lfnt, linewidth=1, size=(650, 350))

    nCoeff = length(pcof)
    # Is there a better approach that avoids re-allocating the working_arrays object?
    wa = Juqbox.Working_Arrays(params, nCoeff)

    custom = 0

    # data file names
    labbasename = casename * "-ctrl-lab"
    ctrlbasename = casename * "-ctrl-rot"

    println("Tikhonov coefficient: tik0 = ", params.tik0)
    println("B-carrier basis: ", params.use_bcarrier)
    println("Number of time steps: ", params.nsteps)

    # evaluate fidelity and unitaryhistory
    objv, unitaryhistory, fidelity = Juqbox.traceobjgrad(pcof, params, wa, true, false);

    # save convergence history
    convname = ""
    if savefiles
        convname = casename * "-conv" * ".png"
    end
    pconv = Juqbox.plot_conv_hist(params, convname)

    # scatter plot of control parameters
    tstring = casename * "-control-vector"
    plcof = scatter(pcof, lab="", title=tstring, xlabel="Index", ylabel="rad/ns")

    guardlev = Juqbox.identify_guard_levels(params, custom)
    forbiddenlev = Juqbox.identify_forbidden_levels(params, custom)

    # make plots of the evolution of probabilities
    pl1 = Juqbox.plotunitary(unitaryhistory, params, guardlev)
    pl3 = Juqbox.plotspecified(unitaryhistory, params, guardlev, forbiddenlev)

    if savefiles
        # save the figure with the state probabilities
        probname = casename * "-prob" * ".png"
        Plots.savefig(pl1,probname)
        println("Saved state population plot on file '", probname, "'");

        # save the figure with the forbidden state
        forbname = casename * "-forb" * ".png"
        Plots.savefig(pl3,forbname)
        println("Saved forbidden state population plot on file '", forbname, "'");
    end

    if params.Nosc == 3 # Generalize to Nosc = 2
        # evaluate marginalized probabilities
        mp = Juqbox.marginalize3(params, unitaryhistory);
        # plot them
        nsteps1=size(unitaryhistory,3)
        tm = range(0, stop = params.T, length = nsteps1);
        # one subfigure for each initial condition
        plotarray = Array{Plots.Plot}(undef, params.N) #empty array for separate plots

        for col in 1:params.N # One plot for each initial condition
            #  params.Nosc == 3
            # Example: Ne[3] = 2, Ne[2] = 3, Ne[1]=4
            #     0,    1,     2,    3,    4,     5,     6,    7,     8,    9,   10,   11,   12,   13,   14,  15,   16,   17,   18,  19,   20,   21,   22,   23  (col-1)
            # 000, 001, 002, 003, 010, 011, 012, 013, 020, 021, 022, 023,  100, 101, 102, 103, 110, 111, 112, 113, 120, 121, 122, 123
            s3 = div((col-1), params.Ne[1]*params.Ne[2])
            s12 = (col-1) % (params.Ne[1]*params.Ne[2])
            s2 = div(s12, params.Ne[1])
            s1 = s12 %  params.Ne[1]
            statestr = string( s3, s2, s1 )

            local titlestr = latexstring("From\\ state\\ |", statestr, "\\rangle")
#            local titlestr = raw"Evolution from state $|" * statestr * raw"\rangle$"
            h = Plots.plot(title = titlestr)

            for row in 1:params.Nt[3]
                labstr = "S-state " * string(row-1)
                Plots.plot!(tm, mp[row,col,:], lab = labstr, xlabel = "Time [ns]", ylabel="Marg. prob.", legend= :outerright)
            end # for row
            plotarray[col] = h
        end # for col
        if params.N <= 2
            plm = Plots.plot(plotarray..., layout = (params.N,1))
        else
            plm = Plots.plot(plotarray..., layout = params.N)
        end
    end

    # final unitary plotted in matrix form
    pluf = plot_final_unitary(unitaryhistory[:,:,end], params, fidelity)

    plen = plot_energy(unitaryhistory, params)

    # Evaluate the ctrl functions on this grid in time
    nplot = round(Int64, params.T*samplerate)
    # is this resolution sufficient for the lab frame ctrl functions so we can get meaningful FFTs?
    td = collect(range(0, stop = params.T, length = nplot+1))

    # Initialize storing of the lab drive
    labdrive = zeros(nplot+1)

    nFFT = length(labdrive)
    dt = td[2] - td[1]
    freq = fftshift( fftfreq(nFFT, 1.0/dt) )

    useMHz = true
    if useMHz
        scalefactor = 1000/(2*pi)
        unitStr = "MHz"
    else
        scalefactor = 1.0
        unitStr = "rad/ns"
    end

    # one subfigure for each control function
    plotarray = Array{Plots.Plot}(undef, params.Ncoupled+ params.Nunc) #empty array for separate plots
    plotarray_fft = Array{Plots.Plot}(undef, params.Ncoupled+ params.Nunc) #empty array for separate plots
    plotarray_fftlog = Array{Plots.Plot}(undef, params.Ncoupled+ params.Nunc) #empty array for separate plots
    plotarray_lab = Array{Plots.Plot}(undef, params.Ncoupled + params.Nunc) #empty array for separate plots

    println("Rotational frequencies: ", params.Rfreq)

    for q=1:params.Ncoupled+params.Nunc
        # evaluate ctrl functions for the q'th Hamiltonian
        pfunc, qfunc = Juqbox.evalctrl(params, pcof, td, q)

        pfunc = scalefactor .* pfunc
        qfunc = scalefactor .* qfunc

        pmax = maximum(abs.(pfunc))
        qmax = maximum(abs.(qfunc))
        # first plot for control function for the symmetric Hamiltonian
        local titlestr = "Rotating frame ctrl - " * string(q) * " Max-p=" *@sprintf("%.3e", pmax) * " Max-q=" *
            @sprintf("%.3e", qmax) * " " * unitStr
        plotarray[q] = Plots.plot(td, pfunc, lab=L"p(t)", title = titlestr, xlabel="Time [ns]",
                                  ylabel=unitStr, legend= :outerright)
        # add in the control function for the anti-symmetric Hamiltonian
        Plots.plot!(td, qfunc, lab=L"q(t)")

        println("Rot. frame ctrl-", q, ": Max-p(t) = ", pmax, " Max-q(t) = ", qmax, " ", unitStr)

        # Corresponding lab frame control
        omq = 2*pi*params.Rfreq[q] # FIX index of Rfreq
        labdrive .= 2*pfunc .* cos.(omq*td) .- 2*qfunc .* sin.(omq*td)

        lmax = maximum(abs.(labdrive))
        local titlestr = "Lab frame ctrl - " * string(q) * " Max=" *@sprintf("%.3e", lmax) * " " * unitStr
        plotarray_lab[q]= Plots.plot(td, labdrive, lab="", title = titlestr, size = (650, 250), xlabel="Time [ns]", ylabel=unitStr)

        println("Lab frame ctrl-", q, " Max amplitude = ", lmax, " ", unitStr)

        # plot the Fourier transform of the control function in the lab frame
        # Fourier transform
        Fdr_lab = fftshift( fft(labdrive) ) / nFFT

        local titlestr = "Spectrum, lab frame ctrl - " * string(q)
        plotarray_fft[q] = Plots.plot(freq, abs.(Fdr_lab), lab="", title = titlestr, size = (650, 350), xlabel="Frequency [GHz]",
                                      ylabel="Amplitude " * unitStr, framestyle = :box) #, grid = :hide

        fmin = 0.5*minimum(params.Rfreq)
        fmax = maximum(params.Rfreq) + 0.5
        xlims!((fmin, fmax))

        # log-scale spectrum
        mag_Fdr_lab = abs.(Fdr_lab)
        plotarray_fftlog[q] = Plots.plot(title = titlestr, xlabel="Frequency [GHz]",
                                         ylabel="Amplitude " * unitStr, yaxis=:log10, framestyle = :box)
        if minimum(mag_Fdr_lab) > 0.0
            Plots.plot!(freq, mag_Fdr_lab, lab="")
            xlims!((fmin, fmax))
        end

        if savefiles
            # Save ctrl functions on file
            pqname = ctrlbasename * "-" * string(q) * ".dat"
            writedlm(pqname, [pfunc, qfunc])
            println("Saved ctrl functions for Hamiltonian #", q, " on file '", pqname, "', samplerate = ", samplerate);

            # save the lab frame ctrl func
            labname = labbasename * "-" * string(q) * ".dat"
            writedlm(labname, labdrive)
            println("Saved lab frame ctrl function on file '", labname, "'", " samplerate = ", samplerate);
        end
    end

    # # Add in uncoupled controls
    # if(params.Nunc >  0)
    #     max_uncoupled = zeros(length(params.Hunc_ops))

    #     for q =1:params.Nunc
    #         qs = 2*params.Ncoupled + (q - 1)*2
    #         qa = qs+1
    #         pfunc = scalefactor .* Juqbox.evalctrl(params, pcof, td, qs)
    #         qfunc = scalefactor .* Juqbox.evalctrl(params, pcof, td, qa)
    #         ffunc = 2*(pfunc .* cos(2*pi*params.Rfreq[q]) .- qfunc .* sin(2*pi*params.Rfreq[q]))
    #         max_uncoupled[q] = maximum(abs.(ffunc))
    #         plotarray_lab[params.Ncoupled + q] =  Plots.plot(td, qfunc, lab="", linewidth = 2, title = "Uncoupled Ctrl Function",
    #                                                          size = (650, 250), xlabel="Time [ns]", ylabel=unitStr)
    #     end
    #     println("Max amplitude uncoupled ctrl = ", maximum(max_uncoupled), unitStr)

    #     # TODO: save uncoupled controls on file
    # end

   # Accumulate all ctrl function sub-plots
   pl2  = Plots.plot(plotarray..., layout = (params.Ncoupled + params.Nunc, 1))
   pl4  = Plots.plot(plotarray_lab..., layout = (params.Ncoupled + params.Nunc, 1))
   pl5  = Plots.plot(plotarray_fft..., layout = (params.Ncoupled + params.Nunc, 1))
   pl6  = Plots.plot(plotarray_fftlog..., layout = (params.Ncoupled + params.Nunc, 1))

   if savefiles
       rotplotname = ctrlbasename * ".png"
       Plots.savefig(pl2, rotplotname)
       println("Saved rotating frame ctrl plot on file '", rotplotname);

       local labplotname = labbasename * ".png"
       Plots.savefig(pl4, labplotname)
       println("Saved lab frame ctrl plot on file '", labplotname);

       fftname = labbasename * "-fft" * ".png"
       fftname2 = labbasename * "-fft-log" * ".png"
       Plots.savefig(pl5, fftname)
       Plots.savefig(pl6, fftname2)
       println("Saved FFT of lab ctrl function on files '", fftname, "' and '", fftname2, "'");
   end

    # final solution matrix
    # ufinal = unitaryhistory[:,:,end]

    # Return an array of plot objects
    return [pl1, pl2, pl3, pl4, pl5, pl6, plcof, pconv, pluf, plen]

end


text:
module SimIndex

export Index, push!, compile!, k_nearest_neighbors, test_error_ratio

using ProgressMeter
using Base.Collections  # For PriorityQueue.
using Iterators         # For drop.
import Distances
import Base.push!       # We overload push! on SimIndex, need Base.push! as well.

type Index{KeyType, ValueType}
    k::Int  # Number of neighbors
    a::Int  # Size of actual neighborset (> k for convergence)
    distance::Function
    items::Dict{KeyType, ValueType}
    index::Dict{KeyType, PriorityQueue{KeyType, Float64}}
    compiled_index::Dict{KeyType, Array{Pair{KeyType, Float64}, 1}}
    dirty::Bool  # Have elements been pushed to the index since a successful compile?
end

function Index(items; k=10, d=Distances.euclidean)
    if !(typeof(items) <: Dict)
        items = [x => x for x in items]
    end
    KeyType = eltype(keys(items))
    ValueType = eltype(values(items))
    index = Dict{KeyType, PriorityQueue{KeyType, Float64}}()
    compiled_index = Dict{KeyType, Array{Pair{KeyType, Float64}, 1}}()
    si = Index{KeyType, ValueType}(k, 2 * k, d, items, index, compiled_index, true)
    compile!(si)
end

function push!(s::Index, value)
    s.items[value] = value
    s.dirty = true
end

function push!(s::Index, key, value)
    s.items[key] = value
    s.dirty = true
end

function compile!(s::Index, delta::Float64=0.05)
    KeyType = eltype(keys(s.items))
    recompile = (length(s.compiled_index) != 0)
    vs = [k for k in keys(s.items)]
    nvs = length(vs)

    # Initialize the neighbor list. If we've never compiled an index before,
    # we'll sample s.a neighbors at random for each vertex's neighbor list. If
    # this is a recompile of an existing index, we'll keep s.k of the best
    # neighbors for each vertex and sample s.a - s.k neighbors.
    if s.a > nvs - 1
        error("k too large: can't sample $(s.a) neighbors from a set of size $(nvs - 1)")
    end
    p = Progress(nvs, 1, "Sampling vertices for initial neighbor lists...")
    for (v, val) in s.items
        if recompile && haskey(s.compiled_index, v)
            q = [x => d for (x,d) in s.compiled_index[v][1:s.k]]
            s.index[v] = PriorityQueue(q, Base.Sort.Reverse)
            for x in sample(vs, s.a - s.k, union(Set(keys(q)), Set([v])))
                enqueue!(s.index[v], x, s.distance(val, s.items[x]))
            end
        else
            sampled = sample(vs, s.a, Set([v]))
            s.index[v] = PriorityQueue(
              [x => s.distance(val, s.items[x])::Float64 for x in sampled],
              Base.Sort.Reverse)
        end
        next!(p)
    end
    empty!(s.compiled_index)

    res = 1.01  # Decrease res to have more granular progress displayed.
    steps = round(Int, ceil(log(res, 2/delta)))  # Number of steps in progress.
    thres = 2.0  # Next milestone for displaying progress.
    ratio = 2.0  # Current progress (c / nvs)
    p = Progress(steps, 1, "Compiling index...")
    while true
        c = 0
        for j=1:nvs
            u = vs[rand(1:nvs)]
            w = rand_key(s.index[rand_key(s.index[u])])
            if u == w
                continue
            end
            d = s.distance(s.items[u], s.items[w])
            for (x,y) in [(u,w), (w,u)]
                _, maxd = peek(s.index[x])
                if maxd > d
                    c += update_priority_queue(s, x, y, d)
                end
            end
        end
        if c == 0 || ratio < delta
            break
        end
        if c / nvs < ratio
            ratio = c / nvs
        end
        while ratio < thres / res
            next!(p)
            thres = thres / res
        end
    end
    finish!(p)
    generate_compiled_index!(s)
end

function rand_key(d)
    first(take(drop(keys(d), rand(1:length(d)) - 1), 1))
end

# Update SimIndex s to add y to x's nearest neighbors if d is less than one of
# x's existing nearest neighbors' distances. Return 1 if y was added to x's
# nearest neighbors, 0 otherwise.
function update_priority_queue(s, x, y, d)
    try enqueue!(s.index[x], y, d)
    catch
        # y is already in s.index[x]
        return 0
    end
    dequeue!(s.index[x])
    return 1
end

# Sample k elements uniformly from the stream xs while avoiding items in the avoid set.
function sample(xs, k, avoid)
    s = Set{eltype(xs)}()
    nx = length(xs)
    while length(s) < k
        x = xs[rand(1:nx)]
        if !in(x, avoid)
            push!(s, x)
        end
    end
    return [x for x in s]
end

function generate_compiled_index!(s::Index)
    KeyType = eltype(keys(s.index))
    for (key, vals) in s.index
        xs = Pair{KeyType, Float64}[]
        while true
            try
                push!(xs, peek(vals))
                dequeue!(vals)
            catch
                break
            end
        end
        reverse!(xs)
        s.compiled_index[key] = xs
    end
    empty!(s.index)
    s.dirty = false
    s
end

function k_nearest_neighbors(s::Index, key; k=s.k)
    if s.dirty
        error("Index needs to be compiled.")
    end
    vals = get(s.compiled_index, key, Union{})
    if vals == Union{}
        return Union{}
    end
    return vals[1:min(length(vals), k)]
end

function test_error_ratio(s::Index, n=50)
    ks = [k for k in keys(s.items)]
    KeyType = eltype(ks)
    ratios = Float64[]
    for i=1:n
        k = ks[rand(1:end)]
        v = s.items[k]
        q = PriorityQueue(Dict{KeyType, Float64}(), Base.Sort.Reverse)
        for (k2, v2) in s.items
            if k2 == k
                continue
            end
            d = s.distance(v, v2)
            enqueue!(q, k2, d)
            if length(q) > s.k
                dequeue!(q)
            end
        end
        # Accumulate key, distance pairs
        xs = Pair{KeyType, Float64}[]
        while true
            try
                push!(xs, peek(q))
                dequeue!(q)
            catch
                break
            end
        end
        reverse!(xs)
        push!(ratios, error_ratio(xs, k_nearest_neighbors(s, k)))
    end
    mean(ratios)
end

function error_ratio(actual, approx)
    if length(actual) != length(approx)
        error("Arrays passed to error_ratio have different lengths")
    end
    xs = Float64[]
    for i=1:length(actual)
        push!(xs, (approx[i].second + eps()) / (actual[i].second + eps()))
    end
    mean(xs)
end

end # module


text:
# === PRIVATE METHODS BELOW HERE ====================================================================================== #
function design(n::Int64)

    # initialize -
    bitStringArray = Array{String,1}()

    # build the bit strings -
    for i = 0:((2^n) - 1)
        stmp = bitstring(i)
        value = stmp[end - n + 1:end]
        push!(bitStringArray, value)
    end

    # return -
    return bitStringArray
end

# ---------------------------------------------------------------------------------------------------- #
# basic -
function _build_basic_agent_strategy(memory::Int64, score::Int64)::VLBasicMinorityGameStrategy

    try

        # initialize -
        size_of_alphabet = 2
        number_of_elements = (size_of_alphabet)^memory
        outcome_vector = Array{Int,1}(undef, number_of_elements)
        actions = [-1,1]

        # generate an outcome vector -
        for index = 1:number_of_elements

            # generate an index -
            r_index = rand(1:2)

            # outcome -
            outcome_vector[index] = actions[r_index]
        end

        # return -
        return VLBasicMinorityGameStrategy(outcome_vector, score)  # default: score = 0 strategy
    catch error
        rethrow(error)
    end
end

function _build_basic_game_agent(numberOfStrategiesPerAgent::Int64,
    agentMemorySize::Int64)::VLBasicMinorityGameAgent

    try

        # build a collection of strategies for this agent -
        # each agent is going to need numberOfStrategiesPerAgent strategies
        local_strategy_cache = Array{VLBasicMinorityGameStrategy,1}(undef, numberOfStrategiesPerAgent)
        for strategy_index = 1:numberOfStrategiesPerAgent

            # build the strategy - initially all strategies will have a score of zero -
            strategy = _build_basic_agent_strategy(agentMemorySize, 0)

            # package -
            local_strategy_cache[strategy_index] = strategy
        end

        # build agent object - use the first strategy as the best -
        return VLBasicMinorityGameAgent(local_strategy_cache; bestStrategy=first(local_strategy_cache))
    catch error
        # just rethrow the error for now ...
        rethrow(error)
    end
end

function _build_basic_game_world(kwargs_dictionary::Dict)

    try

        # in order to create a game world, we need:
        # numberOfAgents::Int64
        # agentMemorySize::Int64
        # gameAgentArray::Array{VLMinorityGameAgent,1}

        # get data from the args -
        numberOfAgents = kwargs_dictionary[:numberOfAgents]
        agentMemorySize = kwargs_dictionary[:memory]
        numberOfStrategiesPerAgent = kwargs_dictionary[:numberOfStrategiesPerAgent]

        # iniialize -
        gameAgentArray = Array{VLBasicMinorityGameAgent,1}(undef, numberOfAgents)

        # we have two of the three, lets create an array of game agents -
        for agent_index = 1:numberOfAgents

            println("Building agent $(agent_index) of $(numberOfAgents)")

            # package -
            gameAgentArray[agent_index] = _build_basic_game_agent(numberOfStrategiesPerAgent, agentMemorySize)
        end

        # build the game world -
        return VLBasicMinorityGameWorld(numberOfAgents, agentMemorySize, gameAgentArray)
    catch error
        # just rethrow the error for now ...
        rethrow(error)
    end
end
# ---------------------------------------------------------------------------------------------------- #

# ---------------------------------------------------------------------------------------------------- #
# thermal -
function _build_thermal_agent_strategy(memory::Int64, score::Int64, probability::Float64)::VLThermalMinorityGameStrategy

    try

        # initialize -
        actions = [-1,1]
        size_of_alphabet = length(actions)
        number_of_elements = (size_of_alphabet)^memory
        outcome_vector = Array{Int,1}(undef, number_of_elements)

        # generate an outcome vector -
        for index = 1:number_of_elements

            # generate an index -
            r_index = rand(1:2)

            # outcome -
            outcome_vector[index] = actions[r_index]
        end

        # return -
        return VLThermalMinorityGameStrategy(outcome_vector, score, probability)  # default: score = 0 strategy
    catch error
        rethrow(error)
    end
end

function _build_thermal_game_agent(numberOfStrategiesPerAgent::Int64, agentMemorySize::Int64)::VLThermalMinorityGameAgent

    try

        # initialize -
        probability = (1.0 / numberOfStrategiesPerAgent)

        # ok, let's build the collection of strategies for this agent, and the ranking array -
        strategyCollection = Array{VLThermalMinorityGameStrategy,1}(undef, numberOfStrategiesPerAgent)
        for strategy_index = 1:numberOfStrategiesPerAgent

            # build a strategy -
            strategyCollection[strategy_index] = _build_thermal_agent_strategy(agentMemorySize, 0, probability)
        end

        # return -
        return VLThermalMinorityGameAgent(strategyCollection, 0)
    catch error
        rethrow(error)
    end
end

function _build_thermal_game_world(kwargs_dictionary::Dict)::VLThermalMinorityGameWorld

    try

        # initialize -
        numberOfAgents = kwargs_dictionary[:numberOfAgents]
        agentMemorySize = kwargs_dictionary[:memory]
        numberOfStrategiesPerAgent = kwargs_dictionary[:numberOfStrategiesPerAgent]
        temperatute = kwargs_dictionary[:temperature]

        # build an array of game agents -
        agentArray = Array{VLThermalMinorityGameAgent,1}(undef, numberOfAgents)
        for agent_index = 1:numberOfAgents

            # build agents, package -
            agentArray[agent_index] = _build_thermal_game_agent(numberOfStrategiesPerAgent, agentMemorySize)
        end

        # build -
        return VLThermalMinorityGameWorld(agentArray, numberOfAgents, agentMemorySize, temperatute)
    catch error
        rethrow(error)
    end
end
# ---------------------------------------------------------------------------------------------------- #

# ---------------------------------------------------------------------------------------------------- #
# grand cannonical -
function _build_gc_agent_strategy(memory::Int64, score::Int64)::VLGCMinorityGameStrategy

    try

        # initialize -
        size_of_alphabet = 2
        number_of_elements = (size_of_alphabet)^memory
        outcome_vector = Array{Int,1}(undef, number_of_elements)
        actions = [-1,1]

        # generate an outcome vector -
        for index = 1:number_of_elements

            # generate an index -
            r_index = rand(1:2)

            # outcome -
            outcome_vector[index] = actions[r_index]
        end

        # return -
        return VLGCMinorityGameStrategy(outcome_vector, score)  # default: score = 0 strategy
    catch error
        rethrow(error)
    end
end

function _build_gc_game_agent(numberOfStrategiesPerAgent::Int64,
    agentMemorySize::Int64)::VLGCMinorityGameAgent

    try

        # build a collection of strategies for this agent -
        # each agent is going to need numberOfStrategiesPerAgent strategies
        local_strategy_cache = Array{VLGCMinorityGameStrategy,1}(undef, numberOfStrategiesPerAgent)
        for strategy_index = 1:numberOfStrategiesPerAgent

            # build the strategy - initially all strategies will have a score of zero -
            strategy = _build_gc_agent_strategy(agentMemorySize, 0)

            # package -
            local_strategy_cache[strategy_index] = strategy
        end

        # build agent object - use the first strategy as the best -
        return VLGCMinorityGameAgent(local_strategy_cache; bestStrategy=first(local_strategy_cache))
    catch error
        rethrow(error)
    end
end

function _build_gc_game_world(kwargs_dictionary::Dict)::VLGCMinorityGameWorld

    try

        # in order to create a game world, we need:
        # numberOfAgents::Int64
        # agentMemorySize::Int64
        # gameAgentArray::Array{VLMinorityGameAgent,1}

        # get data from the args -
        numberOfAgents = kwargs_dictionary[:numberOfAgents]
        agentMemorySize = kwargs_dictionary[:memory]
        numberOfStrategiesPerAgent = kwargs_dictionary[:numberOfStrategiesPerAgent]

        # iniialize -
        gameAgentArray = Array{VLGCMinorityGameAgent,1}(undef, numberOfAgents)

        # we have two of the three, lets create an array of game agents -
        for agent_index = 1:numberOfAgents

            println("Building agent $(agent_index) of $(numberOfAgents)")

            # package -
            gameAgentArray[agent_index] = _build_gc_game_agent(numberOfStrategiesPerAgent, agentMemorySize)
        end

        # build the game world -
        return VLGCMinorityGameWorld(numberOfAgents, agentMemorySize, gameAgentArray)
    catch error
        rethrow(error)
    end
end
# ---------------------------------------------------------------------------------------------------- #

# === PRIVATE METHODS ABOVE HERE ====================================================================================== #

# === PUBLIC METHODS BELOW HERE ======================================================================================= #

# generic -
function build_minority_game_world(game::Union{Type{VLBasicMinorityGameWorld},Type{VLThermalMinorityGameWorld},Type{VLGCMinorityGameWorld}};
    kwargs...)

    try

        # initialize -
        kwargs_dictionary = Dict(kwargs)
        function_dictionarty = Dict()

        # setup function dictionary -
        function_dictionarty[VLBasicMinorityGameWorld] = _build_basic_game_world
        function_dictionarty[VLThermalMinorityGameWorld] = _build_thermal_game_world
        function_dictionarty[VLGCMinorityGameWorld] = _build_gc_game_world

        # ok, call the linked function -
        return function_dictionarty[game](kwargs_dictionary)
    catch error
        rethrow(error)
    end
end

# === PUBLIC METHODS ABOVE HERE ======================================================================================= #


text:
module IJulia

# in the IPython front-end, enable verbose output via IJulia.set_verbose()
verbose = false
function set_verbose(v=true)
    global verbose::Bool = v
end

using ZMQ
using JSON
using Nettle
include("../deps/ipython.jl")

if isdefined(Base, :REPLCompletions)
    using Base.REPLCompletions
else
    using REPLCompletions
end

uuid4() = repr(Base.Random.uuid4())

inited = false

function init(args)
    if length(args) > 0
        global const profile = open(JSON.parse,args[1])
        verbose && println("PROFILE = $profile")
    else
        # generate profile and save
        let port0 = 5678
            global const profile = (String=>Any)[
                "ip" => "127.0.0.1",
                "transport" => "tcp",
                "stdin_port" => port0,
                "control_port" => port0+1,
                "hb_port" => port0+2,
                "shell_port" => port0+3,
                "iopub_port" => port0+4,
                "key" => uuid4()
            ]
            fname = "profile-$(getpid()).json"
            println("connect ipython with --existing $(pwd())/$fname")
            open(fname, "w") do f
                JSON.print(f, profile)
            end
        end
    end
    global hmac
    if isempty(profile["key"])
        hmac(s1,s2,s3,s4) = ""
    else
        signature_scheme = get(profile, "signature_scheme", "hmac-sha256")
        isempty(signature_scheme) && (signature_scheme = "hmac-sha256")
        signature_scheme = split(signature_scheme, "-")
        if signature_scheme[1] != "hmac" || length(signature_scheme) != 2
            error("unrecognized signature_scheme")
        end
        global const hmacstate = HMACState(eval(symbol(uppercase(signature_scheme[2]))),
                                    profile["key"])
        function hmac(s1,s2,s3,s4)
            update!(hmacstate, s1)
            update!(hmacstate, s2)
            update!(hmacstate, s3)
            update!(hmacstate, s4)
            hexdigest!(hmacstate)
        end
    end

    global const ctx = Context()
    global const publish = Socket(ctx, PUB)
    global const raw_input = Socket(ctx, ROUTER)
    global const requests = Socket(ctx, ROUTER)
    global const control = Socket(ctx, ROUTER)
    global const heartbeat = Socket(ctx, REP)
    bind(publish, "$(profile["transport"])://$(profile["ip"]):$(profile["iopub_port"])")
    bind(requests, "$(profile["transport"])://$(profile["ip"]):$(profile["shell_port"])")
    bind(control, "$(profile["transport"])://$(profile["ip"]):$(profile["control_port"])")
    bind(raw_input, "$(profile["transport"])://$(profile["ip"]):$(profile["stdin_port"])")
    bind(heartbeat, "$(profile["transport"])://$(profile["ip"]):$(profile["hb_port"])")

    global const heartbeat_c = cfunction(heartbeat_thread, Void, (Ptr{Void},))
    global const threadid = Array(Int, 128) # sizeof(uv_thread_t) <= 8 on Linux, OSX, Win
    start_heartbeat(heartbeat)

    global const orig_STDIN = STDIN
    global const orig_STDOUT = STDOUT
    global const orig_STDERR = STDERR

    global read_stdin
    global write_stdin
    global read_stdout
    global write_stdout
    global read_stderr
    global write_stderr
    read_stdin, write_stdin = redirect_stdin()
    read_stdout, write_stdout = redirect_stdout()
    read_stderr, write_stderr = redirect_stderr()

    send_status("starting")
    global inited = true
end

include("stdio.jl")
include("msg.jl")
include("history.jl")
include("handlers.jl")
include("heartbeat.jl")

function eventloop(socket)
    try
        while true
            msg = recv_ipython(socket)
            try
                send_status("busy")
                handlers[msg.header["msg_type"]](socket, msg)
            catch e
                # Try to keep going if we get an exception, but
                # send the exception traceback to the front-ends.
                # (Ignore SIGINT since this may just be a user-requested
                #  kernel interruption to interrupt long calculations.)
                if !isa(e, InterruptException)
                    content = pyerr_content(e, "KERNEL EXCEPTION")
                    map(s -> println(orig_STDERR, s), content["traceback"])
                    send_ipython(publish,
                                 execute_msg == nothing ?
                                 Msg([ "pyerr" ],
                                     [ "msg_id" => uuid4(),
                                     "username" => "jlkernel",
                                     "session" => uuid4(),
                                      "msg_type" => "pyerr" ], content) :
                                 msg_pub(execute_msg, "pyerr", content))
                end
            finally
                send_status("idle")
            end
        end
    catch e
        # the IPython manager may send us a SIGINT if the user
        # chooses to interrupt the kernel; don't crash on this
        if isa(e, InterruptException)
            eventloop(socket)
        else
            rethrow()
        end
    end
end

function waitloop()
    @async eventloop(control)
    requests_task = @async eventloop(requests)
    while true
        try
            wait()
        catch e
            # send interrupts (user SIGINT) to the code-execution task
            if isa(e, InterruptException)
                @async Base.throwto(requests_task, e)
            else
                rethrow()
            end
        end
    end
end

export notebook
function notebook(ipython=find_ipython()[1])
    inited && error("IJulia is already running")
    run(`$ipython notebook --profile julia`)
end

end # IJulia