text:
function Base.map(f, A::KeyedArray)
    data = map(f, parent(A))
    KeyedArray(data, map(copy, axiskeys(A)))#, copy(A.meta))
end
for fun in [:map, :map!], (T, S) in [ (:KeyedArray, :KeyedArray),
        (:KeyedArray, :AbstractArray), (:AbstractArray, :KeyedArray),
        (:KeyedArray, :NamedDimsArray), (:NamedDimsArray, :KeyedArray)] # for ambiguities

    @eval function Base.$fun(f, A::$T, B::$S, Cs::AbstractArray...)
        data = $fun(f, keyless(A), keyless(B), keyless.(Cs)...)
        new_keys = unify_keys(keys_or_axes(A), keys_or_axes(B), keys_or_axes.(Cs)...)
        KeyedArray(data, map(copy, new_keys)) # copy sometimes wasteful for map!, but OK.
    end
end

using Base: Generator

function Base.collect(x::Generator{<:KeyedArray})
    data = collect(Generator(x.f, x.iter.data))
    KeyedArray(data, map(copy, axiskeys(x.iter)))#, copy(A.meta))
end
function Base.collect(x::Generator{<:Iterators.Enumerate{<:KeyedArray}})
    data = collect(Generator(x.f, enumerate(x.iter.itr.data)))
    KeyedArray(data, map(copy, axiskeys(x.iter.itr)))
end
for Ts in [(:KeyedArray,), (:KeyedArray, :NamedDimsArray), (:NamedDimsArray, :KeyedArray)]
    @eval function Base.collect(x::Generator{<:Iterators.ProductIterator{<:Tuple{$(Ts...),Vararg{Any}}}})
        data = collect(Generator(x.f, Iterators.product(keyless.(x.iter.iterators)...)))
        all_keys = tuple_flatten(keys_or_axes.(x.iter.iterators)...)
        KeyedArray(data, map(copy, all_keys))
    end
end

tuple_flatten(x::Tuple, ys::Tuple...) = (x..., tuple_flatten(ys...)...)
tuple_flatten() = ()

function Base.mapreduce(f, op, A::KeyedArray; dims=:, kwargs...) # sum, prod, etc
    dims === Colon() && return mapreduce(f, op, parent(A))
    numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    data = mapreduce(f, op, parent(A); dims=numerical_dims, kwargs...)
    new_keys = ntuple(d -> d in numerical_dims ? Base.OneTo(1) : axiskeys(A,d), ndims(A))
    return KeyedArray(data, map(copy, new_keys))#, copy(A.meta))
end

using Statistics
for fun in [:mean, :std, :var] # These don't use mapreduce, but could perhaps be handled better?
    @eval function Statistics.$fun(A::KeyedArray; dims=:)
        dims === Colon() && return $fun(parent(A))
        numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
        data = $fun(parent(A); dims=numerical_dims)
        new_keys = ntuple(d -> d in numerical_dims ? Base.OneTo(1) : axiskeys(A,d), ndims(A))
        return KeyedArray(data, map(copy, new_keys))#, copy(A.meta))
    end
    VERSION >= v"1.3" &&
    @eval function Statistics.$fun(f, A::KeyedArray; dims=:)
        dims === Colon() && return $fun(f, parent(A))
        numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
        data = $fun(f, parent(A); dims=numerical_dims)
        new_keys = ntuple(d -> d in numerical_dims ? Base.OneTo(1) : axiskeys(A,d), ndims(A))
        return KeyedArray(data, map(copy, new_keys))#, copy(A.meta))
    end
end

function Base.dropdims(A::KeyedArray; dims)
    numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    data = dropdims(parent(A); dims=dims)
    new_keys = key_skip(axiskeys(A), numerical_dims...)
    KeyedArray(data, new_keys)#, A.meta)
end

key_skip(tup::Tuple, d, dims...) = key_skip(
    ntuple(n -> n<d ? tup[n] : tup[n+1], length(tup)-1),
    map(n -> n<d ? n : n-1, dims)...)
key_skip(tup::Tuple) = tup

function Base.permutedims(A::KeyedArray, perm)
    numerical_perm = hasnames(A) ? NamedDims.dim(dimnames(A), perm) : perm
    data = permutedims(parent(A), numerical_perm)
    new_keys = ntuple(d -> copy(axiskeys(A, perm[d])), ndims(A))
    KeyedArray(data, new_keys)#, copy(A.meta))
end

if VERSION >= v"1.1"
    # This copies the implementation from Base, except with numerical_dims:
    @inline function Base.eachslice(A::KeyedArray; dims)
        numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
        length(numerical_dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
        dim = first(numerical_dims)
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        inds_before = ntuple(d->(:), dim-1)
        inds_after = ntuple(d->(:), ndims(A)-dim)
        return (view(A, inds_before..., i, inds_after...) for i in axes(A, dim))
    end
end

function Base.mapslices(f, A::KeyedArray; dims)
    numerical_dims = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    data = mapslices(f, parent(A); dims=dims)
    new_keys = ntuple(ndims(A)) do d
        d in dims ? axes(data,d) : copy(axiskeys(A, d))
    end
    KeyedArray(data, new_keys)#, copy(A.meta))
end

for (T, S) in [(:KeyedVecOrMat, :KeyedVecOrMat), # KeyedArray gives ambiguities
    (:KeyedVecOrMat, :AbstractVecOrMat), (:AbstractVecOrMat, :KeyedVecOrMat),
    (:NdaKaVoM, :NdaKaVoM),
    (:NdaKaVoM, :KeyedVecOrMat), (:KeyedVecOrMat, :NdaKaVoM),
    (:NdaKaVoM, :AbstractVecOrMat), (:AbstractVecOrMat, :NdaKaVoM),
    ]

    @eval function Base.vcat(A::$T, B::$S, Cs::AbstractVecOrMat...)
        data = vcat(keyless(A), keyless(B), keyless.(Cs)...)
        new_1 = key_vcat(keys_or_axes(A,1), keys_or_axes(B,1), keys_or_axes.(Cs,1)...)
        new_keys = ndims(A) == 1 ? (new_1,) :
            (new_1, unify_one(keys_or_axes(A,2), keys_or_axes(B,2), keys_or_axes.(Cs,2)...))
        KeyedArray(data, map(copy, new_keys))
    end

    @eval function Base.hcat(A::$T, B::$S, Cs::AbstractVecOrMat...)
        data = hcat(keyless(A), keyless(B), keyless.(Cs)...)
        new_1 = unify_one(keys_or_axes(A,1), keys_or_axes(B,1), keys_or_axes.(Cs,1)...)
        new_2 = ndims(A) == 1 ? axes(data,2) :
            key_vcat(keys_or_axes(A,2), keys_or_axes(B,2), keys_or_axes.(Cs,2)...)
        KeyedArray(data, map(copy, (new_1, new_2)))
    end

end
for (T, S) in [ (:KeyedArray, :KeyedArray),
        (:KeyedArray, :AbstractArray), (:AbstractArray, :KeyedArray),
        (:KeyedArray, :NamedDimsArray), (:NamedDimsArray, :KeyedArray),
        (:NdaKa, :NdaKa),
        (:NdaKa, :KeyedArray), (:KeyedArray, :NdaKa),
        (:NdaKa, :AbstractArray), (:AbstractArray, :NdaKa),
        ]

    @eval function Base.cat(A::$T, B::$S, Cs::AbstractArray...; dims)
        numerical_dims, data = if any(hasnames.((A, B, Cs...)))
            old_names = NamedDims.unify_names_longest(dimnames(A), dimnames(B), dimnames.(Cs)...)
            new_names = NamedDims.expand_dimnames(old_names, dims)
            α = NamedDims.dim(new_names, dims)
            β = cat(keyless(A), keyless(B), keyless.(Cs)...; dims=dims)
            α, β
        else
            α = val_strip(dims)
            β = cat(keyless(A), keyless(B), keyless.(Cs)...; dims=numerical_dims)
            α, β
        end
        new_keys = ntuple(ndims(data)) do d
            if d in numerical_dims
                key_vcat(keys_or_axes(A,d), keys_or_axes(B,d), keys_or_axes.(Cs,d)...)
            else
                unify_one(keys_or_axes(A,d), keys_or_axes(B,d), keys_or_axes.(Cs,d)...)
            end
        end
        KeyedArray(data, map(copy, new_keys)) # , copy(A.meta))
    end

end
val_strip(dims::Val{d}) where {d} = d
val_strip(dims) = dims
key_vcat(a::AbstractVector, b::AbstractVector) = vcat(a,b)
key_vcat(a::Base.OneTo, b::Base.OneTo) = Base.OneTo(a.stop + b.stop)
key_vcat(a,b,cs...) = key_vcat(key_vcat(a,b),cs...)

function Base.sort(A::KeyedArray; dims, kw...)
    dims′ = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    data = sort(parent(A); dims=dims′, kw...)
    # sorts each (say) col independently, thus keys along them loses meaning.
    new_keys = ntuple(d -> d==dims′ ? OneTo(size(A,d)) : axiskeys(A,d), ndims(A))
    KeyedArray(data, map(copy, new_keys)) # , copy(A.meta))
end
function Base.sort(A::KeyedVector; kw...)
    perm = sortperm(parent(A); kw...)
    KeyedArray(parent(A)[perm], (axiskeys(A,1)[perm],)) # , copy(A.meta))
end

function Base.sort!(A::KeyedVector; kw...)
    perm = sortperm(parent(A); kw...)
    permute!(axiskeys(A,1), perm) # error if keys cannot be sorted, could treat like push!
    permute!(parent(A), perm)
    A
end

sort_doc = """
    sortslices(A; dims)
    sortkeys(A; dims=1:ndims(A))

`Base.sortslices` sorts the corresponding keys too, along one dimension.
Calls its own implementation, roughly `p = sortperm(eachslice(A))`,
with default keyword `by=vec` to make this work on slices of any shape.

`sortkeys(A)` instead sorts everything by the keys.
Works along any number of dimensions, by detault all of them.
"""

@doc sort_doc
function Base.sortslices(A::KeyedArray; dims, by=vec, kw...)
    dim′ = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    perms = ntuple(ndims(A)) do d
        d in dim′ || return Colon()
        sortperm(collect(eachslice(parent(A), dims=dim′)); by=by, kw...)
    end
    new_keys = map(getindex, axiskeys(A), perms)
    KeyedArray(keyless(A)[perms...], new_keys) # , copy(A.meta))
end

if VERSION < v"1.1" # defn copied Julia 1.4 Base abstractarraymath.jl:452
    @inline function eachslice(A::AbstractArray; dims)
        length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
        dim = first(dims)
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        inds_before = ntuple(d->(:), dim-1)
        inds_after = ntuple(d->(:), ndims(A)-dim)
        return (view(A, inds_before..., i, inds_after...) for i in axes(A, dim))
    end
end

@doc sort_doc
function sortkeys(A::Union{KeyedArray, NdaKa}; dims=1:ndims(A), kw...)
    dims′ = hasnames(A) ? NamedDims.dim(dimnames(A), dims) : dims
    perms = ntuple(ndims(A)) do d
        d in dims′ || return Colon()
        axiskeys(A,d) isa AbstractUnitRange && return Colon() # avoids OneTo(n) -> 1:n
        sortperm(axiskeys(A,d); kw...)
    end
    new_keys = map(getindex, axiskeys(A), perms)
    KeyedArray(keyless(A)[perms...], new_keys) # , copy(A.meta))
end

Base.filter(f, A::KeyedVector) = getindex(A, map(f, parent(A)))
Base.filter(f, A::KeyedArray) = filter(f, parent(A))

using LinearAlgebra

for (mod, fun, lazy) in [(Base, :permutedims, false),
        (LinearAlgebra, :transpose, true), (LinearAlgebra, :adjoint, true)]
    @eval function $mod.$fun(A::KeyedArray)
        data = $mod.$fun(parent(A))
        new_keys = ndims(A)==1 ? (Base.OneTo(1), axiskeys(A,1)) :
            ndims(data)==1 ? (axiskeys(A,2),) :
            reverse(axiskeys(A))
        KeyedArray(data, $(lazy ? :new_keys : :(map(copy, new_keys))))#, $(lazy ? :(A.meta) : :(copy(A.meta))))
    end
end

Base.reshape(A::KeyedArray, dims::Tuple{Vararg{Int}}) = reshape(parent(A), dims...) # for ambiguities
Base.reshape(A::KeyedArray, dims::Tuple{Vararg{Union{Colon, Int}}}) = reshape(parent(A), dims...)

for fun in [:copy, :deepcopy, :similar, :zero, :one]
    @eval Base.$fun(A::KeyedArray) = KeyedArray($fun(parent(A)), map(copy, axiskeys(A)))
    @eval Base.$fun(A::NdaKa) = NamedDimsArray(KeyedArray(
        $fun(parent(parent(A))),
        map(copy, axiskeys(A))), dimnames(A))
end
Base.similar(A::KeyedArray, T::Type) = KeyedArray(similar(parent(A), T), map(copy, axiskeys(A)))
Base.similar(A::NdaKa, T::Type) = NamedDimsArray(KeyedArray(
    similar(parent(parent(A)), T), map(copy, axiskeys(A))), dimnames(A))
Base.similar(A::KeyedArray, T::Type, dims::Int...) = similar(parent(A), T, dims...)
Base.similar(A::KeyedArray, dims::Int...) = similar(parent(A), dims...)

for fun in [:(==), :isequal, :isapprox]
    for (T, S) in [ (:KeyedArray, :KeyedArray), (:KeyedArray, :NdaKa), (:NdaKa, :KeyedArray) ]
        @eval function Base.$fun(A::$T, B::$S; kw...)
            # Ideally you would pass isapprox(, atol) into unifiable_keys?
            unifiable_keys(axiskeys(A), axiskeys(B)) || return false
            return $fun(keyless(A), keyless(B); kw...)
        end
    end
end

Rlist = [:KeyedMatrix, :KeyedVector,
    :(NdaKa{L,T,2} where {L,T}), :(NdaKa{L,T,1} where {L,T}),
    :(KeyedVector{T} where {T<:Number}), :(NdaKa{L,T,1} where {L,T<:Number}), # ambiguities on 1.5
    ]

Olist = [ :AbstractMatrix, :AbstractVector, :Number,
    :(Adjoint{<:Any,<:AbstractMatrix}), :(Adjoint{<:Any,<:AbstractVector}),
    :(Transpose{<:Any,<:AbstractMatrix}), :(Transpose{<:Any,<:AbstractVector}),
    :(NamedDimsArray{L,T,1} where {L,T}), :(NamedDimsArray{L,T,2} where {L,T}),
    :(Adjoint{<:Number,<:AbstractVector}), # 1.5 problem...
    ]
for (Ts, Ss) in [(Rlist, Rlist), (Rlist, Olist), (Olist, Rlist)]
    for T in Ts, S in Ss # some combinations are errors, later, that's ok

        @eval Base.:*(x::$T, y::$S) = matmul(x,y)
        @eval Base.:\(x::$T, y::$S) = ldiv(x,y)
        @eval Base.:/(x::$T, y::$S) = rdiv(x,y)

    end
end
for (fun, op) in [(:matmul, :*), (:ldiv, :\), (:rdiv, :/)]
    @eval $fun(x::AbstractVecOrMat, y::Number) = KeyedArray($op(keyless(x), y), axiskeys(x))
    @eval $fun(x::Number, y::AbstractVecOrMat) = KeyedArray($op(x, keyless(y)), axiskeys(y))
    @eval $fun(x::AbstractVector, y::AbstractVector) = $op(keyless(x), keyless(y))
end

function matmul(x::AbstractMatrix, y::AbstractVecOrMat)
    data = keyless(x) * keyless(y)
    unify_one(keys_or_axes(x,2), keys_or_axes(y,1)) # just a check, discard these
    if data isa AbstractVecOrMat
        new_keys = (keys_or_axes(x,1), Base.tail(keys_or_axes(y))...)
        KeyedArray(data, map(copy, new_keys))
    else
        data # case V' * V
    end
end
function matmul(x::AbstractVector, y::AbstractMatrix) # used for v * v'
    data = keyless(x) * keyless(y)
    new_keys = (keys_or_axes(x,1), keys_or_axes(y,2))
    KeyedArray(data, map(copy, new_keys))
end

# case of two vectors gives a scalar, caught above.
function ldiv(x::AbstractVecOrMat, y::AbstractVecOrMat)
    data = keyless(x) \ keyless(y)
    unify_one(keys_or_axes(x,1), keys_or_axes(y,1))
    new_keys = (Base.tail(keys_or_axes(x))..., Base.tail(keys_or_axes(y))...)
    KeyedArray(data, map(copy, new_keys))
end
function rdiv(x::AbstractVecOrMat, y::AbstractVecOrMat)
    data = keyless(x) / keyless(y)
    # unify_one(keys_or_axes(x,2), keys_or_axes(y,2)) # not right!
    # new_keys = (tup_head(keys_or_axes(x))..., tup_head(keys_or_axes(y))...)
    # KeyedArray(data, new_keys)
    @warn "/ doesn't preserve keys yet, sorry" maxlog=1
    data
end

tup_head(t::Tuple) = reverse(Base.tail(reverse(t)))

for fun in [:inv, :pinv,
    :det, :logdet, :logabsdet,
    :eigen, :eigvecs, :eigvals, :svd,
    :diag
    ]
    @eval LinearAlgebra.$fun(A::KeyedMatrix) = $fun(parent(A))
end

function LinearAlgebra.cholesky(A::Hermitian{T, <:KeyedArray{T}}; kwargs...) where T
    return cholesky(parent(A); kwargs...)
end

function LinearAlgebra.cholesky(A::KeyedMatrix; kwargs...)
    data = hasnames(A) ? parent(parent(A)) : parent(A)
    return cholesky(data; kwargs...)
end


text:
#
# Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
# Copyright 2015 Oliver Heimlich (oheim@posteo.de)
#
# Original author: Marco Nehmeier (unit tests in libieeep1788)
# Converted into portable ITL format by Oliver Heimlich with minor corrections.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#Language imports

#Test library imports
using FactCheck

#Arithmetic library imports
using ValidatedNumerics

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_empty_test") do
    @fact isempty(∅) --> true
    @fact isempty(Interval(-Inf, Inf)) --> false
    @fact isempty(Interval(1.0, 2.0)) --> false
    @fact isempty(Interval(-1.0, 2.0)) --> false
    @fact isempty(Interval(-3.0, -2.0)) --> false
    @fact isempty(Interval(-Inf, 2.0)) --> false
    @fact isempty(Interval(-Inf, 0.0)) --> false
    @fact isempty(Interval(-Inf, -0.0)) --> false
    @fact isempty(Interval(0.0, Inf)) --> false
    @fact isempty(Interval(-0.0, Inf)) --> false
    @fact isempty(Interval(-0.0, 0.0)) --> false
    @fact isempty(Interval(0.0, -0.0)) --> false
    @fact isempty(Interval(0.0, 0.0)) --> false
    @fact isempty(Interval(-0.0, -0.0)) --> false
end

facts("minimal_empty_dec_test") do
    @fact isempty(DecoratedInterval(∅, trv)) --> true
    @fact isempty(DecoratedInterval(Interval(-Inf, Inf), def)) --> false
    @fact isempty(DecoratedInterval(Interval(1.0, 2.0), com)) --> false
    @fact isempty(DecoratedInterval(Interval(-1.0, 2.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(-3.0, -2.0), dac)) --> false
    @fact isempty(DecoratedInterval(Interval(-Inf, 2.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(-Inf, 0.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(-Inf, -0.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(0.0, Inf), def)) --> false
    @fact isempty(DecoratedInterval(Interval(-0.0, Inf), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(-0.0, 0.0), com)) --> false
    @fact isempty(DecoratedInterval(Interval(0.0, -0.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(0.0, 0.0), trv)) --> false
    @fact isempty(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
end

facts("minimal_entire_test") do
    @fact isentire(∅) --> false
    @fact isentire(Interval(-Inf, Inf)) --> true
    @fact isentire(Interval(1.0, 2.0)) --> false
    @fact isentire(Interval(-1.0, 2.0)) --> false
    @fact isentire(Interval(-3.0, -2.0)) --> false
    @fact isentire(Interval(-Inf, 2.0)) --> false
    @fact isentire(Interval(-Inf, 0.0)) --> false
    @fact isentire(Interval(-Inf, -0.0)) --> false
    @fact isentire(Interval(0.0, Inf)) --> false
    @fact isentire(Interval(-0.0, Inf)) --> false
    @fact isentire(Interval(-0.0, 0.0)) --> false
    @fact isentire(Interval(0.0, -0.0)) --> false
    @fact isentire(Interval(0.0, 0.0)) --> false
    @fact isentire(Interval(-0.0, -0.0)) --> false
end

facts("minimal_entire_dec_test") do
    @fact isentire(DecoratedInterval(∅, trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-Inf, Inf), trv)) --> true
    @fact isentire(DecoratedInterval(Interval(-Inf, Inf), def)) --> true
    @fact isentire(DecoratedInterval(Interval(-Inf, Inf), dac)) --> true
    @fact isentire(DecoratedInterval(Interval(1.0, 2.0), com)) --> false
    @fact isentire(DecoratedInterval(Interval(-1.0, 2.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-3.0, -2.0), dac)) --> false
    @fact isentire(DecoratedInterval(Interval(-Inf, 2.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-Inf, 0.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-Inf, -0.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(0.0, Inf), def)) --> false
    @fact isentire(DecoratedInterval(Interval(-0.0, Inf), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-0.0, 0.0), com)) --> false
    @fact isentire(DecoratedInterval(Interval(0.0, -0.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(0.0, 0.0), trv)) --> false
    @fact isentire(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
end

facts("minimal_nai_dec_test") do
    @fact isnai(DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-Inf, Inf), def)) --> false
    @fact isnai(DecoratedInterval(Interval(-Inf, Inf), dac)) --> false
    @fact isnai(DecoratedInterval(Interval(1.0, 2.0), com)) --> false
    @fact isnai(DecoratedInterval(Interval(-1.0, 2.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-3.0, -2.0), dac)) --> false
    @fact isnai(DecoratedInterval(Interval(-Inf, 2.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-Inf, 0.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-Inf, -0.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(0.0, Inf), def)) --> false
    @fact isnai(DecoratedInterval(Interval(-0.0, Inf), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-0.0, 0.0), com)) --> false
    @fact isnai(DecoratedInterval(Interval(0.0, -0.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(0.0, 0.0), trv)) --> false
    @fact isnai(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
end

facts("minimal_equal_test") do
    @fact Interval(1.0, 2.0) == Interval(1.0, 2.0) --> true
    @fact Interval(1.0, 2.1) == Interval(1.0, 2.0) --> false
    @fact ∅ == ∅ --> true
    @fact ∅ == Interval(1.0, 2.0) --> false
    @fact Interval(-Inf, Inf) == Interval(-Inf, Inf) --> true
    @fact Interval(1.0, 2.4) == Interval(-Inf, Inf) --> false
    @fact Interval(1.0, Inf) == Interval(1.0, Inf) --> true
    @fact Interval(1.0, 2.4) == Interval(1.0, Inf) --> false
    @fact Interval(-Inf, 2.0) == Interval(-Inf, 2.0) --> true
    @fact Interval(-Inf, 2.4) == Interval(-Inf, 2.0) --> false
    @fact Interval(-2.0, 0.0) == Interval(-2.0, 0.0) --> true
    @fact Interval(-0.0, 2.0) == Interval(0.0, 2.0) --> true
    @fact Interval(-0.0, -0.0) == Interval(0.0, 0.0) --> true
    @fact Interval(-0.0, 0.0) == Interval(0.0, 0.0) --> true
    @fact Interval(0.0, -0.0) == Interval(0.0, 0.0) --> true
end

facts("minimal_equal_dec_test") do
    @fact DecoratedInterval(Interval(1.0, 2.0), def) == DecoratedInterval(Interval(1.0, 2.0), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.1), trv) == DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact DecoratedInterval(∅, trv) == DecoratedInterval(∅, trv) --> true
    @fact DecoratedInterval(∅, trv) == DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact DecoratedInterval(∅, trv) == DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), def) == DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.4), trv) == DecoratedInterval(Interval(-Inf, Inf), trv) --> false
    @fact DecoratedInterval(Interval(1.0, Inf), trv) == DecoratedInterval(Interval(1.0, Inf), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.4), def) == DecoratedInterval(Interval(1.0, Inf), trv) --> false
    @fact DecoratedInterval(Interval(-Inf, 2.0), trv) == DecoratedInterval(Interval(-Inf, 2.0), trv) --> true
    @fact DecoratedInterval(Interval(-Inf, 2.4), def) == DecoratedInterval(Interval(-Inf, 2.0), trv) --> false
    @fact DecoratedInterval(Interval(-2.0, 0.0), trv) == DecoratedInterval(Interval(-2.0, 0.0), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, 2.0), def) == DecoratedInterval(Interval(0.0, 2.0), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, -0.0), trv) == DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, 0.0), def) == DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact DecoratedInterval(Interval(0.0, -0.0), trv) == DecoratedInterval(Interval(0.0, 0.0), trv) --> true
end

facts("minimal_subset_test") do
    @fact ∅ ⊆ ∅ --> true
    @fact ∅ ⊆ Interval(0.0, 4.0) --> true
    @fact ∅ ⊆ Interval(-0.0, 4.0) --> true
    @fact ∅ ⊆ Interval(-0.1, 1.0) --> true
    @fact ∅ ⊆ Interval(-0.1, 0.0) --> true
    @fact ∅ ⊆ Interval(-0.1, -0.0) --> true
    @fact ∅ ⊆ Interval(-Inf, Inf) --> true
    @fact Interval(0.0, 4.0) ⊆ ∅ --> false
    @fact Interval(-0.0, 4.0) ⊆ ∅ --> false
    @fact Interval(-0.1, 1.0) ⊆ ∅ --> false
    @fact Interval(-Inf, Inf) ⊆ ∅ --> false
    @fact Interval(0.0, 4.0) ⊆ Interval(-Inf, Inf) --> true
    @fact Interval(-0.0, 4.0) ⊆ Interval(-Inf, Inf) --> true
    @fact Interval(-0.1, 1.0) ⊆ Interval(-Inf, Inf) --> true
    @fact Interval(-Inf, Inf) ⊆ Interval(-Inf, Inf) --> true
    @fact Interval(1.0, 2.0) ⊆ Interval(1.0, 2.0) --> true
    @fact Interval(1.0, 2.0) ⊆ Interval(0.0, 4.0) --> true
    @fact Interval(1.0, 2.0) ⊆ Interval(-0.0, 4.0) --> true
    @fact Interval(0.1, 0.2) ⊆ Interval(0.0, 4.0) --> true
    @fact Interval(0.1, 0.2) ⊆ Interval(-0.0, 4.0) --> true
    @fact Interval(-0.1, -0.1) ⊆ Interval(-4.0, 3.4) --> true
    @fact Interval(0.0, 0.0) ⊆ Interval(-0.0, -0.0) --> true
    @fact Interval(-0.0, -0.0) ⊆ Interval(0.0, 0.0) --> true
    @fact Interval(-0.0, 0.0) ⊆ Interval(0.0, 0.0) --> true
    @fact Interval(-0.0, 0.0) ⊆ Interval(0.0, -0.0) --> true
    @fact Interval(0.0, -0.0) ⊆ Interval(0.0, 0.0) --> true
    @fact Interval(0.0, -0.0) ⊆ Interval(-0.0, 0.0) --> true
end

facts("minimal_subset_dec_test") do
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(0.0, 4.0), trv) --> true
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(-0.0, 4.0), def) --> true
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(-0.1, 1.0), trv) --> true
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(-0.1, 0.0), trv) --> true
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(-0.1, -0.0), trv) --> true
    @fact DecoratedInterval(∅, trv) ⊆ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) ⊆ DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(Interval(-0.0, 4.0), def) ⊆ DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(Interval(-0.1, 1.0), trv) ⊆ DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ⊆ DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) ⊆ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, 4.0), trv) ⊆ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(-0.1, 1.0), trv) ⊆ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ⊆ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ⊆ DecoratedInterval(Interval(1.0, 2.0), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ⊆ DecoratedInterval(Interval(0.0, 4.0), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), def) ⊆ DecoratedInterval(Interval(-0.0, 4.0), def) --> true
    @fact DecoratedInterval(Interval(0.1, 0.2), trv) ⊆ DecoratedInterval(Interval(0.0, 4.0), trv) --> true
    @fact DecoratedInterval(Interval(0.1, 0.2), trv) ⊆ DecoratedInterval(Interval(-0.0, 4.0), def) --> true
    @fact DecoratedInterval(Interval(-0.1, -0.1), trv) ⊆ DecoratedInterval(Interval(-4.0, 3.4), trv) --> true
    @fact DecoratedInterval(Interval(0.0, 0.0), trv) ⊆ DecoratedInterval(Interval(-0.0, -0.0), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, -0.0), trv) ⊆ DecoratedInterval(Interval(0.0, 0.0), def) --> true
    @fact DecoratedInterval(Interval(-0.0, 0.0), trv) ⊆ DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact DecoratedInterval(Interval(-0.0, 0.0), trv) ⊆ DecoratedInterval(Interval(0.0, -0.0), trv) --> true
    @fact DecoratedInterval(Interval(0.0, -0.0), def) ⊆ DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact DecoratedInterval(Interval(0.0, -0.0), trv) ⊆ DecoratedInterval(Interval(-0.0, 0.0), trv) --> true
end

facts("minimal_less_test") do
    @fact <=(∅, ∅) --> true
    @fact ∅ ≤ ∅ --> true
    @fact <=(Interval(1.0, 2.0), ∅) --> false
    @fact Interval(1.0, 2.0) ≤ ∅ --> false
    @fact <=(∅, Interval(1.0, 2.0)) --> false
    @fact ∅ ≤ Interval(1.0, 2.0) --> false
    @fact <=(Interval(-Inf, Inf), Interval(-Inf, Inf)) --> true
    @fact Interval(-Inf, Inf) ≤ Interval(-Inf, Inf) --> true
    @fact <=(Interval(1.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact Interval(1.0, 2.0) ≤ Interval(-Inf, Inf) --> false
    @fact <=(Interval(0.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact Interval(0.0, 2.0) ≤ Interval(-Inf, Inf) --> false
    @fact <=(Interval(-0.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact Interval(-0.0, 2.0) ≤ Interval(-Inf, Inf) --> false
    @fact <=(Interval(-Inf, Inf), Interval(1.0, 2.0)) --> false
    @fact Interval(-Inf, Inf) ≤ Interval(1.0, 2.0) --> false
    @fact <=(Interval(-Inf, Inf), Interval(0.0, 2.0)) --> false
    @fact Interval(-Inf, Inf) ≤ Interval(0.0, 2.0) --> false
    @fact <=(Interval(-Inf, Inf), Interval(-0.0, 2.0)) --> false
    @fact Interval(-Inf, Inf) ≤ Interval(-0.0, 2.0) --> false
    @fact <=(Interval(0.0, 2.0), Interval(0.0, 2.0)) --> true
    @fact Interval(0.0, 2.0) ≤ Interval(0.0, 2.0) --> true
    @fact <=(Interval(0.0, 2.0), Interval(-0.0, 2.0)) --> true
    @fact Interval(0.0, 2.0) ≤ Interval(-0.0, 2.0) --> true
    @fact <=(Interval(0.0, 2.0), Interval(1.0, 2.0)) --> true
    @fact Interval(0.0, 2.0) ≤ Interval(1.0, 2.0) --> true
    @fact <=(Interval(-0.0, 2.0), Interval(1.0, 2.0)) --> true
    @fact Interval(-0.0, 2.0) ≤ Interval(1.0, 2.0) --> true
    @fact <=(Interval(1.0, 2.0), Interval(1.0, 2.0)) --> true
    @fact Interval(1.0, 2.0) ≤ Interval(1.0, 2.0) --> true
    @fact <=(Interval(1.0, 2.0), Interval(3.0, 4.0)) --> true
    @fact Interval(1.0, 2.0) ≤ Interval(3.0, 4.0) --> true
    @fact <=(Interval(1.0, 3.5), Interval(3.0, 4.0)) --> true
    @fact Interval(1.0, 3.5) ≤ Interval(3.0, 4.0) --> true
    @fact <=(Interval(1.0, 4.0), Interval(3.0, 4.0)) --> true
    @fact Interval(1.0, 4.0) ≤ Interval(3.0, 4.0) --> true
    @fact <=(Interval(-2.0, -1.0), Interval(-2.0, -1.0)) --> true
    @fact Interval(-2.0, -1.0) ≤ Interval(-2.0, -1.0) --> true
    @fact <=(Interval(-3.0, -1.5), Interval(-2.0, -1.0)) --> true
    @fact Interval(-3.0, -1.5) ≤ Interval(-2.0, -1.0) --> true
    @fact <=(Interval(0.0, 0.0), Interval(-0.0, -0.0)) --> true
    @fact Interval(0.0, 0.0) ≤ Interval(-0.0, -0.0) --> true
    @fact <=(Interval(-0.0, -0.0), Interval(0.0, 0.0)) --> true
    @fact Interval(-0.0, -0.0) ≤ Interval(0.0, 0.0) --> true
    @fact <=(Interval(-0.0, 0.0), Interval(0.0, 0.0)) --> true
    @fact Interval(-0.0, 0.0) ≤ Interval(0.0, 0.0) --> true
    @fact <=(Interval(-0.0, 0.0), Interval(0.0, -0.0)) --> true
    @fact Interval(-0.0, 0.0) ≤ Interval(0.0, -0.0) --> true
    @fact <=(Interval(0.0, -0.0), Interval(0.0, 0.0)) --> true
    @fact Interval(0.0, -0.0) ≤ Interval(0.0, 0.0) --> true
    @fact <=(Interval(0.0, -0.0), Interval(-0.0, 0.0)) --> true
    @fact Interval(0.0, -0.0) ≤ Interval(-0.0, 0.0) --> true
end

facts("minimal_less_dec_test") do
    @fact <=(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(∅, trv)) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ≤ DecoratedInterval(∅, trv) --> false
    @fact <=(DecoratedInterval(∅, trv), DecoratedInterval(Interval(1.0, 2.0), def)) --> false
    @fact DecoratedInterval(∅, trv) ≤ DecoratedInterval(Interval(1.0, 2.0), def) --> false
    @fact <=(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(∅, trv)) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ≤ DecoratedInterval(∅, trv) --> false
    @fact <=(DecoratedInterval(∅, trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> false
    @fact DecoratedInterval(∅, trv) ≤ DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact <=(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> true
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ≤ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact <=(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), def) ≤ DecoratedInterval(Interval(-Inf, Inf), trv) --> false
    @fact <=(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact DecoratedInterval(Interval(0.0, 2.0), trv) ≤ DecoratedInterval(Interval(-Inf, Inf), trv) --> false
    @fact <=(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact DecoratedInterval(Interval(-0.0, 2.0), trv) ≤ DecoratedInterval(Interval(-Inf, Inf), trv) --> false
    @fact <=(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ≤ DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact <=(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(0.0, 2.0), def)) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ≤ DecoratedInterval(Interval(0.0, 2.0), def) --> false
    @fact <=(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-0.0, 2.0), trv)) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ≤ DecoratedInterval(Interval(-0.0, 2.0), trv) --> false
    @fact <=(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(0.0, 2.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 2.0), trv) ≤ DecoratedInterval(Interval(0.0, 2.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-0.0, 2.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 2.0), trv) ≤ DecoratedInterval(Interval(-0.0, 2.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(0.0, 2.0), def), DecoratedInterval(Interval(1.0, 2.0), def)) --> true
    @fact DecoratedInterval(Interval(0.0, 2.0), def) ≤ DecoratedInterval(Interval(1.0, 2.0), def) --> true
    @fact <=(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> true
    @fact DecoratedInterval(Interval(-0.0, 2.0), trv) ≤ DecoratedInterval(Interval(1.0, 2.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ≤ DecoratedInterval(Interval(1.0, 2.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), def)) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) ≤ DecoratedInterval(Interval(3.0, 4.0), def) --> true
    @fact <=(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(1.0, 3.5), trv) ≤ DecoratedInterval(Interval(3.0, 4.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(1.0, 4.0), trv) ≤ DecoratedInterval(Interval(3.0, 4.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(-2.0, -1.0), trv), DecoratedInterval(Interval(-2.0, -1.0), trv)) --> true
    @fact DecoratedInterval(Interval(-2.0, -1.0), trv) ≤ DecoratedInterval(Interval(-2.0, -1.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(-3.0, -1.5), trv), DecoratedInterval(Interval(-2.0, -1.0), trv)) --> true
    @fact DecoratedInterval(Interval(-3.0, -1.5), trv) ≤ DecoratedInterval(Interval(-2.0, -1.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 0.0), trv) ≤ DecoratedInterval(Interval(-0.0, -0.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), def)) --> true
    @fact DecoratedInterval(Interval(-0.0, -0.0), trv) ≤ DecoratedInterval(Interval(0.0, 0.0), def) --> true
    @fact <=(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv)) --> true
    @fact DecoratedInterval(Interval(-0.0, 0.0), trv) ≤ DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, -0.0), trv)) --> true
    @fact DecoratedInterval(Interval(-0.0, 0.0), trv) ≤ DecoratedInterval(Interval(0.0, -0.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(0.0, -0.0), def), DecoratedInterval(Interval(0.0, 0.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, -0.0), def) ≤ DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact <=(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, -0.0), trv) ≤ DecoratedInterval(Interval(-0.0, 0.0), trv) --> true
end

facts("minimal_precedes_test") do
    @fact precedes(∅, Interval(3.0, 4.0)) --> true
    @fact precedes(Interval(3.0, 4.0), ∅) --> true
    @fact precedes(∅, ∅) --> true
    @fact precedes(Interval(1.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact precedes(Interval(0.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact precedes(Interval(-0.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact precedes(Interval(-Inf, Inf), Interval(1.0, 2.0)) --> false
    @fact precedes(Interval(-Inf, Inf), Interval(-Inf, Inf)) --> false
    @fact precedes(Interval(1.0, 2.0), Interval(3.0, 4.0)) --> true
    @fact precedes(Interval(1.0, 3.0), Interval(3.0, 4.0)) --> true
    @fact precedes(Interval(-3.0, -1.0), Interval(-1.0, 0.0)) --> true
    @fact precedes(Interval(-3.0, -1.0), Interval(-1.0, -0.0)) --> true
    @fact precedes(Interval(1.0, 3.5), Interval(3.0, 4.0)) --> false
    @fact precedes(Interval(1.0, 4.0), Interval(3.0, 4.0)) --> false
    @fact precedes(Interval(-3.0, -0.1), Interval(-1.0, 0.0)) --> false
    @fact precedes(Interval(0.0, 0.0), Interval(-0.0, -0.0)) --> true
    @fact precedes(Interval(-0.0, -0.0), Interval(0.0, 0.0)) --> true
    @fact precedes(Interval(-0.0, 0.0), Interval(0.0, 0.0)) --> true
    @fact precedes(Interval(-0.0, 0.0), Interval(0.0, -0.0)) --> true
    @fact precedes(Interval(0.0, -0.0), Interval(0.0, 0.0)) --> true
    @fact precedes(Interval(0.0, -0.0), Interval(-0.0, 0.0)) --> true
end

facts("minimal_precedes_dec_test") do
    @fact precedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), def)) --> true
    @fact precedes(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv)) --> true
    @fact precedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv)) --> true
    @fact precedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(1.0, 3.0), trv), DecoratedInterval(Interval(3.0, 4.0), def)) --> true
    @fact precedes(DecoratedInterval(Interval(-3.0, -1.0), def), DecoratedInterval(Interval(-1.0, 0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(-3.0, -1.0), trv), DecoratedInterval(Interval(-1.0, -0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(-3.0, -0.1), trv), DecoratedInterval(Interval(-1.0, 0.0), trv)) --> false
    @fact precedes(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), def)) --> true
    @fact precedes(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(-0.0, 0.0), def), DecoratedInterval(Interval(0.0, -0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv)) --> true
    @fact precedes(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv)) --> true
end

facts("minimal_interior_test") do
    @fact ∅ ⪽ ∅ --> true
    @fact interior(∅, ∅) --> true
    @fact ∅ ⪽ Interval(0.0, 4.0) --> true
    @fact interior(∅, Interval(0.0, 4.0)) --> true
    @fact Interval(0.0, 4.0) ⪽ ∅ --> false
    @fact interior(Interval(0.0, 4.0), ∅) --> false
    @fact Interval(-Inf, Inf) ⪽ Interval(-Inf, Inf) --> true
    @fact interior(Interval(-Inf, Inf), Interval(-Inf, Inf)) --> true
    @fact Interval(0.0, 4.0) ⪽ Interval(-Inf, Inf) --> true
    @fact interior(Interval(0.0, 4.0), Interval(-Inf, Inf)) --> true
    @fact ∅ ⪽ Interval(-Inf, Inf) --> true
    @fact interior(∅, Interval(-Inf, Inf)) --> true
    @fact Interval(-Inf, Inf) ⪽ Interval(0.0, 4.0) --> false
    @fact interior(Interval(-Inf, Inf), Interval(0.0, 4.0)) --> false
    @fact Interval(0.0, 4.0) ⪽ Interval(0.0, 4.0) --> false
    @fact interior(Interval(0.0, 4.0), Interval(0.0, 4.0)) --> false
    @fact Interval(1.0, 2.0) ⪽ Interval(0.0, 4.0) --> true
    @fact interior(Interval(1.0, 2.0), Interval(0.0, 4.0)) --> true
    @fact Interval(-2.0, 2.0) ⪽ Interval(-2.0, 4.0) --> false
    @fact interior(Interval(-2.0, 2.0), Interval(-2.0, 4.0)) --> false
    @fact Interval(-0.0, -0.0) ⪽ Interval(-2.0, 4.0) --> true
    @fact interior(Interval(-0.0, -0.0), Interval(-2.0, 4.0)) --> true
    @fact Interval(0.0, 0.0) ⪽ Interval(-2.0, 4.0) --> true
    @fact interior(Interval(0.0, 0.0), Interval(-2.0, 4.0)) --> true
    @fact Interval(0.0, 0.0) ⪽ Interval(-0.0, -0.0) --> false
    @fact interior(Interval(0.0, 0.0), Interval(-0.0, -0.0)) --> false
    @fact Interval(0.0, 4.4) ⪽ Interval(0.0, 4.0) --> false
    @fact interior(Interval(0.0, 4.4), Interval(0.0, 4.0)) --> false
    @fact Interval(-1.0, -1.0) ⪽ Interval(0.0, 4.0) --> false
    @fact interior(Interval(-1.0, -1.0), Interval(0.0, 4.0)) --> false
    @fact Interval(2.0, 2.0) ⪽ Interval(-2.0, -1.0) --> false
    @fact interior(Interval(2.0, 2.0), Interval(-2.0, -1.0)) --> false
end

facts("minimal_interior_dec_test") do
    @fact DecoratedInterval(∅, trv) ⪽ DecoratedInterval(Interval(0.0, 4.0), trv) --> true
    @fact interior(DecoratedInterval(∅, trv), DecoratedInterval(Interval(0.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 4.0), def) ⪽ DecoratedInterval(∅, trv) --> false
    @fact interior(DecoratedInterval(Interval(0.0, 4.0), def), DecoratedInterval(∅, trv)) --> false
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) ⪽ DecoratedInterval(∅, trv) --> false
    @fact interior(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(∅, trv)) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ⪽ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact interior(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) ⪽ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact interior(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> true
    @fact DecoratedInterval(∅, trv) ⪽ DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact interior(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> true
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) ⪽ DecoratedInterval(Interval(0.0, 4.0), trv) --> false
    @fact interior(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) --> false
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) ⪽ DecoratedInterval(Interval(0.0, 4.0), trv) --> false
    @fact interior(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), def) ⪽ DecoratedInterval(Interval(0.0, 4.0), trv) --> true
    @fact interior(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(0.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(-2.0, 2.0), trv) ⪽ DecoratedInterval(Interval(-2.0, 4.0), def) --> false
    @fact interior(DecoratedInterval(Interval(-2.0, 2.0), trv), DecoratedInterval(Interval(-2.0, 4.0), def)) --> false
    @fact DecoratedInterval(Interval(-0.0, -0.0), trv) ⪽ DecoratedInterval(Interval(-2.0, 4.0), trv) --> true
    @fact interior(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(-2.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 0.0), def) ⪽ DecoratedInterval(Interval(-2.0, 4.0), trv) --> true
    @fact interior(DecoratedInterval(Interval(0.0, 0.0), def), DecoratedInterval(Interval(-2.0, 4.0), trv)) --> true
    @fact DecoratedInterval(Interval(0.0, 0.0), trv) ⪽ DecoratedInterval(Interval(-0.0, -0.0), trv) --> false
    @fact interior(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
    @fact DecoratedInterval(Interval(0.0, 4.4), trv) ⪽ DecoratedInterval(Interval(0.0, 4.0), trv) --> false
    @fact interior(DecoratedInterval(Interval(0.0, 4.4), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) --> false
    @fact DecoratedInterval(Interval(-1.0, -1.0), trv) ⪽ DecoratedInterval(Interval(0.0, 4.0), def) --> false
    @fact interior(DecoratedInterval(Interval(-1.0, -1.0), trv), DecoratedInterval(Interval(0.0, 4.0), def)) --> false
    @fact DecoratedInterval(Interval(2.0, 2.0), def) ⪽ DecoratedInterval(Interval(-2.0, -1.0), trv) --> false
    @fact interior(DecoratedInterval(Interval(2.0, 2.0), def), DecoratedInterval(Interval(-2.0, -1.0), trv)) --> false
end

facts("minimal_strictLess_test") do
    @fact ∅ < ∅ --> true
    @fact Interval(1.0, 2.0) < ∅ --> false
    @fact ∅ < Interval(1.0, 2.0) --> false
    @fact Interval(-Inf, Inf) < Interval(-Inf, Inf) --> true
    @fact Interval(1.0, 2.0) < Interval(-Inf, Inf) --> false
    @fact Interval(-Inf, Inf) < Interval(1.0, 2.0) --> false
    @fact Interval(1.0, 2.0) < Interval(1.0, 2.0) --> false
    @fact Interval(1.0, 2.0) < Interval(3.0, 4.0) --> true
    @fact Interval(1.0, 3.5) < Interval(3.0, 4.0) --> true
    @fact Interval(1.0, 4.0) < Interval(3.0, 4.0) --> false
    @fact Interval(0.0, 4.0) < Interval(0.0, 4.0) --> false
    @fact Interval(-0.0, 4.0) < Interval(0.0, 4.0) --> false
    @fact Interval(-2.0, -1.0) < Interval(-2.0, -1.0) --> false
    @fact Interval(-3.0, -1.5) < Interval(-2.0, -1.0) --> true
end

facts("minimal_strictLess_dec_test") do
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) < DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(∅, trv) < DecoratedInterval(Interval(1.0, 2.0), def) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), def) < DecoratedInterval(∅, trv) --> false
    @fact DecoratedInterval(∅, trv) < DecoratedInterval(Interval(1.0, 2.0), def) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) < DecoratedInterval(Interval(-Inf, Inf), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) < DecoratedInterval(Interval(-Inf, Inf), trv) --> false
    @fact DecoratedInterval(Interval(-Inf, Inf), trv) < DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) < DecoratedInterval(Interval(1.0, 2.0), trv) --> false
    @fact DecoratedInterval(Interval(1.0, 2.0), trv) < DecoratedInterval(Interval(3.0, 4.0), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 3.5), def) < DecoratedInterval(Interval(3.0, 4.0), trv) --> true
    @fact DecoratedInterval(Interval(1.0, 4.0), trv) < DecoratedInterval(Interval(3.0, 4.0), def) --> false
    @fact DecoratedInterval(Interval(0.0, 4.0), trv) < DecoratedInterval(Interval(0.0, 4.0), def) --> false
    @fact DecoratedInterval(Interval(-0.0, 4.0), def) < DecoratedInterval(Interval(0.0, 4.0), trv) --> false
    @fact DecoratedInterval(Interval(-2.0, -1.0), def) < DecoratedInterval(Interval(-2.0, -1.0), def) --> false
    @fact DecoratedInterval(Interval(-3.0, -1.5), trv) < DecoratedInterval(Interval(-2.0, -1.0), trv) --> true
end

facts("minimal_strictPrecedes_test") do
    @fact strictprecedes(∅, Interval(3.0, 4.0)) --> true
    @fact strictprecedes(Interval(3.0, 4.0), ∅) --> true
    @fact strictprecedes(∅, ∅) --> true
    @fact strictprecedes(Interval(1.0, 2.0), Interval(-Inf, Inf)) --> false
    @fact strictprecedes(Interval(-Inf, Inf), Interval(1.0, 2.0)) --> false
    @fact strictprecedes(Interval(-Inf, Inf), Interval(-Inf, Inf)) --> false
    @fact strictprecedes(Interval(1.0, 2.0), Interval(3.0, 4.0)) --> true
    @fact strictprecedes(Interval(1.0, 3.0), Interval(3.0, 4.0)) --> false
    @fact strictprecedes(Interval(-3.0, -1.0), Interval(-1.0, 0.0)) --> false
    @fact strictprecedes(Interval(-3.0, -0.0), Interval(0.0, 1.0)) --> false
    @fact strictprecedes(Interval(-3.0, 0.0), Interval(-0.0, 1.0)) --> false
    @fact strictprecedes(Interval(1.0, 3.5), Interval(3.0, 4.0)) --> false
    @fact strictprecedes(Interval(1.0, 4.0), Interval(3.0, 4.0)) --> false
    @fact strictprecedes(Interval(-3.0, -0.1), Interval(-1.0, 0.0)) --> false
end

facts("minimal_strictPrecedes_dec_test") do
    @fact strictprecedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact strictprecedes(DecoratedInterval(Interval(3.0, 4.0), def), DecoratedInterval(∅, trv)) --> true
    @fact strictprecedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact strictprecedes(DecoratedInterval(Interval(3.0, 4.0), def), DecoratedInterval(∅, trv)) --> true
    @fact strictprecedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact strictprecedes(DecoratedInterval(Interval(1.0, 3.0), def), DecoratedInterval(Interval(3.0, 4.0), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-3.0, -1.0), trv), DecoratedInterval(Interval(-1.0, 0.0), def)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-3.0, -0.0), def), DecoratedInterval(Interval(0.0, 1.0), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-3.0, 0.0), trv), DecoratedInterval(Interval(-0.0, 1.0), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), def)) --> false
    @fact strictprecedes(DecoratedInterval(Interval(-3.0, -0.1), trv), DecoratedInterval(Interval(-1.0, 0.0), trv)) --> false
end

facts("minimal_disjoint_test") do
    @fact isdisjoint(∅, Interval(3.0, 4.0)) --> true
    @fact isdisjoint(Interval(3.0, 4.0), ∅) --> true
    @fact isdisjoint(∅, ∅) --> true
    @fact isdisjoint(Interval(3.0, 4.0), Interval(1.0, 2.0)) --> true
    @fact isdisjoint(Interval(0.0, 0.0), Interval(-0.0, -0.0)) --> false
    @fact isdisjoint(Interval(0.0, -0.0), Interval(-0.0, 0.0)) --> false
    @fact isdisjoint(Interval(3.0, 4.0), Interval(1.0, 7.0)) --> false
    @fact isdisjoint(Interval(3.0, 4.0), Interval(-Inf, Inf)) --> false
    @fact isdisjoint(Interval(-Inf, Inf), Interval(1.0, 7.0)) --> false
    @fact isdisjoint(Interval(-Inf, Inf), Interval(-Inf, Inf)) --> false
end

facts("minimal_disjoint_dec_test") do
    @fact isdisjoint(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), def)) --> true
    @fact isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv)) --> true
    @fact isdisjoint(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), trv)) --> true
    @fact isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv)) --> true
    @fact isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(Interval(1.0, 2.0), def)) --> true
    @fact isdisjoint(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
    @fact isdisjoint(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv)) --> false
    @fact isdisjoint(DecoratedInterval(Interval(3.0, 4.0), def), DecoratedInterval(Interval(1.0, 7.0), def)) --> false
    @fact isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
    @fact isdisjoint(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 7.0), trv)) --> false
    @fact isdisjoint(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) --> false
end
# FactCheck.exitstatus()

ext: jl


text:
using Test, MadNLP, MadNLPKrylov, MadNLPTests

testset = [
    [
        "Iterative",
        ()->MadNLP.Optimizer(
            iterative_solver=MadNLPKrylov,
            print_level=MadNLP.ERROR),
        []
    ],
]

@testset "MadNLPKrylov test" begin
    for (name,optimizer_constructor,exclude) in testset
        test_madnlp(name,optimizer_constructor,exclude)
    end
end


text:
# Author: Cedric St-Jean
# with some code taken from NBInclude.jl, by Steven G. Johnson
__precompile__()

module NBTesting

include("scrub_stderr.jl")

using JSON

export nbtest, is_testing, @testing_noop

is_skip(line) = startswith(line, "#NBSKIP") || startswith(line, "# NBSKIP")

const last_test_times = []
clear_test_times!() = empty!(last_test_times)
macro store_time(expr)
    esc(quote
        push!($NBTesting.last_test_times, ($(Expr(:quote, expr)), @elapsed($expr)))
        end)
end

const current_section = fill("")

"""     nbtranslate(path; outfile_name=..., verbose=5)

Translates the given .ipynb file into a .jl file, for testing.

 - All code following #NBSKIP within a cell will be ignored
 - `outfile_name` defaults to the name of the .ipynb file, preceded by `NBTest_`, and
   with a .jl extension.
 - The code will be wrapped inside a module called `NBTest_[Notebook name]`.
 - All headers that start with N pound signs (#) will be turned into print statements,
whenever `N <= verbose` (so the higher `verbose`, the more titles are printed)
"""
function nbtranslate(path::AbstractString;
                     outfile_name=joinpath(dirname(path),
                                           "NBTest_" * splitext(basename(path))[1]*".jl"),
                     verbose=5, time_cells=false)
    _, ext = splitext(path)
    @assert ext == ".ipynb" "nbtranslate only accepts notebook files (.ipynb)"

    # sleep a bit to process file requests from other nodes
    nprocs()>1 && sleep(0.005)
    nb = open(JSON.parse, path, "r")

    # check for an acceptable notebook:
    nb["nbformat"] == 4 || error("unrecognized notebook format ", nb["nbformat"])
    lang = lowercase(nb["metadata"]["language_info"]["name"])
    lang == "julia" || error("notebook is for unregognized language $lang")

    shell_or_help = r"^\s*[;?]" # pattern for shell command or help

    module_name = splitext(basename(outfile_name))[1]
    open(outfile_name, "w") do outfile
        write(outfile, string("module ", module_name, "\n"))
        write(outfile, "using NBTesting\n")  # to avoid ambiguity warning for is_testing
        write(outfile, "is_testing = true\n\n")
        write(outfile, "NBTesting.clear_test_times!()\n\n")
        counter = 0
        for cell in nb["cells"]
            if cell["cell_type"] == "code" && !isempty(cell["source"])
                counter += 1
                s = join(cell["source"])
                isempty(strip(s)) && continue # Jupyter doesn't number empty cells
                ismatch(shell_or_help, s) && continue
                cellnum = string(counter)

                lines = split(s, "\n")
                if !is_skip(lines[1])
                    if time_cells
                        lines = vcat(["nbtest_temp_time = Base.time_ns()"],
                                     lines,
                                     ["push!(NBTesting.last_test_times, ($counter, ((Base.time_ns)() - nbtest_temp_time)/1.e9))"])
                    end
                    write(outfile, string("# Cell #", cellnum, "\n"))
                    for line in lines
                        if is_skip(line) break end
                        write(outfile, line)
                        write(outfile, "\n")
                    end
                    write(outfile, "\n")
                end
            elseif cell["cell_type"] == "markdown" && !isempty(cell["source"])
                first_line = split(join(cell["source"]), "\n")[1]
                n_pound = findfirst(x->x!='#', first_line) - 1
                if 1 <= n_pound
                    if n_pound <= verbose
                        write(outfile, "# ----------------------------------- \n")
                        write(outfile, "println(\"$first_line\"); ")
                        # Since the point of verbosity is partly to report on progress,
                        # we flush STDOUT.
                        write(outfile, "flush(STDOUT)\n")
                    else
                        write(outfile, first_line)
                    end
                    write(outfile, "NBTesting.current_section[] = \"$first_line\"\n\n")
                end
            end
        end
        write(outfile, "end  # module")
    end
    return outfile_name
end

""" `is_testing` is true when called within `nbtest()`, and false otherwise. """
is_testing = false

"""     nbtest(path; outfile_name=..., verbose=5, keep_module=false, time_cells=false)

Translates the given .ipynb file into a .jl file for testing, then executes the file.

 - All code following #NBSKIP within a cell will be ignored
 - `outfile_name` defaults to the name of the .ipynb file, with a .jl extension.
 - NBTest wraps the notebook code inside a module called `NBTest_[Notebook name]`.
 - All headers that start with N pound signs (#) will be turned into print statements,
whenever `N <= verbose` (so the higher `verbose`, the more titles are printed)
 - If `time_cells=true`, each cell's running time is saved as `(Cell#, runtime)` in
`NBTesting.last_test_times`. See the generated .jl file to match the `Cell#`.
 - If `keep_module==true`, the second time that `nbtest(path)` is run for the same path,
the code will be evaluated in the same `NBTest_...` module instead of creating a new
module with the same name.
"""
function nbtest(path::AbstractString; verbose=5, keep_module=false, kwargs...)
    fname = nbtranslate(path; verbose=verbose, kwargs...)
    if verbose > 0
        info("Testing $path"); flush(STDERR) end
    current_section[] = ""  # just to make sure it's not misleading
    if keep_module
        return include_in_module(fname)
    else
        return include(fname)
    end
end

function include_in_module(nb_file::String)
    # I make a lot of effort to be able to run the code with `include`
    tmp_nb_file = splitext(nb_file)[1] * "_tmp.jl"
    local mod
    try
        mod_sym = open(tmp_nb_file, "w") do out
            println(out, "# Comment line so that the line numbers of this temp file lines up with the real one")
            lines = split(open(readstring, nb_file, "r"), '\n')
            for line in lines[2:end-1]
                write(out, line)
                write(out, "\n")
            end
            module_, mod_name = split(lines[1], " ")
            @assert module_ == "module"
            Symbol(mod_name)
        end
        if isdefined(Main, mod_sym)
            mod = eval(Main, mod_sym)
        else
            # module doesn't exist; just incude the file normally
            return include(nb_file)
        end
        eval(mod, :($NBTesting.sinclude($tmp_nb_file)))
        return mod
    finally
        rm(tmp_nb_file)
    end
end

noop(args...; kwargs...) = nothing

"""    @testing_noop fun1 fun2 ...

This macro doesn't do anything under normal execution, but when it is run by `nbtest`,
it turns the given function names into no-ops. It's primarily meant for disabling output
functions like `plot` during testing. For example:

    using Plots
    @testing_noop plot plot! vline!
"""
macro testing_noop(funs::Symbol...)
    esc(quote
        if is_testing
            $([:(@eval $f = $NBTesting.noop) for f in funs]...)
        end
    end)
end

end # module


text:
#================
Stats
    Module
================#
@doc """
      |====== Lathe.stats ======\n
      | ~~~~~~~~~~ Base ~~~~~~~~~~~\n
      |_____stats.mean(array)\n
      |_____stats.mode(array)\n
      |_____stats.variance(array)\n
      |_____stats.confiints(data,confidence = .95)\n
      |_____stats.standarderror(array)\n
      |_____stats.std(data)\n
      |_____stats.correlationcoeff(x,y)\n
      |_____stats.z(array)\n
      |_____stats.firstquar(array)\n
      |_____stats.secondquar(array)\n
      |_____stats.thirdquar(array)\n
      | ~~~~~~~~~~ Inferential ~~~~~~~~~~~\n
      |_____stats.independent_t(sample,general)\n
      |_____stats.spearman(var1,var2)\n
      |_____stats.pearson(x,y)\n
      |_____stats.chisqu(array)\n
      |_____stats.sign(array)\n
      |_____stats.f_test(sample,general)\n
      | ~~~~~~~~~~ Bayesian ~~~~~~~~~~~\n
      |_____stats.bay_ther(p,a,b)\n
      |_____stats.cond_prob(p,a,b)\n
      | ~~~~~~~~~~ Distributions ~~~~~~~~~~~\n
      |=== DISTRIBUTIONS HAVE NOW BEEN MOVED!\n
      |||||Now you can just use binomial_dist() or normal() for example\n
      |====== validate ======\n
      |____________/ Metrics ___________\n
      |_____validate.mae(actual,pred)\n
      |_____validate.r2(actual,pred)\n
       """
module stats
include("Distributions.jl")
#<----Mean---->
@doc """
      Calculates the mean of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      mean = Lathe.stats.mean(array)\n
      println(mean)\n
        10
       """
function mean(array)
    observations = length(array)
    average = sum(array)/observations
    return(average)
end
#<----Median---->
@doc """
      Calculates the median (numerical center) of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      median = Lathe.stats.median(array)\n
      println(median)\n
        10
       """
function median(array)
    n = length(array)
    half = n / 2
    current = 1
    sorted = sort!(array,rev = false)
    median = 0
    for i in sorted
        if current >= half
            median = i
        else
        current = current + 1
    end
    end
    return(median)
end
#<----Mode---->
@doc """
      Gives the digit most common in a given array\n
      --------------------\n
      array = [5,10,15,15,10,5,10]\n
      mode = Lathe.stats.mode(array)\n
      println(mode)\n
        10
       """
function mode(array)
    m = findmax(array)
    return(m)
end
#<----Variance---->
@doc """
      Gives the variance of an array..\n
      --------------------\n
      array = [5,10,15]\n
      variance = Lathe.stats.variance(array)\n
       """
function variance(array)
    me = mean(array)
    sq = sum(array) - me
    squared_mean = sq ^ 2
    return(squared_mean)
end
#<----Confidence Intervals---->
@doc """
      Returns the confidence intervals of given data.\n
      --------------------\n
      array = [5,10,15]\n
      confidence = .98\n
      low, high = Lathe.stats.confints(array,confidence)\n
      --------------------\n
      PARAMETERS\n
      confidence: Confidence is a float percentage representing the level of
      confidence required in your test. The confidence metric is used
      exclusively for calculating the interval.
       """
function confiints(data, confidence=.95)
    mean = mean(data)
    std = std(data)
    stderr = standarderror(data)
#    interval = stderr * scs.t.ppf((1 + confidence) / 2.0, n-1)
#    return (mean-interval, mean+interval)
end
#<----Standard Error---->
@doc """
      Calculates the Standard Error of an array.\n
      --------------------\n
      array = [5,10,15]\n
      ste = Lathe.stats.standarderror(array)\n
       """
function standarderror(data)
    std = std(data)
    sample = length(data)
    ste = (std/sqrt(sample))
    return(ste)
end
#<----Standard Deviation---->
@doc """
      Calculates the Standard Deviation of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      std = Lathe.stats.std(array)\n
       """
function std(array3)
    m = mean(array3)
    [i = (i-m) ^ 2 for i in array3]
    m = mean(array3)
    m = sqrt(Complex(m))
    return(m)
end
#<---- correlation Coefficient --->
@doc """
      Calculates the Correlation Coeffiecient of between two features\n
      --------------------\n
      x = [5,10,15]\n
      y = [5,10,15]\n
      r = Lathe.stats.correlationcoeff(x,y)\n
       """
function correlationcoeff(x,y)
    n = length(x)
    yl = length(y)
    if n != yl
        throw(ArgumentError("The array shape does not match!"))
    end
    xy = x .* y
    sx = sum(x)
    sy = sum(y)
    sxy = sum(xy)
    x2 = x .^ 2
    y2 = y .^ 2
    sx2 = sum(x2)
    sy2 = sum(y2)
    r = ((n*sxy) - (sx * sy)) / (sqrt((((n*sx2)-(sx^2)) * ((n*sy2)-(sy^2)))))
    return(r)
end
#<----Z score---->
@doc """
      Calculates the Z score of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      z = Lathe.stats.z(array)\n
       """
function z(array)
    x̄ = mean(array)
    σ = std(array)
    return map(x -> (x - x̄) / σ, array)
end
#<----Quartiles---->
# - First
@doc """
      Returns the point in an array located at 25 percent of the sorted data.\n
      --------------------\n
      array = [5,10,15]\n
      q1 = Lathe.stats.firstquar(array)\n
       """
function firstquar(array)
    m = median(array)
    q1 = array * .5
    return(q1)
end
# - Second(median)
@doc """
      Returns the point in an array located at 50 percent of the sorted data.
      The second quartile is also known as the median, or the middle of the sorted data.\n
      --------------------\n
      array = [5,10,15]\n
      q2 = Lathe.stats.secondquar(array)\n
       """
function secondquar(array)
    m = median(array)
    return(m)
end
# - Third
@doc """
      Returns the point in an array located at 75 percent of the sorted data.\n
      --------------------\n
      array = [5,10,15]\n
      q2 = Lathe.stats.secondquar(array)\n
       """
function thirdquar(array)
    q = median(array)
    q3 = q * 1.5
    return(q3)
end
# <---- Rank ---->
@doc """
      Ranks indices in an array based on quantitative weight (count of the
      numbers) and returns a new array of the ranks of each column. This
      function is made primarily for the Wilcox Rank-Sum test.\n
      --------------------\n
      array = [5,10,15]\n
      q2 = Lathe.stats.secondquar(array)\n
       """
function getranks(array,rev = false)
    sortedar = sort(array,rev=rev)
    num = 1
    list = []
    for i in sortedar
        append!(list,i)
        num = num + 1
    end
    return(list)
end
#-------Inferential-----------__________
#<----T Test---->
# - Independent
@doc """
      Performs an independent T test.\n
      --------------------\n
      sample = [5,10,15]
      general = [15,25,35]\n
      t = Lathe.stats.independent_t(sample,general)\n
       """
function independent_t(sample,general)
    sampmean = mean(sample)
    genmean = mean(general)
    samples = length(sample)
    m = genmean
    [i = (i-m) ^ 2 for i in general]
    m = mean(general)
    m = sqrt(m)
    std = m
    t = (sampmean - genmean) / (std / sqrt(samples))
    return(t)
end
# - Paired
#@doc """
#    THIS FUNCTION IS NOT YET WRITTEN\n
#      Paired T (Dependent T) is a T-test that doesn't require a sample.\n
#      --------------------\n
#      array = [5,10,15]\n
#      q2 = Lathe.stats.paired_t(var1,var2)\n
#       """ ->
function paired_t(var1,var2)
    d = var1 .- var2
    d̄ = mean(x)
end
#<---- Correlations ---->
# - Spearman
@doc """
      Returns a probability using a Spearman correlation.\n
      --------------------\n
      var1 = [5,10,15]\n
      var2 = [5,4,3,5,6]\n
      r = Lathe.stats.spearman(var1,var2)\n
       """ ->
function spearman(var1,var2)
    rgX = getranks(var1)
    rgY = getranks(var2)
    ρ = rgX*rgY / (std(rgX)*std(rgY))
    return(ρ)
end
# - Pearson
@doc """
      Returns a probability using a Pearson correlation.\n
      --------------------\n
      x = [5,10,15]\n
      y = [5,4,3,5,6]\n
      r = Lathe.stats.spearman(x,y)\n
       """
function pearson(x,y)
    sx = std(x)
    sy = std(y)
    x̄ = mean(x)
    ȳ = mean(x)
    [i = (i-x̄) / sx for i in x]
    [i = (i-ȳ) / sy for i in y]
    n1 = n-1
    mult = x .* y
    sq = sum(mult)
    corrcoff = sq / n1
    return(corrcoff)
end
#<---- Chi-Square ---->
@doc """
      Returns a probability using a chi squared distribution.\n
      --------------------\n
      var1 = [5,10,15]\n
      var2 = [5,4,3,5,6]\n
      p = Lathe.stats.chisq(var1,var2)\n
       """
function chisq(var1,var2)
    chistat(obs, exp) = (obs - exp)^2/exp
    return chistat.(x, e) |> sum
end
#<---- Wilcoxon ---->
# - Wilcoxon Rank-Sum Test
function wilcoxrs(var1,var2)
    #Hash
end
@doc """
      FUNCTION NOT YET WRITTEN\n
      Wilcox Sum Rank Tests are used to determine a probability with ranks\n
      --------------------\n
      var1 = [5,10,15]\n
      var2 = [19,25,30]\n
      p = Lathe.stats.wilcoxsr(var1,var2)\n
       """
function wilcoxsr(var1,var2)
    #Hash
end
#<---- Sign Test ---->
@doc """
      The Sign test determines correlation through negative and positive
      placement with binomial distribution.\n
      --------------------\n
      var1 = [5,10,15]\n
      var2 = [19,25,30]\n
      p = Lathe.stats.sign(var1,var2)\n
       """
function sign(var1,var2)
    sets = var1 .- var2
    positives = []
    negatives = []
    zeros = []
    for i in sets
        if i == 0
            zeros.append(i)
        elseif i > 0
            positives.append(i)
        elseif i < 0
            negatives.append(i)
        end
    end
    totalpos = length(positives)
    totallen = length(sets)
    ans = binomial_dist(positives,totallen)
    return(ans)
end
#<---- F-Test---->
@doc """
      An F test returns a probability of correlation, and is used similarly
      to a T test.\n
      --------------------\n
      array = [5,10,15]\n
      r = Lathe.stats.anova(array)\n
       """
function f_test(sample,general)
    totvariance = variance(general)
    sampvar = variance(sample)
    f =  sampvar / totvariance
    return(f)
end
#-------Bayesian--------------___________
#<----Bayes Theorem---->
#P = prob, A = prior, B = Evidence,
function bay_ther(p,a,b)
    psterior = (p*(b|a) * p*(a)) / (p*b)
    return(psterior)
end
@doc """
      Performs Bayesian Conditional Probability and returns probability.\n
      --------------------\n
      prob = .50\n
      prior = .20\n
      evidence = .30\n
      p = Lathe.stats.cond_prob(prob,prior,evidence)\n
       """
function cond_prob(p,a,b)
    psterior = bay_ther(p,a,b)
    cond = p*(a|b)
    return(cond)
end

#-------Model Metrics--------____________
## <---- Mean Absolute Error ---->
@doc """
      Mean absolute error (MAE) subtracts two arrays and averages the
      difference.\n
      --------------------\n
      array = [5,10,15]\n
      r = Lathe.stats.anova(array)\n
       """
function mae(actual,pred)
    l = length(actual)
    lp = length(pred)
    if l != lp
        throw(ArgumentError("The array shape does not match!"))
    end
    result = actual-pred
    maeunf = stats.mean(result)
    if maeunf < 0
        maeunf = (maeunf - maeunf) - maeunf
    end
    return(maeunf)
end
# <---- Mean Squared Error ---->
@doc """
      Mean Square error (MSE) subtracts two arrays, squares the
      difference, and averages the result\n
      --------------------\n
      array = [5,10,15]\n
      r = Lathe.stats.anova(array)\n
       """
function mse(y,ŷ)
    diff = y .- ŷ
    diff = diff .^ 2
    Σdiff = sum(diff)
    return(Σdiff)
end
# <---- R Squared ---->
@doc """
      R squared is the correlation coefficient of regression, and is found
      by squaring the correlation coefficient.\n
      --------------------\n
      array = [5,10,15]\n
      r = Lathe.stats.anova(array)\n
       """
function r2(actual,pred)
    l = length(actual)
    lp = length(pred)
    if l != lp
        throw(ArgumentError("The array shape does not match!"))
    end
    r = Lathe.stats.correlationcoeff(actual,pred)
    rsq = r^2
    rsq = rsq * 100
    return(rsq)
end
function catacc(yhat,testy)
    n = length(yhat)
    if n != length(testy)
        throw(ArgumentError("The array shape does not match!"))
    end
    score = 0
    for (i,w) in zip(yhat,testy)
        if i == w
            score += 1
        end
    end
    acc = score / n
end

#---------------------------
end


text:
######################
# Redundancy removal #
######################

# Redundancy
export detecthlinearity!, detectvlinearity!, dim
export isredundant, removevredundancy!, removehredundancy!, gethredundantindices, getvredundantindices

"""
    detecthlinearity!(p::VRep)

Detects all the hyperplanes contained in the H-representation and remove all redundant hyperplanes.

## Examples
The representation
```julia
h = HalfSpace([1, 1], 1]) ∩ HalfSpace([-1, -1], -1)
```
contains the hyperplane `HyperPlane([1, 1], 1)`.
"""
detecthlinearity!(p::HRep) = error("detecthlinearity! not implemented for $(typeof(p))")
detecthlinearity!(p::Polyhedron) = sethrep!(p, removeduplicates(hrep(p)))

"""
    detectvlinearity!(p::VRep)

Detects all the lines contained in the V-representation and remove all redundant lines.

## Examples
The representation
```julia
v = conichull([1, 1], [-1, -1])
```
contains the line `Line([1, 1])`.
"""
detectvlinearity!(p::VRep) = error("detectvlinearity! not implemented for $(typeof(p))")
detectvlinearity!(p::Polyhedron) = setvrep!(p, removeduplicates(vrep(p)))
function detectvlinearity!(p::VPolytope) end # No ray so no line

"""
    dim(h::HRep, current=false)

Returns the dimension of the affine hull of the polyhedron.
That is the number of non-redundant hyperplanes that define it.
If `current` is `true` then it simply returns the dimension according the current number of hyperplanes, assuming that the H-linearity has already been detected.
Otherwise, it first calls [`detecthlinearity!`](@ref).
"""
function dim(h::HRep, current=false)
    if !current
        detecthlinearity!(h)
    end
    fulldim(h) - nhyperplanes(h)
end

"""
    isredundant(p::Rep, idx::Index; strongly=false)

Return a `Bool` indicating whether the element with index `idx` can be removed without changing the polyhedron represented by `p`.
If `strongly` is `true`,
* if `idx` is an H-representation element `h`, it returns `true` only if no V-representation element of `p` is in the hyperplane of `h`.
* if `idx` is a V-representation element `v`, it returns `true` only if `v` is in the relative interior of `p`.
"""
function isredundant end

"""
    removehredundancy!(p::HRep)

Removes the elements of the H-representation of `p` that can be removed without changing the polyhedron represented by `p`. That is, it only keeps the halfspaces corresponding to facets of the polyhedron.
"""
function removehredundancy! end
function removehredundancy!(p::Polyhedron)
    detectvlinearity!(p)
    detecthlinearity!(p)
    sethrep!(p, removehredundancy(hrep(p), vrep(p)))
end

"""
    removevredundancy!(p::VRep)

Removes the elements of the V-representation of `p` that can be removed without changing the polyhedron represented by `p`. That is, it only keeps the extreme points and rays. This operation is often called "convex hull" as the remaining points are the extreme points of the convex hull of the initial set of points.
"""
function removevredundancy! end
function removevredundancy!(p::Polyhedron)
    detecthlinearity!(p)
    detectvlinearity!(p)
    setvrep!(p, removevredundancy(vrep(p), hrep(p)))
end

function _filter(f, it)
    # FIXME returns a Vector{Any}
    #collect(Iterators.filter(f, it)) # filter does not implement length so we need to collect
    ret = eltype(it)[]
    for el in it
        if f(el)
            push!(ret, el)
        end
    end
    ret
end
function removevredundancy(vrepit::VIt, hrep::HRep; nl=nlines(hrep), kws...)
    _filter(v -> !isredundant(hrep, v; nl=nl, kws...), vrepit)
end

# Remove redundancy in the V-representation using the H-representation
# There shouldn't be any duplicates in hrep for this to work
function removevredundancy(vrep::VRep, hrep::HRep; kws...)
    nl = nlines(vrep)
    typeof(vrep)(FullDim(vrep),
                 removevredundancy.(vreps(vrep), hrep; nl=nl, kws...)...)::typeof(vrep)
end

function removehredundancy(hrepit::HIt, vrep::VRep; strongly=true, d=dim(vrep))
    _filter(h -> !isredundant(vrep, h, strongly=strongly, d=d), hrepit)
end

# Remove redundancy in the H-representation using the V-representation
# There shouldn't be any duplicates in vrep for this to work
function removehredundancy(hrep::HRep, vrep::VRep; strongly=true)
    R = BitSet()
    d = dim(hrep, true) # TODO dim(hrep)
    typeof(hrep)(FullDim(hrep),
                 removehredundancy.(hreps(hrep), vrep,
                                    strongly=strongly, d=d)...)
end

#function gethredundantindices(hrep::HRep; strongly=false, solver=Polyhedra.solver(hrep))
#    red = BitSet()
#    for (i, h) in enumerate(hreps(hrep))
#        if ishredundant(hrep, h; strongly=strongly, solver=solver)
#            push!(red, i)
#        end
#    end
#    red
#end
#function getvredundantindices(vrep::VRep; strongly = true, solver=Polyhedra.solver(vrep))
#    red = BitSet()
#    for (i, v) in enumerate(vreps(vrep))
#        if isvredundant(vrep, v; strongly=strongly, solver=solver)
#            push!(red, i)
#        end
#    end
#    red
#end

# V-redundancy
# If p is an H-representation, nl needs to be given otherwise if p is a Polyhedron, it can be asked to p.
# TODO nlines should be the number of non-redundant lines so something similar to dim
function isredundant(p::HRep{T}, v::Union{AbstractVector, Line, Ray}; strongly = true, nl::Int=nlines(p), solver=Polyhedra.solver(p)) where {T}
    # v is in every hyperplane otherwise it would not be valid
    hcount = nhyperplanes(p) + count(h -> v in hyperplane(h), halfspaces(p))
    strong = (isray(v) ? fulldim(p)-1 : fulldim(p)) - nl
    hcount < (strongly ? strong : min(strong, 1))
end
# A line is never redundant but it can be a duplicate
isredundant(p::HRep{T}, v::Line; strongly = true, nl::Int=nlines(p), solver=Polyhedra.solver(p)) where {T} = false

# H-redundancy
# If p is a V-representation, nl needs to be given otherwise if p is a Polyhedron, it can be asked to p.
function isredundant(p::VRep{T}, h::HRepElement; strongly = true, d::Int=dim(p), solver=Polyhedra.solver(p)) where {T}
    checkvconsistency(p)
    hp = hyperplane(h)
    pcount = count(p -> p in hp, points(p))
    # every line is in h, otherwise it would not be valid
    rcount = nlines(p) + count(r -> r in hp, rays(p))
    pcount < min(d, 1) || (strongly && pcount + rcount < d)
end
# An hyperplane is never redundant but it can be a duplicate
isredundant(p::VRep{T}, h::HyperPlane; strongly = true, d::Int=dim(p), solver=Polyhedra.solver(p)) where {T} = false

# H-redundancy
#function ishredundantaux(p::HRep, a, β, strongly, solver)
#    sol = MPB.linprog(-a, p, solver)
#    if sol.status == :Unbounded
#        false
#    elseif sol.status == :Optimal
#        if _gt(sol.objval, β)
#            false
#        elseif _geq(sol.objval, β)
#            if strongly
#                false
#            else
#                true
#            end
#        else
#            true
#        end
#    end
#end
#function ishredundant(p::Rep, h::HRepElement; strongly = false, solver=Polyhedra.solver(p))
#    if islin(h)
#        sol = ishredundantaux(p, h.a, h.β, strongly, solver)
#        if !sol[1]
#            sol
#        else
#            ishredundantaux(p, -h.a, -h.β, strongly, solver)
#        end
#    else
#        ishredundantaux(p, h.a, h.β, strongly, solver)
#    end
#end

######################
# Duplicates removal #
######################
export removeduplicates

"""
    removeduplicates(rep::Representation)

Removes the duplicates in the Representation.

* In an H-representation, it removes the redundant hyperplanes and it remove an halfspace when it is equal to another halfspace in the affine hull.
  For instance, `HalfSpace([1, 1], 1)` is equal to `HalfSpace([1, 0], 0)` in the affine hull generated by `HyperPlane([0, 1], 1])`.
* In a V-representation, it removes the redundant lines and it remove a point (resp. ray) when it is equal to another point (resp. ray) in the line hull.
  For instance, in the line hull generated by `Line([0, 1])`, `[1, 1]` is equal to `[1, 0]` and `Ray([2, 2])` is equal to `Ray([1, 0])`.
"""
function removeduplicates end

# H-duplicates
# Separate function so that it is compiled with a concrete type for p
#function vpupdatedup!(aff, points, sympoints, p::SymPoint)
#    found = false
#    # sympoints are treated before points so there shouldn't be any
#    @assert isempty(points)
##    for (i, q) in enumerate(points)
##        if (coord(p) - q) in aff || (coord(p) + q) in aff
##            found = true
##            deleteat!(points, i)
##            push!(sympoints, p)
##            break
##        end
##    end
#    if !found && !any(sp -> (coord(sp) - coord(p)) in aff || (coord(sp) + coord(p)) in aff, sympoints)
#        push!(sympoints, p)
#    end
#end
function vpupdatedup!(aff, points, p::AbstractVector)
    if !any(point -> (point - p) in aff, points)
        push!(points, p)
    end
end
#function vrupdatedup!(rays, lines, l::Line)
#    if !any(isapprox.(lines, [l]))
#        for (i, r) in enumerate(rays)
#            if line(r) ≈ l
#                found = true
#                deleteat!(rays, i)
#                break
#            end
#        end
#        push!(lines, l)
#    end
#end
function vrupdatedup!(aff::VLinearSpace, rays::Vector{<:Ray}, r)
    r = remproj(r, aff)
    if !isapproxzero(r) && !any(ray -> remproj(ray, aff) ≈ r, rays)
        l = line(r)
        found = false
        for (i, s) in enumerate(rays)
            if line(s) ≈ l
                deleteat!(rays, i)
                found = true
                break
            end
        end
        if found
            convexhull!(aff, l)
        else
            push!(rays, r)
        end
        found
    else
        false
    end
end
function premovedups(vrep::VRepresentation, aff::VLinearSpace)
    ps = pointtype(vrep)[]
    for p in points(vrep)
        vpupdatedup!(aff, ps, p)
    end
    tuple(ps)
end
function removeduplicates(vrep::VPolytope)
    typeof(vrep)(FullDim(vrep), premovedups(vrep, emptyspace(vrep))...)
end
function removeduplicates(vrep::VRepresentation)
    aff = linespace(vrep, true)
    newlin = true
    rs = raytype(vrep)[]
    while newlin
        newlin = false
        empty!(rs)
        for r in rays(vrep)
            newlin |= vrupdatedup!(aff, rs, r)
        end
    end
    typeof(vrep)(FullDim(vrep), premovedups(vrep, aff)..., aff.lines, rs)
end

# H-duplicates
#function hupdatedup!(hp, hs, h::HyperPlane)
#    if !any(isapprox.(hp, [h]))
#        for (i, s) in enumerate(hs)
#            if hyperplane(s) ≈ h
#                deleteat!(hs, i)
#                break # There should be no duplicate in hp so no need to continue
#            end
#        end
#        push!(hp, h)
#    end
#end
function hupdatedup!(aff::HAffineSpace, hss, h::HalfSpace)
    h = remproj(h, aff)
    if !isapproxzero(h) && !any(hs -> _isapprox(remproj(hs, aff), h), hss)
        hp = hyperplane(h)
        found = false
        for (i, hs) in enumerate(hss)
            # TODO Not enough, e.g.
            # x <= 1
            # y <= 1
            # x + y >= 2
            if hyperplane(hs) ≈ hp
                deleteat!(hss, i)
                found = true
                break # There should be no duplicate in hp so no need to continue
            end
        end
        if found
            intersect!(aff, hp)
        else
            push!(hss, h)
        end
        found
    else
        false
    end
end
function removeduplicates(hrep::HRepresentation{T}) where {T}
    aff = affinehull(hrep, true)
    newlin = true
    hs = halfspacetype(hrep)[]
    while newlin
        newlin = false
        aff = removeduplicates(aff)
        empty!(hs)
        for h in halfspaces(hrep)
            newlin |= hupdatedup!(aff, hs, h)
        end
    end
    typeof(hrep)(FullDim(hrep), aff.hyperplanes, hs)
end


text:
# Sequentialize.jl
#
# Sequentialization for recurrent networks
# The current Flux.jl setup does not allow to feed sequences as a matrix to a recurrent
# network. The constructions in this file help to unwrap a matrix to a sequence, apply
# an RNN and then wrap things up again.

"""
    HiddenRecur(cell)

`HiddenRecur` takes a recurrent cell and makes it stateful, managing the hidden state
in the background. As opposed to `Recur`, it returns the both the hidden state and the
cell state.

See also: [`Flux.Recur`](@ref)
"""
mutable struct HiddenRecur{T,S}
    cell::T
    state::S
end

HiddenRecur(m) = HiddenRecur(m, m.state0)

function (m::HiddenRecur)(xs...)
    h, y = m.cell(m.state, xs...)
    m.state = h
    return h  # collect(h)
end

Flux.@functor HiddenRecur
trainable(a::HiddenRecur) = (a.cell,)

Base.show(io::IO, m::HiddenRecur) = print(io, "HiddenRecur(", m.cell, ")")
Flux.reset!(m::HiddenRecur) = (m.state = m.cell.state0)

"""
    Seq(RNN)

`Seq` takes a recurrent neural network and "sequentializes" it, i.e. when `Seq(RNN)` is
called with a matrix of input features over a certain time interval, the recurrent neural
net is fed with a sequence of inputs, and results are transformed back to matrix form.
"""
mutable struct Seq{T}
    chain::T
    state
end

Seq(chain) = Seq(chain, [0.0f0])
(l::Seq)(x) = l(l.chain.state, l.chain, x)
getbuffersize(chain, x) = getbuffersize(typeof(chain), chain.state, x)

# CPU Arrays
function (l::Seq)(::Array, _, x)
	l.state = Align(map(l.chain, Slices(x, True(), False())), 1)
	return l.state
end
function (l::Seq)(::Array, ::Flux.Recur, x)
	l.state = Align(map(l.chain, Slices(x, True(), False())), 1)
	return l.state
end
function (l::Seq)(::Tuple, ::Flux.Recur, x)
	l.state = Align(map(l.chain, Slices(x, True(), False())), 1)
	return l.state
end
function (l::Seq)(::Tuple, _, x)
	tuples = map(l.chain, Slices(x, True(), False()))
	l.state = [Align(map(x -> dropdims(x[i], dims=2), tuples), 1) for i in 1:length(l.chain.state)]
	return l.state
end

# GPU Arrays
getbuffersize(::Type{<:Union{Flux.Recur,StackedLSTMCell}},
 				state::AbstractArray, x) = ((length(state), size(x, 2)), 1)
getbuffersize(::Type{<:Union{Flux.Recur,StackedLSTMCell}},
				state::Tuple, x) = ((length(state[1]), size(x, 2)), 1)
getbuffersize(::Type{<:Union{HiddenRecur}},
				state::Tuple, x) = ((length(state[1]), size(x, 2)), length(state))

function writebuffer(chain::HiddenRecur, x)
	sizeval, numhidden = getbuffersize(chain, x)
	out = Flux.Zygote.Buffer(x, numhidden * sizeval[1], sizeval[2])
	for i = 1:sizeval[2]
	  out[:,i] = cat(chain(x[:,i])...; dims=1)
	end
	output = copy(out)
	return [output[i * sizeval[1] + 1:(i + 1) * sizeval[1],:] for i = 0:numhidden - 1]
end

function writebuffer(chain, x)
	sizeval, numhidden = getbuffersize(chain, x)
	out = Flux.Zygote.Buffer(x, sizeval)
	for i = 1:sizeval[2]
	  out[:,i] = chain(x[:,i])
	end
	return copy(out)
end

function (l::Seq)(::Flux.CUDA.CuArray, _, x)
	l.state = writebuffer(l.chain, x)
	return l.state
end

function (l::Seq)(::Tuple{Flux.CUDA.CuArray,Flux.CUDA.CuArray}, _, x)
	l.state = writebuffer(l.chain, x)
	return l.state
end

function Base.show(io::IO, l::Seq)
    print(io, "Seq(", l.chain)
    print(io, ")")
end

function Flux.reset!(m::Seq)
    Flux.reset!(m.chain)
    # m.state = [0.0f0]  # this return value makes Zygote crash - why?
end

Flux.@functor Seq
Flux.trainable(m::Seq) = Flux.trainable(m.chain)

"""
    SeqSkip(RNNCell, skiplength::Integer)

`SeqSkip` takes a recurrent neural network cell and "sequentializes" it, i.e. when
it is called with a matrix of input features over a certain time interval, the recurrent
neural net is fed with a sequence of inputs, and results are transformed back to matrix
form. In addition, the hidden state from `skiplength` timesteps ago is used instead of the
current hidden state. This structure combines functionality of `Recur` in that it makes a
recurrent neural network cell stateful, as well as `Seq` in that it feeds matrices of
input features as elements of a time series.

See also: [`Seq`](@ref), [`Flux.Recur`](@ref)
"""
mutable struct SeqSkip{T,P}
	cell::T
	p::P  # skiplength
	state  # current cell state
	fullstate  # cell states over time series
end

SeqSkip(m, p) = SeqSkip(m, p, m.state0, m.state0)

getbuffersize(::Type{<:SeqSkip}, state::AbstractArray, x) = ((length(state), size(x, 2)), 1)
getbuffersize(::Type{<:SeqSkip}, state::Tuple, x) = ((length(state[1]), size(x, 2)), length(state))

writebuffer(l::SeqSkip,x) = writebuffer(typeof(l.state), l, x)

function writebuffer(::Type{<:Union{AbstractArray}}, l::SeqSkip, x)
	sizeval, numhidden = getbuffersize(l, x)
	out = Flux.Zygote.Buffer(x, sizeval)
	h, y = l.cell(l.state, x[:,1])
	out[:,1] = h
	l.state = h
	for i = 2:sizeval[2]
	  h, y = l.cell(out[:,max(1, i - l.p)], x[:,i])
	  l.state = h
	  out[:,i] = h
	end
	l.fullstate = copy(out)
	return l.fullstate
end

function writebuffer(::Type{<:Tuple}, l::SeqSkip, x)
	sizeval, numhidden = getbuffersize(l, x)
	out = Flux.Zygote.Buffer(x, numhidden * sizeval[1], sizeval[2])
	h, y = l.cell(l.state, x[:,1])
	out[:,1] = vcat(h...)
	l.state = h
	for i = 2:sizeval[2]
		hidden = Tuple([out[j * sizeval[1] + 1:(j + 1) * sizeval[1],max(1, i - l.p)] for j = 0:numhidden - 1 ])
		h, y = l.cell(hidden, x[:,i])
		l.state = h
		out[:,i] = vcat(h...)
	end
	l.fullstate = copy(out)
	return l.fullstate[1:sizeval[1],:]
end
# TO DO: For HiddenRecur-type we need writebuffer that outputs `output`

(l::SeqSkip)(x) = writebuffer(l, x)

Flux.@functor SeqSkip

Base.show(io::IO, m::SeqSkip) = print(io, "SeqSkip(", m.cell, ")")
Flux.reset!(m::SeqSkip) = (m.state = m.cell.state0)

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
export AdjointOperator

"""
`AdjointOperator(A::AbstractOperator)`

Shorthand constructor:

`'(A::AbstractOperator)`

Returns the adjoint operator of `A`.

```julia
julia> AdjointOperator(DFT(10))
ℱᵃ  ℂ^10 -> ℝ^10

julia> [DFT(10); DCT(10)]'
[ℱ;ℱc]ᵃ  ℂ^10  ℝ^10 -> ℝ^10
```
"""
struct AdjointOperator{T <: AbstractOperator} <: AbstractOperator
	A::T
	function AdjointOperator(A::T) where {T<:AbstractOperator}
		is_linear(A) == false && error("Cannot transpose a nonlinear operator. You might use `jacobian`")
		new{T}(A)
	end
end

# Constructors

AdjointOperator(L::AdjointOperator) = L.A

# Properties

size(L::AdjointOperator) = size(L.A,2), size(L.A,1)

domainType(L::AdjointOperator) = codomainType(L.A)
codomainType(L::AdjointOperator) = domainType(L.A)

fun_name(L::AdjointOperator)  = fun_name(L.A)*"ᵃ"

is_linear(L::AdjointOperator) = is_linear(L.A)
is_null(L::AdjointOperator) = is_null(L.A)
is_eye(L::AdjointOperator) = is_eye(L.A)
is_diagonal(L::AdjointOperator) = is_diagonal(L.A)
is_AcA_diagonal(L::AdjointOperator) = is_AAc_diagonal(L.A)
is_AAc_diagonal(L::AdjointOperator) = is_AcA_diagonal(L.A)
is_orthogonal(L::AdjointOperator) = is_orthogonal(L.A)
is_invertible(L::AdjointOperator) = is_invertible(L.A)
is_full_row_rank(L::AdjointOperator) = is_full_column_rank(L.A)
is_full_column_rank(L::AdjointOperator) = is_full_row_rank(L.A)

diag(L::AdjointOperator) = diag(L.A)
diag_AcA(L::AdjointOperator) = diag_AAc(L.A)
diag_AAc(L::AdjointOperator) = diag_AcA(L.A)


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    RegularGrid(dims, origin, spacing)

A regular grid with dimensions `dims`, lower left corner at `origin`
and cell spacing `spacing`. The three arguments must have the same length.

    RegularGrid(start, finish, dims=dims)

Alternatively, construct a regular grid from a `start` point (lower left)
to a `finish` point (upper right), optionally passing `dims`.

    RegularGrid{T}(dims)
    RegularGrid{T}(dim1, dim2, ...)

Finally, a regular grid can be constructed by only passing the dimensions
`dims` as a tuple, or by passing each dimension `dim1`, `dim2`, ... separately.
In this case, the origin and spacing default to (0,0,...) and (1,1,...).

## Examples

Create a 3D grid with 100x100x50 locations:

```julia
julia> RegularGrid{Float64}(100,100,50)
```

Create a 2D grid with 100x100 locations and origin at (10.,20.) units:

```julia
julia> RegularGrid((100,100),(10.,20.),(1.,1.))
```

Create a 1D grid from -1 to 1 with 100 locations:

```julia
julia> RegularGrid((-1.,),(1.,), dims=(100,))
```
"""
struct RegularGrid{T<:Real,N} <: AbstractDomain{T,N}
  dims::Dims{N}
  origin::NTuple{N,T}
  spacing::NTuple{N,T}

  function RegularGrid{T,N}(dims, origin, spacing) where {N,T<:Real}
    @assert all(dims .> 0) "dimensions must be positive"
    @assert all(spacing .> 0) "spacing must be positive"
    new(dims, origin, spacing)
  end
end

RegularGrid(dims::Dims{N}, origin::NTuple{N,T}, spacing::NTuple{N,T}) where {N,T<:Real} =
  RegularGrid{T,N}(dims, origin, spacing)

RegularGrid(start::NTuple{N,T}, finish::NTuple{N,T}; dims::Dims{N}=ntuple(i -> 100, N)) where {N,T<:Real} =
  RegularGrid{T,N}(dims, start, ntuple(i -> (finish[i]-start[i])/(dims[i]-1), N))

RegularGrid{T}(dims::Dims{N}) where {N,T<:Real} =
  RegularGrid{T,N}(dims, (zeros(T,N)...,), (ones(T,N)...,))

RegularGrid{T}(dims::Vararg{Int,N}) where {N,T<:Real} = RegularGrid{T}(dims)

npoints(grid::RegularGrid) = prod(grid.dims)

Base.size(grid::RegularGrid) = grid.dims
origin(grid::RegularGrid) = grid.origin
spacing(grid::RegularGrid) = grid.spacing

function coordinates!(buff::AbstractVector{T}, grid::RegularGrid{T,N},
                      location::Int) where {N,T<:Real}
  intcoords = CartesianIndices(grid.dims)[location]
  for i in 1:N
    @inbounds buff[i] = grid.origin[i] + (intcoords[i] - 1)*grid.spacing[i]
  end
end

function nearestlocation(grid::RegularGrid{T,N}, coords::AbstractVector{T}) where {N,T<:Real}
  sz = size(grid)
  or = origin(grid)
  sp = spacing(grid)

  units = ntuple(i -> @inbounds(return round(Int, (coords[i] - or[i]) / sp[i])), N)
  icoords = ntuple(i -> @inbounds(return units[i] + 1), N) # 1-based indexing

  # make sure integer coordinates lie inside of the grid
  c = ntuple(i -> @inbounds(return clamp(icoords[i], 1, sz[i])), N)

  if N == 1
    @inbounds return c[1]
  elseif N == 2
    @inbounds return c[1] + sz[1]*(c[2]-1)
  elseif N == 3
    @inbounds return c[1] + sz[1]*(c[2]-1) + sz[1]*sz[2]*(c[3]-1)
  else # higher dimensions
    LinearIndices(sz)[c...]
  end
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, grid::RegularGrid{T,N}) where {N,T<:Real}
  dims = join(grid.dims, "×")
  print(io, "$dims RegularGrid{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", grid::RegularGrid{T,N}) where {N,T<:Real}
  println(io, "RegularGrid{$T,$N}")
  println(io, "  dimensions: ", grid.dims)
  println(io, "  origin:     ", grid.origin)
  print(  io, "  spacing:    ", grid.spacing)
end


text:
# Plots an random example of a prior and the resulting information distribution

using SubmodularMaximization
using PyPlot

close("all")

steps = 100
grid_size = 10
num_observations = 3

grid = Grid(grid_size, grid_size)
sensor = RangingSensor()

target_state = random_state(grid)

# precomputation
histogram_filter = Filter(grid, target_state)

# Create some fake observations and update the filter
for ii = 1:num_observations
  robot_state = random_state(grid)
  range_observation = generate_observation(grid, sensor, robot_state,
                                           target_state)
  measurement_update!(histogram_filter, robot_state, get_states(grid), sensor,
                      grid, range_observation)
end

@show entropy(histogram_filter)

# Visualize the new prior
visualize_filter(histogram_filter, show_colorbar = false)
title("Prior")

# Information gain on the prior
information = Array(get_data(histogram_filter))
for ii = 1:length(information)
  state = SubmodularMaximization.index_to_state(grid, ii)

  information[ii] = finite_horizon_information(grid, histogram_filter, sensor,
                                               [state]).reward
end

# Visualize the information gain
figure()
limits = [0.5, grid_size + 0.5, 0.5, grid_size + 0.5]
image = imshow(information', cmap="viridis", vmin=0.0,
               extent=limits,
               interpolation="nearest", origin="lower")
#make_tight_colorbar(image)
title("Information gain")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
function HDF5.hdf5_type_id(::Type{T}) where T <: NamedTuple
    type_id = h5t_create(H5T_COMPOUND, sizeof(T))
    for i in 1:fieldcount(T)
        fname = string(fieldname(T, i))
        offset = fieldoffset(T, i)
        ftype = fieldtype(T, i)
        fid = hdf5_type_id(ftype)
        h5t_insert(type_id, fname, offset, fid)
        ftype <: MaxLenString && h5t_close(fid)
    end
    type_id
end

Base.zero(::Type{T}) where T <: NamedTuple = T(Base.zero(fieldtype(T, s)) for s in fieldnames(T))
Base.zero(::T) where T <: NamedTuple = T(Base.zero(fieldtype(T, s)) for s in fieldnames(T))

@h5bitslike NamedTuple