text:
using Test

function test_geometry()
    println(extents(Bounds(SA_F64[0, -1, 0], SA_F64[1, 2, 3])))
    true
end

function test_ideal_gas_eos()
    physical_eos = @show AdiabaticIdealGasEOS(
        0.1,
        PHYSICAL;
        dim = DIM1,
        atomicity = MONATOMIC,
        units = SI,
    )
    norm_eos = @show AdiabaticIdealGasEOS(2.1, NORMALIZED)
    true
end

function test_gaussian_kernel()
    W1 = @show GaussianKernel{1}()
    W2 = @show GaussianKernel{2}()
    W3 = @show GaussianKernel{3}()
    q, h = 0.7, 0.2
    @show W1(q, h)
    @show W2(q, h)
    @show W3(q, h)
    true
end

function test_initial_piecewise_uniform_mass_distr()
    @show PiecewiseUniformMassDistribution(Bounds(0, 3), [1, 2], [1, 2, 1], 10)()
    true
end

function test_initial_static_velocity_distr()
    @show StaticVelocityDistribution{2}()(
        fill(zeros(SVectorF{2}), 10),
        fill(1.0, 10),
    )
    true
end

function test_initial_uniform_energy_distr()
    @show UniformEnergyDistribution{2}(0.3)(
        fill(zeros(SVectorF{2}), 10),
        fill(1.0, 10),
        AdiabaticIdealGasEOS(2.0, NORMALIZED),
    )
    true
end

function test_initial_piecewise_uniform_energy_distr()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )
    data = mass_distribution()
    @show energy_distribution(
        data.positions,
        data.mass_densities,
        AdiabaticIdealGasEOS(2.1, NORMALIZED),
    )
    true
end

function test_initial_isothermal_energy_distr()
    @show IsothermalEnergyDistribution{2}(0.3)(
        fill(zeros(SVectorF{2}), 10),
        fill(1.0, 10),
        AdiabaticIdealGasEOS(2.5, NORMALIZED),
    )
    true
end

function test_velocity()
    positions = fill(zeros(SVectorF{2}), 10)
    mass_densities = fill(1.0, 10)
    velocity_distribution = StaticVelocityDistribution{2}()
    velocities = @show Velocities(
        positions,
        mass_densities,
        velocity_distribution,
    )
    accelerations = @show initderivatives(velocities)
    accelerations[2] = SA_F64[10, 5]
    initderivatives!(accelerations, velocities)
    @show accelerations
    true
end

function test_adaptive_mass_component()
    mass_distribution =
        PiecewiseUniformMassDistribution(Bounds(0, 3), [1, 2], [1, 2, 1], 10)
    kernel = GaussianKernel{1}()
    mass_component, positions = @show AdaptiveMass(mass_distribution, kernel)
    @show get_kernel_widths(mass_component)
    updatemasses!(mass_component, positions, kernel)
    @show get_kernel_widths(mass_component)
    true
end

function test_standard_energy_component()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    mass_component, positions =
        AdaptiveMass(mass_distribution, GaussianKernel{1}())
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )
    energy_component = @show StandardEnergy(
        positions,
        get_mass_densities(mass_component),
        energy_distribution,
        AdiabaticIdealGasEOS(2.5, NORMALIZED),
    )
    dudt = @show initderivatives(energy_component)
    dudt[2] = 5.0
    initderivatives!(dudt, energy_component)
    @show dudt
    true
end

function test_isothermal_energy_component()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    mass_component, positions =
        AdaptiveMass(mass_distribution, GaussianKernel{1}())
    energy_distribution = IsothermalEnergyDistribution{1}(4.5)
    energy_component = @show IsothermalEnergy(
        positions,
        get_mass_densities(mass_component),
        energy_distribution,
        AdiabaticIdealGasEOS(2.5, NORMALIZED),
    )
    true
end

function test_pressure()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    eos = AdiabaticIdealGasEOS(2.5, NORMALIZED)
    mass_component, positions =
        AdaptiveMass(mass_distribution, GaussianKernel{1}())
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )
    energy_component = StandardEnergy(
        positions,
        get_mass_densities(mass_component),
        energy_distribution,
        eos,
    )
    pressures = @show Pressures(mass_component, energy_component, eos)
    updatepressures!(pressures, mass_component, energy_component, eos)
    @show pressures
    true
end

function test_standard_diffusion_component()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    eos = AdiabaticIdealGasEOS(2.5, NORMALIZED)
    mass_component, positions =
        AdaptiveMass(mass_distribution, GaussianKernel{1}())
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )
    energy_component = StandardEnergy(
        positions,
        get_mass_densities(mass_component),
        energy_distribution,
        eos,
    )
    diffusion_component = @show StandardDiffusion(α = 1.0, β = 2.0, ε = 1e-3)
    updatediffusion!(diffusion_component, energy_component, eos)
    @show diffusion_component
    true
end

function test_euler_stepper()
    exterior_boundary_positions = Bounds(0, 3)
    interior_boundary_positions = [1, 2]
    mass_densities = [1, 2, 1]
    specific_energies = [0.1, 0.2, 0.1]
    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        10,
    )
    eos = AdiabaticIdealGasEOS(2.5, NORMALIZED)
    W = GaussianKernel{1}()
    mass_component, positions = AdaptiveMass(mass_distribution, W)
    mass_densities = get_mass_densities(mass_component)
    velocity_distribution = StaticVelocityDistribution{1}()
    velocities = Velocities(positions, mass_densities, velocity_distribution)
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )
    energy_component =
        StandardEnergy(positions, mass_densities, energy_distribution, eos)
    diffusion_component = StandardDiffusion()
    stepper = @show EulerStepper(
        positions,
        velocities,
        mass_component,
        energy_component,
        diffusion_component,
        W,
        eos,
    )
    step!(stepper, 0.1)
    @show stepper
    true
end

function test_simulation_builder()
    builder = @show SimulationBuilder(
        UniformMassDistribution(Bounds(SA_F64[0, 0], SA_F64[1, 1]), 1.0, 3^2),
        UniformEnergyDistribution{2}(1),
    ) + NoDiffusion
    @show build(builder)
    true
end

function runtests()

    @testset "geometry.jl" begin
        @test test_geometry()
    end

    @testset "eos.jl" begin
        @test test_ideal_gas_eos()
    end

    @testset "kernels.jl" begin
        @test test_gaussian_kernel()
    end

    @testset "initial/mass/piecewise_uniform.jl" begin
        @test test_initial_piecewise_uniform_mass_distr()
    end

    @testset "initial/velocity/static.jl" begin
        @test test_initial_static_velocity_distr()
    end

    @testset "initial/energy/uniform.jl" begin
        @test test_initial_uniform_energy_distr()
    end

    @testset "initial/energy/piecewise_uniform.jl" begin
        @test test_initial_piecewise_uniform_energy_distr()
    end

    @testset "initial/energy/isothermal.jl" begin
        @test test_initial_isothermal_energy_distr()
    end

    @testset "velocity.jl" begin
        @test test_velocity()
    end

    @testset "mass/adaptive.jl" begin
        @test test_adaptive_mass_component()
    end

    @testset "energy/standard.jl" begin
        @test test_standard_energy_component()
    end

    @testset "energy/isothermal.jl" begin
        @test test_isothermal_energy_component()
    end

    @testset "pressure.jl" begin
        @test test_pressure()
    end

    @testset "diffusion/standard.jl" begin
        @test test_standard_diffusion_component()
    end

    @testset "stepping/euler.jl" begin
        @test test_euler_stepper()
    end

    @testset "simulation.jl" begin
        @test test_simulation_builder()
    end
end
export runtests


text:
using IntervalSets
using Test
using Dates
using Statistics
import Statistics: mean
using Random

import IntervalSets: Domain, endpoints, closedendpoints, TypedEndpointsInterval

struct MyClosedUnitInterval <: TypedEndpointsInterval{:closed,:closed,Int} end
endpoints(::MyClosedUnitInterval) = (0,1)
Base.promote_rule(::Type{MyClosedUnitInterval}, ::Type{ClosedInterval{T}}) where T =
        ClosedInterval{T}

struct MyUnitInterval <: AbstractInterval{Int}
    isleftclosed::Bool
    isrightclosed::Bool
end
endpoints(::MyUnitInterval) = (0,1)
closedendpoints(I::MyUnitInterval) = (I.isleftclosed,I.isrightclosed)

struct IncompleteInterval <: AbstractInterval{Int} end

@testset "IntervalSets" begin
    if VERSION >= v"1.1"
        # Julia 1.0 defines getindex(a::GenericArray, i...) in Test,
        # which could cause an ambiguity with getindex(A::AbstractArray, ::EllipsisNotation.Ellipsis)
        @test isempty(detect_ambiguities(IntervalSets, Base, Core))
    end

    @test ordered(2, 1) == (1, 2)
    @test ordered(1, 2) == (1, 2)
    @test ordered(Float16(1), 2) == (1, 2)

    @testset "Basic Closed Sets" begin
        @test_throws ArgumentError :a .. "b"
        I = 0..3
        @test I === ClosedInterval(0,3) === ClosedInterval{Int}(0,3) ===
                 Interval(0,3)
        @test string(I) == "0..3"
        @test @inferred(UnitRange(I)) === 0:3
        @test @inferred(range(I)) === 0:3
        @test @inferred(UnitRange{Int16}(I)) === Int16(0):Int16(3)
        @test @inferred(ClosedInterval(0:3)) === I
        @test @inferred(ClosedInterval{Float64}(0:3)) === 0.0..3.0
        @test @inferred(ClosedInterval(Base.OneTo(3))) === 1..3
        J = 3..2
        K = 5..4
        L = 3 ± 2
        M = @inferred(ClosedInterval(2, 5.0))
        @test string(M) == "2.0..5.0"
        N = @inferred(ClosedInterval(UInt8(255), 300))
        O = @inferred(CartesianIndex(1, 2, 3, 4) ± 2)
        @test O == (-1..3, 0..4, 1..5, 2..6)

        @test eltype(I) == Int
        @test eltype(M) == Float64

        @test !isempty(I)
        @test isempty(J)
        @test J == K
        @test I != K
        @test I == I
        @test J == J
        @test L == L
        @test isequal(I, I)
        @test isequal(J, K)

        @test typeof(leftendpoint(M)) == typeof(rightendpoint(M)) && typeof(leftendpoint(M)) == Float64
        @test typeof(leftendpoint(N)) == typeof(rightendpoint(N)) && typeof(leftendpoint(N)) == Int
        @test @inferred(endpoints(M)) === (2.0,5.0)
        @test @inferred(endpoints(N)) === (255,300)

        @test maximum(I) === 3
        @test minimum(I) === 0
        @test extrema(I) === (0, 3)

        @test 2 in I
        @test issubset(1..2, 0.5..2.5)

        @test @inferred(I ∪ L) == ClosedInterval(0, 5)
        @test @inferred(I ∩ L) == ClosedInterval(1, 3)
        @test isempty(J ∩ K)
        @test isempty((2..5) ∩ (7..10))
        @test isempty((1..10) ∩ (7..2))
        A = Float16(1.1)..Float16(1.234)
        B = Float16(1.235)..Float16(1.3)
        C = Float16(1.236)..Float16(1.3)
        D = Float16(1.1)..Float16(1.236)
        @test D ∪ B == Float16(1.1)..Float16(1.3)
        @test D ∪ C == Float16(1.1)..Float16(1.3)
        @test_throws(ArgumentError, A ∪ C)

        @test 1.5 ∉ 0..1
        @test 1.5 ∉ 2..3
        @test_broken (0..1) ∪ (2..3) ≠ 0..3
        # even though A and B contain all Float32s between their extrema,
        # union should not return an interval as there exists a Float64
        # inbetween
        x32 = nextfloat(rightendpoint(A))
        x64 = nextfloat(Float64(rightendpoint(A)))
        @test x32 ∉ A
        @test x32 ∈ B
        @test x64 ∉ A
        @test x64 ∉ B
        # these tests
        @test_broken x32 ∈ A ∪ B
        @test_broken x64 ∉ A ∪ B

        @test J ⊆ L
        @test (L ⊆ J) == false
        @test K ⊆ I
        @test ClosedInterval(1, 2) ⊆ I
        @test I ⊆ I
        @test (ClosedInterval(7, 9) ⊆ I) == false
        @test I ⊇ I
        @test I ⊇ ClosedInterval(1, 2)

        @test hash(1..3) == hash(1.0..3.0)

        @test width(I) == 3
        @test width(J) == 0

        @test width(ClosedInterval(3,7)) ≡ 4
        @test width(ClosedInterval(4.0,8.0)) ≡ 4.0

        @test mean(0..1) == 0.5

        @test promote(1..2, 1.0..2.0) === (1.0..2.0, 1.0..2.0)

        @test duration(1..2) == 2
        # duration deliberately not defined for non-integer intervals
        @test_throws MethodError duration(1.2..2.4)
    end

    @testset "Day interval" begin
        A = Date(1990, 1, 1); B = Date(1990, 3, 1)
        @test width(ClosedInterval(A, B)) == Dates.Day(59)
        @test width(ClosedInterval(B, A)) == Dates.Day(0)
        @test isempty(ClosedInterval(B, A))
        @test duration(ClosedInterval(A, B)) ≡ 60
        @test duration(ClosedInterval(B, A)) ≡ 0
    end

    @testset "Convert" begin
        I = 0..3
        @test @inferred(convert(ClosedInterval{Float64}, I))         ===
                @inferred(convert(AbstractInterval{Float64}, I))     ===
                @inferred(convert(Domain{Float64}, I))  ===
                @inferred(ClosedInterval{Float64}(I))                ===
                @inferred(convert(TypedEndpointsInterval{:closed,:closed,Float64},I)) ===
                0.0..3.0
        @test @inferred(convert(ClosedInterval, I))                  ===
                @inferred(convert(Interval, I))                      ===
                @inferred(ClosedInterval(I))                         ===
                @inferred(Interval(I))                               ===
                @inferred(convert(AbstractInterval, I))              ===
                @inferred(convert(Domain, I))           ===
                @inferred(convert(TypedEndpointsInterval{:closed,:closed}, I)) ===
                @inferred(convert(TypedEndpointsInterval{:closed,:closed,Int}, I)) ===
                @inferred(convert(ClosedInterval{Int}, I)) === I
        @test_throws InexactError convert(OpenInterval, I)
        @test_throws InexactError convert(Interval{:open,:closed}, I)
        @test_throws InexactError convert(Interval{:closed,:open}, I)
        @test !(convert(ClosedInterval{Float64}, I) === 0..3)
        @test ClosedInterval{Float64}(1,3) === 1.0..3.0
        @test ClosedInterval(0.5..2.5) === 0.5..2.5
        @test ClosedInterval{Int}(1.0..3.0) === 1..3
        J = OpenInterval(I)
        @test_throws InexactError convert(ClosedInterval, J)
        @test @inferred(convert(OpenInterval{Float64}, J))         ===
                @inferred(convert(AbstractInterval{Float64}, J))     ===
                @inferred(convert(Domain{Float64}, J)) ===
                @inferred(OpenInterval{Float64}(J))                === OpenInterval(0.0..3.0)
        @test @inferred(convert(OpenInterval, J))                ===
                @inferred(convert(Interval, J))                      ===
                @inferred(convert(AbstractInterval, J))              ===
                @inferred(convert(Domain, J))           ===
                @inferred(OpenInterval(J))                          ===
                @inferred(OpenInterval{Int}(J)) ===
                @inferred(convert(OpenInterval{Int},J)) === OpenInterval(J)
        J = Interval{:open,:closed}(I)
        @test_throws InexactError convert(Interval{:closed,:open}, J)
        @test @inferred(convert(Interval{:open,:closed,Float64}, J))         ===
                @inferred(convert(AbstractInterval{Float64}, J))     ===
                @inferred(convert(Domain{Float64}, J)) ===
                @inferred(Interval{:open,:closed,Float64}(J))                === Interval{:open,:closed}(0.0..3.0)
        @test @inferred(convert(Interval{:open,:closed}, J))                ===
                @inferred(convert(Interval, J))                      ===
                @inferred(convert(AbstractInterval, J))              ===
                @inferred(convert(Domain, J))           ===
                @inferred(Interval{:open,:closed}(J))                          === Interval{:open,:closed}(J)
        J = Interval{:closed,:open}(I)
        @test_throws InexactError convert(Interval{:open,:closed}, J)
        @test @inferred(convert(Interval{:closed,:open,Float64}, J))         ===
                @inferred(convert(AbstractInterval{Float64}, J))     ===
                @inferred(convert(Domain{Float64}, J)) ===
                @inferred(Interval{:closed,:open,Float64}(J))                === Interval{:closed,:open}(0.0..3.0)
        @test @inferred(convert(Interval{:closed,:open}, J))                ===
                @inferred(convert(Interval, J))                      ===
                @inferred(convert(AbstractInterval, J))              ===
                @inferred(convert(Domain, J))           ===
                @inferred(Interval{:closed,:open}(J))                          === Interval{:closed,:open}(J)

        @test 1.0..2.0 === 1.0..2 === 1..2.0 === ClosedInterval{Float64}(1..2) ===
                Interval(1.0,2.0)

        @test convert(AbstractInterval, 1.0) ==
                convert(AbstractInterval{Float64}, 1) ==
                convert(TypedEndpointsInterval{:closed,:closed}, 1.0) ==
                convert(TypedEndpointsInterval{:closed,:closed,Float64}, 1) ==
                convert(ClosedInterval, 1.0) ==
                convert(ClosedInterval{Float64}, 1) ==
                1.0..1.0

        @test promote_type(Interval{:closed,:open,Float64}, Interval{:closed,:open,Int}) ===
                        Interval{:closed,:open,Float64}
    end

    @testset "Interval tests" begin
        for T in (Float32,Float64,BigFloat)
            d = zero(T) .. one(T)
            @test T(0.5) ∈ d
            @test T(1.1) ∉ d
            @test 0.5f0 ∈ d
            @test 1.1f0 ∉ d
            @test BigFloat(0.5) ∈ d
            @test BigFloat(1.1) ∉ d
            @test leftendpoint(d) ∈ d
            @test BigFloat(leftendpoint(d)) ∈ d
            @test nextfloat(leftendpoint(d)) ∈ d
            @test nextfloat(BigFloat(leftendpoint(d))) ∈ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test rightendpoint(d) ∈ d
            @test BigFloat(rightendpoint(d)) ∈ d
            @test nextfloat(rightendpoint(d)) ∉ d
            @test nextfloat(BigFloat(rightendpoint(d))) ∉ d
            @test prevfloat(rightendpoint(d)) ∈ d
            @test prevfloat(rightendpoint(d)) ∈ d

            @test leftendpoint(d) == zero(T)
            @test rightendpoint(d) == one(T)
            @test minimum(d) == infimum(d) == leftendpoint(d)
            @test maximum(d) == supremum(d) == rightendpoint(d)

            @test IntervalSets.isclosedset(d)
            @test !IntervalSets.isopenset(d)
            @test IntervalSets.isleftclosed(d)
            @test !IntervalSets.isleftopen(d)
            @test !IntervalSets.isrightopen(d)
            @test IntervalSets.isrightclosed(d)

            @test convert(AbstractInterval, d) ≡ d
            @test convert(AbstractInterval{T}, d) ≡ d
            @test convert(IntervalSets.Domain, d) ≡ d
            @test convert(IntervalSets.Domain{T}, d) ≡ d

            d = OpenInterval(zero(T) .. one(T))
            @test IntervalSets.isopenset(d)
            @test !IntervalSets.isclosedset(d)
            @test IntervalSets.isopenset(d)
            @test !IntervalSets.isclosedset(d)
            @test !IntervalSets.isleftclosed(d)
            @test IntervalSets.isleftopen(d)
            @test IntervalSets.isrightopen(d)
            @test !IntervalSets.isrightclosed(d)
            @test leftendpoint(d) ∉ d
            @test BigFloat(leftendpoint(d)) ∉ d
            @test nextfloat(leftendpoint(d)) ∈ d
            @test nextfloat(BigFloat(leftendpoint(d))) ∈ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test rightendpoint(d) ∉ d
            @test BigFloat(rightendpoint(d)) ∉ d
            @test nextfloat(rightendpoint(d)) ∉ d
            @test nextfloat(BigFloat(rightendpoint(d))) ∉ d
            @test prevfloat(rightendpoint(d)) ∈ d
            @test prevfloat(rightendpoint(d)) ∈ d
            @test infimum(d) == leftendpoint(d)
            @test supremum(d) == rightendpoint(d)
            @test_throws ArgumentError minimum(d)
            @test_throws ArgumentError maximum(d)

            @test isempty(OpenInterval(1,1))

            d = Interval{:open,:closed}(zero(T) .. one(T))
            @test !IntervalSets.isopenset(d)
            @test !IntervalSets.isclosedset(d)
            @test !IntervalSets.isleftclosed(d)
            @test IntervalSets.isleftopen(d)
            @test !IntervalSets.isrightopen(d)
            @test IntervalSets.isrightclosed(d)
            @test leftendpoint(d) ∉ d
            @test BigFloat(leftendpoint(d)) ∉ d
            @test nextfloat(leftendpoint(d)) ∈ d
            @test nextfloat(BigFloat(leftendpoint(d))) ∈ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test prevfloat(BigFloat(leftendpoint(d))) ∉ d
            @test rightendpoint(d) ∈ d
            @test BigFloat(rightendpoint(d)) ∈ d
            @test nextfloat(rightendpoint(d)) ∉ d
            @test nextfloat(BigFloat(rightendpoint(d))) ∉ d
            @test prevfloat(rightendpoint(d)) ∈ d
            @test prevfloat(BigFloat(rightendpoint(d))) ∈ d
            @test infimum(d) == leftendpoint(d)
            @test maximum(d) == supremum(d) == rightendpoint(d)
            @test_throws ArgumentError minimum(d)

            d = Interval{:closed,:open}(zero(T) .. one(T))
            @test !IntervalSets.isopenset(d)
            @test !IntervalSets.isclosedset(d)
            @test IntervalSets.isleftclosed(d)
            @test !IntervalSets.isleftopen(d)
            @test IntervalSets.isrightopen(d)
            @test !IntervalSets.isrightclosed(d)
            @test leftendpoint(d) ∈ d
            @test BigFloat(leftendpoint(d)) ∈ d
            @test nextfloat(leftendpoint(d)) ∈ d
            @test nextfloat(BigFloat(leftendpoint(d))) ∈ d
            @test prevfloat(leftendpoint(d)) ∉ d
            @test prevfloat(BigFloat(leftendpoint(d))) ∉ d
            @test rightendpoint(d) ∉ d
            @test BigFloat(rightendpoint(d)) ∉ d
            @test nextfloat(rightendpoint(d)) ∉ d
            @test nextfloat(BigFloat(rightendpoint(d))) ∉ d
            @test prevfloat(rightendpoint(d)) ∈ d
            @test prevfloat(BigFloat(rightendpoint(d))) ∈ d
            @test infimum(d) == minimum(d) == leftendpoint(d)
            @test supremum(d) == rightendpoint(d)
            @test_throws ArgumentError maximum(d)

            # - empty interval
            @test isempty(one(T) .. zero(T))
            @test zero(T) ∉ one(T) .. zero(T)

            d = one(T) .. zero(T)
            @test_throws ArgumentError minimum(d)
            @test_throws ArgumentError maximum(d)
            @test_throws ArgumentError infimum(d)
            @test_throws ArgumentError supremum(d)
        end
    end

    @testset "Issubset" begin
        I = 0..3
        J = 1..2
        @test J ⊆ I
        @test I ⊈ J
        @test OpenInterval(J) ⊆ I
        @test OpenInterval(I) ⊈ J
        @test J ⊆ OpenInterval(I)
        @test I ⊈ OpenInterval(J)
        @test OpenInterval(J) ⊆ OpenInterval(I)
        @test OpenInterval(I) ⊈ OpenInterval(J)
        @test Interval{:closed,:open}(J) ⊆ OpenInterval(I)
        @test Interval{:open,:closed}(J) ⊆ OpenInterval(I)
        @test Interval{:open,:closed}(J) ⊆ Interval{:open,:closed}(I)
        @test OpenInterval(I) ⊈ OpenInterval(J)

        @test Interval{:closed,:open}(J) ⊆ I
        @test I ⊈ Interval{:closed,:open}(J)

        @test I ⊆ I
        @test OpenInterval(I) ⊆ I
        @test Interval{:open,:closed}(I) ⊆ I
        @test Interval{:closed,:open}(I) ⊆ I
        @test I ⊈ OpenInterval(I)
        @test I ⊈ Interval{:open,:closed}(I)
        @test I ⊈ Interval{:closed,:open}(I)

        @test Interval{:closed,:open}(I) ⊆ Interval{:closed,:open}(I)
        @test Interval{:open,:closed}(I) ⊈ Interval{:closed,:open}(I)

        @test !isequal(I, OpenInterval(I))
        @test !(I == OpenInterval(I))
    end

    @testset "Union and intersection" begin
        for T in (Float32,Float64)
            i1 = zero(T) .. one(T)
            i2 = one(T)/3 .. one(T)/2
            i3 = one(T)/2 .. 2*one(T)
            i4 = T(2) .. T(3)
            i5 = nextfloat(one(T)) .. 2one(T)
            i_empty = one(T) ..zero(T)

            # - union of completely overlapping intervals
            @test i1 ∪ i2 ≡ i2 ∪ i1 ≡ i1
            @test Interval{:open,:closed}(i1) ∪ Interval{:open,:closed}(i2) ≡
                  Interval{:open,:closed}(i2) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)
            @test Interval{:closed,:open}(i1) ∪ Interval{:closed,:open}(i2) ≡
                  Interval{:closed,:open}(i2) ∪ Interval{:closed,:open}(i1) ≡ Interval{:closed,:open}(i1)
            @test OpenInterval(i1) ∪ OpenInterval(i2) ≡
                    OpenInterval(i2) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test i1 ∪ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(i2) ∪ i1 ≡ i1
            @test i1 ∪ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∪ i1 ≡ i1
            @test i1 ∪ OpenInterval(i2) ≡ OpenInterval(i2) ∪ i1 ≡ i1
            @test OpenInterval(i1) ∪ i2 ≡ i2 ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test OpenInterval(i1) ∪ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(i2) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test OpenInterval(i1) ∪ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test Interval{:open,:closed}(i1) ∪ OpenInterval(i2) ≡ OpenInterval(i2) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)
            @test Interval{:open,:closed}(i1) ∪ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)

            # - intersection of completely overlapping intervals
            @test i1 ∩ i2 ≡ i2 ∩ i1 ≡ i2
            @test Interval{:open,:closed}(i1) ∩ Interval{:open,:closed}(i2) ≡
                  Interval{:open,:closed}(i2) ∩ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i2)
            @test Interval{:closed,:open}(i1) ∩ Interval{:closed,:open}(i2) ≡
                  Interval{:closed,:open}(i2) ∩ Interval{:closed,:open}(i1) ≡ Interval{:closed,:open}(i2)
            @test OpenInterval(i1) ∩ OpenInterval(i2) ≡
                    OpenInterval(i2) ∩ OpenInterval(i1) ≡ OpenInterval(i2)
            @test i1 ∩ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(i2) ∩ i1 ≡  Interval{:open,:closed}(i2)
            @test i1 ∩ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∩ i1 ≡ Interval{:closed,:open}(i2)
            @test i1 ∩ OpenInterval(i2) ≡ OpenInterval(i2) ∩ i1 ≡ OpenInterval(i2)
            @test OpenInterval(i1) ∩ i2 ≡ i2 ∩ OpenInterval(i1) ≡ i2
            @test OpenInterval(i1) ∩ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(i2) ∩ OpenInterval(i1) ≡ Interval{:open,:closed}(i2)
            @test OpenInterval(i1) ∩ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∩ OpenInterval(i1) ≡ Interval{:closed,:open}(i2)
            @test Interval{:open,:closed}(i1) ∩ OpenInterval(i2) ≡ OpenInterval(i2) ∩ Interval{:open,:closed}(i1) ≡ OpenInterval(i2)
            @test Interval{:open,:closed}(i1) ∩ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(i2) ∩ Interval{:open,:closed}(i1) ≡ Interval{:closed,:open}(i2)

            # - union of partially overlapping intervals
            d = zero(T) .. 2*one(T)
            @test i1 ∪ i3 ≡ i3 ∪ i1 ≡ d
            @test Interval{:open,:closed}(i1) ∪ Interval{:open,:closed}(i3) ≡
                  Interval{:open,:closed}(i3) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(d)
            @test Interval{:closed,:open}(i1) ∪ Interval{:closed,:open}(i3) ≡
                  Interval{:closed,:open}(i3) ∪ Interval{:closed,:open}(i1) ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i1) ∪ OpenInterval(i3) ≡
                  OpenInterval(i3) ∪ OpenInterval(i1) ≡ OpenInterval(d)
            @test i1 ∪ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∪ i1 ≡ d
            @test i1 ∪ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∪ i1 ≡ Interval{:closed,:open}(d)
            @test i1 ∪ OpenInterval(i3) ≡ OpenInterval(i3) ∪ i1 ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i1) ∪ i3 ≡ i3 ∪ OpenInterval(i1) ≡ Interval{:open,:closed}(d)
            @test OpenInterval(i1) ∪ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∪ OpenInterval(i1) ≡ Interval{:open,:closed}(d)
            @test OpenInterval(i1) ∪ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∪ OpenInterval(i1) ≡ OpenInterval(d)
            @test Interval{:open,:closed}(i1) ∪ OpenInterval(i3) ≡ OpenInterval(i3) ∪ Interval{:open,:closed}(i1) ≡ OpenInterval(d)
            @test Interval{:open,:closed}(i1) ∪ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∪ Interval{:open,:closed}(i1) ≡ OpenInterval(d)

            # - intersection of partially overlapping intervals
            d = one(T)/2 .. one(T)
            @test i1 ∩ i3 ≡ i3 ∩ i1 ≡ d
            @test Interval{:open,:closed}(i1) ∩ Interval{:open,:closed}(i3) ≡
                  Interval{:open,:closed}(i3) ∩ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(d)
            @test Interval{:closed,:open}(i1) ∩ Interval{:closed,:open}(i3) ≡
                  Interval{:closed,:open}(i3) ∩ Interval{:closed,:open}(i1) ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i1) ∩ OpenInterval(i3) ≡
                    OpenInterval(i3) ∩ OpenInterval(i1) ≡ OpenInterval(d)
            @test i1 ∩ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∩ i1 ≡ Interval{:open,:closed}(d)
            @test i1 ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ i1 ≡ d
            @test i1 ∩ OpenInterval(i3) ≡ OpenInterval(i3) ∩ i1 ≡ Interval{:open,:closed}(d)
            @test OpenInterval(i1) ∩ i3 ≡ i3 ∩ OpenInterval(i1) ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i1) ∩ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∩ OpenInterval(i1) ≡ OpenInterval(d)
            @test OpenInterval(i1) ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ OpenInterval(i1) ≡ Interval{:closed,:open}(d)
            @test Interval{:open,:closed}(i1) ∩ OpenInterval(i3) ≡ OpenInterval(i3) ∩ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(d)
            @test Interval{:open,:closed}(i1) ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ Interval{:open,:closed}(i1) ≡ d

            # - union of barely overlapping intervals
            d = one(T)/3 .. 2*one(T)
            @test i2 ∪ i3 ≡ i3 ∪ i2 ≡ d
            @test Interval{:open,:closed}(i2) ∪ Interval{:open,:closed}(i3) ≡
                  Interval{:open,:closed}(i3) ∪ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(d)
            @test Interval{:closed,:open}(i2) ∪ Interval{:closed,:open}(i3) ≡
                  Interval{:closed,:open}(i3) ∪ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(d)
            @test_throws ArgumentError OpenInterval(i2) ∪ OpenInterval(i3)
            @test i2 ∪ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∪ i2 ≡ d
            @test i2 ∪ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∪ i2 ≡ Interval{:closed,:open}(d)
            @test i2 ∪ OpenInterval(i3) ≡ OpenInterval(i3) ∪ i2 ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i2) ∪ i3 ≡ i3 ∪ OpenInterval(i2) ≡ Interval{:open,:closed}(d)
            @test_throws ArgumentError OpenInterval(i2) ∪ Interval{:open,:closed}(i3)
            @test Interval{:open,:closed}(i2) ∪ OpenInterval(i3) ≡ OpenInterval(i3) ∪ Interval{:open,:closed}(i2) ≡ OpenInterval(d)
            @test Interval{:open,:closed}(i2) ∪ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∪ Interval{:open,:closed}(i2) ≡ OpenInterval(d)

            # - intersection of barely overlapping intervals
            d = one(T)/2 .. one(T)/2
            @test i2 ∩ i3 ≡ i3 ∩ i2 ≡ d
            @test Interval{:open,:closed}(i2) ∩ Interval{:open,:closed}(i3) ≡
                  Interval{:open,:closed}(i3) ∩ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(d)
            @test Interval{:closed,:open}(i2) ∩ Interval{:closed,:open}(i3) ≡
                  Interval{:closed,:open}(i3) ∩ Interval{:closed,:open}(i2) ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i2) ∩ OpenInterval(i3) ≡
                    OpenInterval(i3) ∩ OpenInterval(i2) ≡ OpenInterval(d)
            @test i2 ∩ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∩ i2 ≡ Interval{:open,:closed}(d)
            @test i2 ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ i2 ≡ d
            @test i2 ∩ OpenInterval(i3) ≡ OpenInterval(i3) ∩ i2 ≡ Interval{:open,:closed}(d)
            @test OpenInterval(i2) ∩ i3 ≡ i3 ∩ OpenInterval(i2) ≡ Interval{:closed,:open}(d)
            @test OpenInterval(i2) ∩ Interval{:open,:closed}(i3) ≡ Interval{:open,:closed}(i3) ∩ OpenInterval(i2) ≡ OpenInterval(d)
            @test OpenInterval(i2) ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ OpenInterval(i2) ≡ Interval{:closed,:open}(d)
            @test Interval{:open,:closed}(i2) ∩ OpenInterval(i3) ≡ OpenInterval(i3) ∩ Interval{:open,:closed}(i2) ≡ Interval{:open,:closed}(d)
            @test Interval{:open,:closed}(i2) ∩ Interval{:closed,:open}(i3) ≡ Interval{:closed,:open}(i3) ∩ Interval{:open,:closed}(i2) ≡ d

            # - intersection of custom intervals
            @test intersect(MyUnitInterval(true,true), MyUnitInterval(false,false)) == OpenInterval(0,1)
            @test intersect(MyUnitInterval(true,true), OpenInterval(0,1)) == OpenInterval(0,1)

            # - union of non-overlapping intervals
            @test_throws ArgumentError i1 ∪ i4
            @test_throws ArgumentError i4 ∪ i1
            @test_throws ArgumentError OpenInterval(i1) ∪ i4
            @test_throws ArgumentError i1 ∪ OpenInterval(i4)
            @test_throws ArgumentError Interval{:closed,:open}(i1) ∪ i4
            @test_throws ArgumentError Interval{:closed,:open}(i1) ∪ OpenInterval(i4)

            # - union of almost-overlapping intervals
            @test_throws ArgumentError i1 ∪ i5
            @test_throws ArgumentError i5 ∪ i1
            @test_throws ArgumentError OpenInterval(i1) ∪ i5
            @test_throws ArgumentError i1 ∪ OpenInterval(i5)
            @test_throws ArgumentError Interval{:closed,:open}(i1) ∪ i5
            @test_throws ArgumentError Interval{:closed,:open}(i1) ∪ OpenInterval(i5)

            # - intersection of non-overlapping intervals
            @test isempty(i1 ∩ i4)
            @test isempty(i4 ∩ i1)
            @test isempty(OpenInterval(i1) ∩ i4)
            @test isempty(i1 ∩ OpenInterval(i4))
            @test isempty(Interval{:closed,:open}(i1) ∩ i4)

            # - intersection of almost-overlapping intervals
            @test isempty(i1 ∩ i5)
            @test isempty(i5 ∩ i1)
            @test isempty(OpenInterval(i1) ∩ i5)
            @test isempty(i1 ∩ OpenInterval(i5))
            @test isempty(Interval{:closed,:open}(i1) ∩ i5)

            # - union of interval with empty
            @test i1 ∪ i_empty ≡ i_empty ∪ i1 ≡ i1
            @test Interval{:open,:closed}(i1) ∪ Interval{:open,:closed}(i_empty) ≡
                  Interval{:open,:closed}(i_empty) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)
            @test Interval{:closed,:open}(i1) ∪ Interval{:closed,:open}(i_empty) ≡
                  Interval{:closed,:open}(i_empty) ∪ Interval{:closed,:open}(i1) ≡ Interval{:closed,:open}(i1)
            @test OpenInterval(i1) ∪ OpenInterval(i_empty) ≡
                    OpenInterval(i_empty) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test i1 ∪ Interval{:open,:closed}(i_empty) ≡ Interval{:open,:closed}(i_empty) ∪ i1 ≡ i1
            @test i1 ∪ Interval{:closed,:open}(i_empty) ≡ Interval{:closed,:open}(i_empty) ∪ i1 ≡ i1
            @test i1 ∪ OpenInterval(i_empty) ≡ OpenInterval(i_empty) ∪ i1 ≡ i1
            @test OpenInterval(i1) ∪ i_empty ≡ i_empty ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test OpenInterval(i1) ∪ Interval{:open,:closed}(i_empty) ≡ Interval{:open,:closed}(i_empty) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test OpenInterval(i1) ∪ Interval{:closed,:open}(i_empty) ≡ Interval{:closed,:open}(i_empty) ∪ OpenInterval(i1) ≡ OpenInterval(i1)
            @test Interval{:open,:closed}(i1) ∪ OpenInterval(i_empty) ≡ OpenInterval(i_empty) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)
            @test Interval{:open,:closed}(i1) ∪ Interval{:closed,:open}(i_empty) ≡ Interval{:closed,:open}(i_empty) ∪ Interval{:open,:closed}(i1) ≡ Interval{:open,:closed}(i1)

            # - intersection of interval with empty
            @test isempty(i1 ∩ i_empty)
            @test isempty(i_empty ∩ i1)
            @test isempty(OpenInterval(i1) ∩ i_empty)
            @test isempty(i1 ∩ OpenInterval(i_empty))
            @test isempty(Interval{:closed,:open}(i1) ∩ i_empty)

            # - test matching endpoints
            @test (0..1) ∪ OpenInterval(0..1) ≡ OpenInterval(0..1) ∪ (0..1) ≡  0..1
            @test Interval{:open,:closed}(0..1) ∪ OpenInterval(0..1) ≡
                    OpenInterval(0..1) ∪ Interval{:open,:closed}(0..1) ≡
                    Interval{:open,:closed}(0..1)
            @test Interval{:closed,:open}(0..1) ∪ OpenInterval(0..1) ≡
                    OpenInterval(0..1) ∪ Interval{:closed,:open}(0..1) ≡
                    Interval{:closed,:open}(0..1)
        end
    end

    @testset "Empty" begin
        for T in (Float32,Float64)
            @test isempty(Interval{:open,:open}(zero(T),zero(T)))
            @test zero(T) ∉ Interval{:open,:open}(zero(T),zero(T))
            @test isempty(Interval{:open,:closed}(zero(T),zero(T)))
            @test zero(T) ∉ Interval{:open,:closed}(zero(T),zero(T))
            @test isempty(Interval{:closed,:open}(zero(T),zero(T)))
            @test zero(T) ∉ Interval{:closed,:open}(zero(T),zero(T))
        end
    end

    @testset "Custom intervals" begin
        I = MyUnitInterval(true,true)
        @test eltype(I) == eltype(typeof(I)) == Int
        @test leftendpoint(I) == 0
        @test rightendpoint(I) == 1
        @test isleftclosed(I)
        @test !isleftopen(I)
        @test isrightclosed(I)
        @test !isrightopen(I)
        @test ClosedInterval(I) === convert(ClosedInterval, I) ===
                ClosedInterval{Int}(I) === convert(ClosedInterval{Int}, I)  ===
                convert(Interval, I) === Interval(I) === 0..1
        @test_throws InexactError convert(OpenInterval, I)
        I = MyUnitInterval(false,false)
        @test leftendpoint(I) == 0
        @test rightendpoint(I) == 1
        @test !isleftclosed(I)
        @test !isrightclosed(I)
        @test OpenInterval(I) === convert(OpenInterval, I) ===
                OpenInterval{Int}(I) === convert(OpenInterval{Int}, I)  ===
                convert(Interval, I) === Interval(I) === OpenInterval(0..1)
        I = MyUnitInterval(false,true)
        @test leftendpoint(I) == 0
        @test rightendpoint(I) == 1
        @test isleftclosed(I) == false
        @test isrightclosed(I) == true
        @test Interval{:open,:closed}(I) === convert(Interval{:open,:closed}, I) ===
                Interval{:open,:closed,Int}(I) === convert(Interval{:open,:closed,Int}, I)  ===
                convert(Interval, I) === Interval(I) === Interval{:open,:closed}(0..1)
        I = MyUnitInterval(true,false)
        @test leftendpoint(I) == 0
        @test rightendpoint(I) == 1
        @test isleftclosed(I) == true
        @test isrightclosed(I) == false
        @test Interval{:closed,:open}(I) === convert(Interval{:closed,:open}, I) ===
                Interval{:closed,:open,Int}(I) === convert(Interval{:closed,:open,Int}, I)  ===
                convert(Interval, I) === Interval(I) === Interval{:closed,:open}(0..1)
        @test convert(AbstractInterval, I) === convert(AbstractInterval{Int}, I) === I
    end

    @testset "Custom typed endpoints interval" begin
        I = MyClosedUnitInterval()
        @test leftendpoint(I) == 0
        @test rightendpoint(I) == 1
        @test isleftclosed(I) == true
        @test isrightclosed(I) == true
        @test ClosedInterval(I) === convert(ClosedInterval, I) ===
                ClosedInterval{Int}(I) === convert(ClosedInterval{Int}, I)  ===
                convert(Interval, I) === Interval(I) === 0..1
        @test_throws InexactError convert(OpenInterval, I)
        @test I ∩ I === 0..1
        @test I ∩ (0.0..0.5) === 0.0..0.5
    end

    @testset "Missing endpoints" begin
        @test ismissing(2 in 1..missing)
        @test_broken ismissing(2 in missing..1)  # would be fixed by julialang#31171
    end

    @testset "in" begin
        @test in(0.1, 0.0..1.0) == true
        @test in(0.0, 0.0..1.0) == true
        @test in(1.1, 0.0..1.0) == false
        @test in(0.0, nextfloat(0.0)..1.0) == false
    end

    @testset "issubset" begin
        @test issubset(Interval{:closed,:closed}(1,2), Interval{:closed,:closed}(1,2)) == true
        @test issubset(Interval{:closed,:closed}(1,2), Interval{:open  ,:open  }(1,2)) == false
        @test issubset(Interval{:closed,:open  }(1,2), Interval{:open  ,:open  }(1,2)) == false
        @test issubset(Interval{:open  ,:closed}(1,2), Interval{:open  ,:open  }(1,2)) == false
        @test issubset(Interval{:closed,:closed}(1,2), Interval{:closed,:closed}(1,prevfloat(2.0))) == false
        @test issubset(Interval{:closed,:open  }(1,2), Interval{:open  ,:open  }(prevfloat(1.0),2)) == true
    end

    @testset "missing in" begin
        @test ismissing(missing in 0..1)
        @test !(missing in 1..0)
        @test ismissing(missing in OpenInterval(0, 1))
        @test ismissing(missing in Interval{:closed, :open}(0, 1))
        @test ismissing(missing in Interval{:open, :closed}(0, 1))
    end

    @testset "complex in" begin
        @test 0+im ∉ 0..2
        @test 0+0im ∈ 0..2
        @test 0+eps()im ∉ 0..2

        @test 0+im ∉ OpenInterval(0,2)
        @test 0+0im ∉ OpenInterval(0,2)
        @test 1+0im ∈ OpenInterval(0,2)
        @test 1+eps()im ∉ OpenInterval(0,2)

        @test 0+im ∉ Interval{:closed,:open}(0,2)
        @test 0+0im ∈ Interval{:closed,:open}(0,2)
        @test 1+0im ∈ Interval{:closed,:open}(0,2)
        @test 1+eps()im ∉ Interval{:closed,:open}(0,2)

        @test 0+im ∉ Interval{:open,:closed}(0,2)
        @test 0+0im ∉ Interval{:open,:closed}(0,2)
        @test 1+0im ∈ Interval{:open,:closed}(0,2)
        @test 1+eps()im ∉ Interval{:open,:closed}(0,2)
    end

    @testset "closedendpoints" begin
        @test closedendpoints(0..1) == closedendpoints(MyClosedUnitInterval()) == (true,true)
        @test closedendpoints(Interval{:open,:closed}(0,1)) == (false,true)
        @test closedendpoints(Interval{:closed,:open}(0,1)) == (true,false)
        @test closedendpoints(OpenInterval(0,1)) == (false,false)
    end

    @testset "OneTo" begin
        @test_throws ArgumentError Base.OneTo{Int}(0..5)
        @test_throws ArgumentError Base.OneTo(0..5)
        @test Base.OneTo(1..5) == Base.OneTo{Int}(1..5) == Base.OneTo(5)
        @test Base.Slice(1..5) == Base.Slice{UnitRange{Int}}(1..5) == Base.Slice(1:5)
    end

    @testset "IteratorSize" begin
        @test Base.IteratorSize(ClosedInterval) == Base.SizeUnknown()
    end

    @testset "IncompleteInterval" begin
        I = IncompleteInterval()
        @test eltype(I) === Int
        @test_throws ErrorException endpoints(I)
        @test_throws ErrorException closedendpoints(I)
    end

    include("findall.jl")
end

ext: jl


text:
module IntervalOptimization

using IntervalArithmetic
using ForwardDiff
using StaticArrays
using DataStructures

const gradient = ForwardDiff.gradient
const ∇ = gradient

include("basic_functions.jl")
include("branch_and_bound.jl")

end


text:
n=rand(1:10)
H = Hankel(collect(-n:n))
@test diag(Matrix(H)) == collect(-n:2:n)


text:
using LinearAlgebra
using Random

function RASJulia(X::Array{Float64,2},bw::Float64,lambda::Float64,c::Float64,t::Float64,epsilon::Float64,updating::Bool,nbFF::Int64)

N::Int64 = size(X,1);

#  STEP 1: RANDOM FOURIER FEATURES
B = Random_Fourier(X,bw,nbFF);
Id = Diagonal(ones(nbFF,nbFF));
b = cholesky(B'*B+N*lambda*Id);
F = (b.U)'\(B'); # P = F'*F
B = 0; b = 0; Id = 0;

lev_score::Array{Float64,2} = (sum(F.^2,dims=1));

## STEP 2: Run the sampling algorithm

idS = zeros(Int64,N,1);
weights = zeros(N,1);
nmbSamples::Int64 = 0;
epsilon::Float64 = 1e-10;

PS = Array{Float64,2}(undef,N,nbFF);
R = Array{Float64,2}(undef,1,1);

p_i::Float64 = 0;

for i = 1:N

    if mod(i,1000) == 1
    	println("Iteration: ",i)
        println("Number of landmarks sampled: ", nmbSamples);
    end

    if nmbSamples == 0
        p_i = min( 1, (1+t)*c*(1/epsilon)*lev_score[i]);

        if rand(1)[1] < p_i
            nmbSamples = 1;
            idS[nmbSamples] = i;
            weights[nmbSamples] = 1/sqrt(p_i);
            R[1,1] = (lev_score[i] .+ epsilon*weights[1].^(-2))^(-1);
            PS[:,1] = F'*F[:,idS[nmbSamples]];
        end
    else

        @views PS_i = PS[i,1:nmbSamples];
        temp = sum((PS_i'*R).*PS_i',dims=2);

        @inbounds p_i = min( 1, (1+t)*c*(1/epsilon)*max(0,lev_score[i] - temp[1]) );

        if rand(1)[1] < p_i

            nmbSamples = nmbSamples + 1;
            @inbounds idS[nmbSamples] = i;
            weights[nmbSamples] = 1/sqrt(p_i);

	    @views F_sliced = F[:,idS[nmbSamples]];
            PS[:,nmbSamples] = (F_sliced'*F)';#F'*F_sliced;

             if updating #Update inverse

                 @inbounds newA = PS[i,1:nmbSamples];
                 newA[end,end] = newA[end,end]+ epsilon*weights[nmbSamples].^(-2);
                 R = updateLS(R,newA);
             else
                @inbounds SPS = F[:,idS[1:nmbSamples]]'*F[:,idS[1:nmbSamples]];
                R = (SPS + diagm(0=>vec(epsilon*weights[1:nmbSamples].^(-2))))\diagm(0=>vec(ones(nmbSamples,1)));
             end
        end
    end
end

return idS, nmbSamples

end

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

function  gaussianRandomFeatures(D::Int64, p::Int64,gamma)
    # sample random Fourier directions and angles
    Omega = sqrt(2*gamma)*randn(p,D); # RVs defining RFF transform
    beta = rand(1,D)*2*pi;
    return Omega, beta
end

function transformFeatures( X::Array{Float64,2}, Omega::Array{Float64,2}, beta::Array{Float64,2} )

    D = size(Omega,2);
    Z = cos.(X*Omega .+ beta)*sqrt(2/D);

    return Z
end

function Random_Fourier(X::Array{Float64,2},bw::Float64,nmbFeatures::Int64)

n, p = size(X);
gamma = 1/(2*bw^2);
Omega, beta  = gaussianRandomFeatures(nmbFeatures,p,gamma);
Z = transformFeatures(X, Omega, beta ); # K = Z*Z';
return Z

end

function updateLS(K_inv::Array{Float64,2},newA::Array{Float64,1})

# Init
 a = newA[1:end-1];
 alpha = newA[end];

# Update K_inv
 B1 = K_inv .+ (K_inv*a*a'*K_inv)/(alpha - a'*K_inv*a);
 B2 = -K_inv*a/(alpha - a'*K_inv*a);
 B3 = -a'*K_inv/(alpha - a'*K_inv*a);
 B4 = 1/(alpha-a'*K_inv*a);
 K_inv = [B1 B2;B3 B4];

return K_inv

end


text:
module JuliennedArrays

import Base: axes, getindex, Int, setindex!, size
using Base: @inline, @propagate_inbounds, @pure, tail

export Slices, Align
export True, False

function is_in(::Needle, ::Needle, _...) where {Needle}
    True()
end
function is_in(needle, _, straws...)
    is_in(needle, straws...)
end
function is_in(_)
    False()
end

function in_unrolled(straws, needle1, needles...)
    is_in(needle1, straws...), in_unrolled(straws, needles...)...
end
function in_unrolled(_)
    ()
end

@pure function as_vals(them::Int...)
    map(Val, them)
end

abstract type TypedBool end
"""
    struct True

Typed `true`
"""
struct True <: TypedBool end
"""
    struct False

Typed `false`
"""
struct False <: TypedBool end

@inline function untyped(::True)
    true
end
@inline function untyped(::False)
    false
end

function not(::False)
    True()
end
function not(::True)
    False()
end

function getindex_unrolled(::Tuple{}, ::Tuple{})
    ()
end
function getindex_unrolled(into, switches)
    next = getindex_unrolled(tail(into), tail(switches))
    if untyped(first(switches))
        (first(into), next...)
    else
        next
    end
end

function setindex_unrolled(::Tuple{}, _, ::Tuple{})
    ()
end
function setindex_unrolled(old, new, switches)
    if untyped(first(switches))
        first(new), setindex_unrolled(tail(old), tail(new), tail(switches))...
    else
        first(old), setindex_unrolled(tail(old), new, tail(switches))...
    end
end

###
# Slices
###
struct Slices{Item,Dimensions,Whole,Alongs} <: AbstractArray{Item,Dimensions}
    whole::Whole
    alongs::Alongs
end
function Slices{Item,Dimensions}(whole::Whole, alongs::Alongs) where {Item,Dimensions,Whole,Alongs}
    Slices{Item,Dimensions,Whole,Alongs}(whole, alongs)
end

function axes(slices::Slices)
    getindex_unrolled(axes(slices.whole), map(not, slices.alongs))
end
function size(slices::Slices)
    map(length, axes(slices))
end

function slice_index(slices, indices)
    setindex_unrolled(axes(slices.whole), indices, map(not, slices.alongs))
end
@inline @propagate_inbounds function getindex(
    slices::Slices{Item,Dimensions},
    indices::Vararg{Int,Dimensions},
) where {Item,Dimensions}
    view(slices.whole, slice_index(slices, indices)...)
end
@inline @propagate_inbounds function setindex!(
    slices::Slices{Item,Dimensions},
    value,
    indices::Vararg{Int,Dimensions},
) where {Item,Dimensions}
    slices.whole[slice_index(slices, indices)...] = value
end

function axis_or_1(switch, axis)
    if untyped(switch)
        axis
    else
        1
    end
end

"""
    Slices(whole, alongs::TypedBool...)

Slice `whole` into `view`s.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be replaced with `:` when slicing.

```jldoctest slices_bools
julia> using JuliennedArrays

julia> whole = [1 2; 3 4];

julia> slices = Slices(whole, False(), True())
2-element Slices{SubArray{$Int, 1}, 1}:
 [1, 2]
 [3, 4]

julia> slices[1] == whole[1, :]
true

julia> slices[1] = [2, 1];

julia> whole
2×2 Matrix{$Int}:
 2  1
 3  4

julia> larger = rand(5, 5, 5);

julia> larger_slices = Slices(larger, True(), False(), False());

julia> size(first(larger_slices))
(5,)
```

You must include one along for each dimension of `whole`.

```jldoctest slices_bools
julia> using Test: @test_throws

julia> @test_throws MethodError Slices(whole, True());
```
"""
function Slices(
    whole::AbstractArray{<:Any,Dimensions},
    alongs::Vararg{TypedBool,Dimensions},
) where {Dimensions}
    Slices{
        typeof(@inbounds view(whole, map(axis_or_1, alongs, axes(whole))...)),
        length(getindex_unrolled(alongs, map(not, alongs))),
    }(
        whole,
        alongs,
    )
end

"""
    Slices(whole, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be replaced with `:` when slicing.

```jldoctest slices_ints
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{$Int}, 2, 2, 2) with eltype $Int:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> s = Slices(input, 1, 3)
2-element Slices{SubArray{$Int, 2}, 1}:
 [1 5; 2 6]
 [3 7; 4 8]

julia> map(sum, s)
2-element Vector{$Int}:
 14
 22
```

No along should be greater than the number of dimensions of `whole`.

```jldoctest slices_ints
julia> Slices(input, 4)
ERROR: ArgumentError: All alongs values (4,) should be less than or equal to 3
[...]
```

You can infer the result if the dimensions are constant.

```jldoctest slices_ints
julia> using Test: @inferred

julia> slices_1_3(x) = Slices(x, 1, 3);

julia> @inferred slices_1_3(input)
2-element Slices{SubArray{$Int, 2}, 1}:
 [1 5; 2 6]
 [3 7; 4 8]
```
"""
@inline function Slices(
    whole::AbstractArray{Item,Dimensions},
    alongs::Int...,
) where {Item,Dimensions}
    value_alongs = as_vals(alongs...)
    value_dimensions = ntuple(Val, Dimensions)
    if untyped(is_in(False(), in_unrolled(value_dimensions, value_alongs...)...))
        throw(
            ArgumentError(
                "All alongs values $(alongs) should be less than or equal to $(Dimensions)",
            ),
        )
    end
    Slices(whole, in_unrolled(value_alongs, value_dimensions...)...)
end

function Base.showarg(io::IO, ::Slices{Item,Dimensions}, toplevel) where {Item,Dimensions}
    print(
        io,
        "Slices{",
        basetype(Item),
        "{",
        eltype(Item),
        ", ",
        ndims(Item),
        "}, ",
        Dimensions,
        "}",
    )
end

# This is to be added to Julia (maybe under a different name)
# Follow https://github.com/JuliaLang/julia/issues/35543 for progress
function basetype(Item::Type)
    Base.typename(Item).wrapper
end
function basetype(Item)
    basetype(typeof(Item))
end

@inline function Int(::True)
    1
end
@inline function Int(::False)
    0
end

function match_dimensions(has_dimensions, used_dimensions::Int)
    if used_dimensions !== has_dimensions
        throw(
            DimensionMismatch(
                "$used_dimensions of $has_dimensions inner dimensions are used",
            ),
        )
    end
end

function check_dimensions(has_dimensions, ::Vararg{Int,Dimensions}) where {Dimensions}
    match_dimensions(has_dimensions, Dimensions)
end

function check_dimensions(has_dimensions, alongs::Vararg{TypedBool})
    match_dimensions(has_dimensions, mapreduce(Int, +, alongs))
end

function check_sizes(slices, alongs...)
    if !isempty(slices)
        first_slice = first(slices)
        first_axes = axes(first_slice)
        check_dimensions(ndims(first_slice), alongs...)
        for slice in slices
            if axes(slice) != first_axes
                throw(
                    ArgumentError(
                        "Axes of $slice does not match the axes of the first slice: $first_axes",
                    ),
                )
            end
        end
    end
end

###
# Align
###
struct Align{Item,Dimensions,Sliced,Alongs} <: AbstractArray{Item,Dimensions}
    slices::Sliced
    alongs::Alongs
end

function Align{Item,Dimensions}(
    slices::Sliced,
    alongs::NTuple{Dimensions,TypedBool},
) where {Item,Dimensions,Sliced}
    Align{Item,Dimensions,Sliced,typeof(alongs)}(slices, alongs)
end

function axes(aligned::Align)
    setindex_unrolled(
        setindex_unrolled(aligned.alongs, axes(aligned.slices), map(not, aligned.alongs)),
        axes(first(aligned.slices)),
        aligned.alongs,
    )
end
function size(aligned::Align)
    map(length, axes(aligned))
end

function split_indices(aligned, indices)
    getindex_unrolled(indices, map(not, aligned.alongs)),
    getindex_unrolled(indices, aligned.alongs)
end
@inline @propagate_inbounds function getindex(
    aligned::Align{Item,Dimensions},
    indices::Vararg{Int,Dimensions},
) where {Item,Dimensions}
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...]
end
@inline @propagate_inbounds function setindex!(
    aligned::Align{Item,Dimensions},
    value,
    indices::Vararg{Int,Dimensions},
) where {Item,Dimensions}
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...] = value
end

"""
    Align(slices, alongs::TypedBool...)

`Align` an array of arrays, all with the same size.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be taken up by the inner arrays. Inverse of [`Slices`](@ref).

```jldoctest align_bools
julia> using JuliennedArrays

julia> slices = [[1, 2], [3, 4]];

julia> aligned = Align(slices, False(), True())
2×2 Align{$Int, 2} with eltype $Int:
 1  2
 3  4

julia> aligned[1, :] == slices[1]
true

julia> aligned[1, 1] = 0;

julia> slices
2-element Vector{Vector{$Int}}:
 [0, 2]
 [3, 4]
```

Will throw an error if you try to align slices with different axes.
Use `@inbounds` to skip this check.

```jldoctest align_bools
julia> unequal = [[1], [1, 2]];

julia> Align(unequal, False(), True())
ERROR: ArgumentError: Axes of [1, 2] does not match the axes of the first slice: (Base.OneTo(1),)
[...]
```
"""
@inline @propagate_inbounds function Align(
    slices::AbstractArray{<:AbstractArray{Item,InnerDimensions},OuterDimensions},
    alongs::TypedBool...,
) where {Item,InnerDimensions,OuterDimensions}
    @boundscheck check_sizes(slices, alongs...)
    Align{Item,OuterDimensions + InnerDimensions}(slices, alongs)
end

"""
    Align(slices, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be taken up by the inner arrays.

```jldoctest align_ints
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{$Int}, 2, 2, 2) with eltype $Int:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> slices = Slices(input, 1, 3)
2-element Slices{SubArray{$Int, 2}, 1}:
 [1 5; 2 6]
 [3 7; 4 8]

julia> Align(slices, 1, 3)
2×2×2 Align{$Int, 3} with eltype $Int:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8
```

You must include one `along` for each inner dimension.

```jldoctest align_ints
julia> Align(slices, 1)
ERROR: DimensionMismatch("1 of 2 inner dimensions are used")
[...]
```

Julia can infer the result if `alongs` is constant.

```jldoctest align_ints
julia> using Test: @inferred

julia> align_1_3(x) = Align(x, 1, 3);

julia> @inferred align_1_3(slices)
2×2×2 Align{$Int, 3} with eltype $Int:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8
```
"""
@inline @propagate_inbounds function Align(
    slices::AbstractArray{<:AbstractArray{Item,InnerDimensions},OuterDimensions},
    alongs::Int...,
) where {Item,InnerDimensions,OuterDimensions}
    @boundscheck check_sizes(slices, alongs...)
    Align(
        slices,
        in_unrolled(
            as_vals(alongs...),
            ntuple(Val, InnerDimensions + OuterDimensions)...,
        )...,
    )
end

function Base.showarg(io::IO, ::Align{Item,Dimensions}, toplevel) where {Item,Dimensions}
    print(io, "Align{", Item, ", ", Dimensions, "}")
    toplevel && print(io, " with eltype ", Item)
end

end # module


text:
#=
   Urn puzzle in Turing.jl

   From Berlin Bayesians
   https://app.slack.com/client/TFPMSKW3F/CFQHMRD6K/thread/CFQHMRD6K-1623812230.000500
   """
   Three urns, first blindly take one from the first and put in the second,
   then blindly take one from the second and put in the third and then blindly pick
   one from the third. How likely is it to pick a black one in the last step?

   [
     Urns

     1: white black black
     2: white white black
     3: white black

   ]
   """

   Probability of drawing a black in the last step:
   Distributions of variable (num:0)
   0.00000 =>    5273  (0.527300)
   1.00000 =>    4727  (0.472700)

   (Cf the WebPPL model's exact solution: 0.47222222222222227)

   Cf ~/webppl/urn_puzzle
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model urn_puzzle() = begin
    white = 1
    black = 2

    urn1 = [white,black,black]
    urn2 = [white,white,black]
    urn3 = [white,black]

    # From urn 1 to urn 2
    # Note: We have to do this in two steps, i.e.
    # first pick some integer and then fetch that p1'th ball
    # from the urn
    p1 ~ DiscreteUniform(1,length(urn1))
    urn2b = push!(urn2,urn1[p1])

    # From urn 2 to urn 3
    p2 ~ DiscreteUniform(1,length(urn2b))
    urn3b = push!(urn3,urn2b[p2])

    # Pick a ball from urn 3
    p3 ~ DiscreteUniform(1,length(urn3b))
    ball = urn3b[p3]

    # Is it black?
    pblack ~ Dirac(ball==black)

end

model = urn_puzzle()

num_chns = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)
#

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:pblack)
show_var_dist_pct(chns,:ball)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
###########################################
# Define three modules with different
# but compatible data strucures
#############
module ModFoo

struct Foo
    foo::Vector
end

export Foo
end

#############
module ModBar

struct Bar
    bar::Vector
end

export Bar
end

#############
module ModBaz

struct Baz
    baz::Vector
end

export Baz
end

#############################################
# The Rosetta module stores code which knows
# how to convert between the data structures of the
# different modules without drawing in all the depedencies

# If this would be a package, dependecies only would be used
# in the unit tests, not while using this package

module Rosetta

# Dummy module to define default
module NoModule end

# Heuristic methods to determine
# a particular module
isModFoo(mod)=isdefined(mod,:Foo)
isModBar(mod)=isdefined(mod,:Bar)

# Concrete conversion methods
convert_foo_to_bar(foo,FooModule,BarModule)=BarModule.Bar(foo.foo)
convert_bar_to_foo(bar,BarModule,FooModule)=FooModule.Foo(bar.bar)

# Generic conversion method - the one "exported" from Rosetta
function Base.convert(newtype, x; From=NoModule, To=NoModule)
    if isModFoo(From)&&isModBar(To)
        convert_foo_to_bar(x,From,To)
    elseif isModBar(From)&&isModFoo(To)
        convert_bar_to_foo(x,From,To)
    else
        throw( ArgumentError("undefined conversion from $(typeof(x)) to $(newtype)") )
    end
end

end

#################################
# Sample code (using locally defined modules here)
using .ModFoo
using .ModBar
using .ModBaz

using .Rosetta

function testrosetta()
    foo=Foo(rand(100))
    bar=convert(Bar, foo,From=ModFoo, To=ModBar)
    foo1=convert(Foo, bar, From=ModBar, To=ModFoo)
    foo1===foo || println("transitivity error")
    try
        baz=convert(Baz,foo)
    catch e
        println("undefined conversion:")
        println(e)
    end
end


text:
f(x, y) = 2*x + 3*y
∂ₓf(x, y) = ForwardDiff.derivative(x->f(x, y), x)


text:
function init!(
    crr::CorefResolver,
    params::InitParams,
    )

    df_pair = crr.crs.esp_dict[params.entity_type]
    space_pair = crr |> main
    for (samp_n, refspace, df) in
        zip(
            [params.left_sample, params.right_sample],
            [space_pair.erspace1, space_pair.erspace2],
            [df_pair.df1, df_pair.df2]
            )
        init!(refspace, df, params.entity_type, samp_n)
    end
end

function init!(
    refspace::EntityRefSpace,
    df::DataFrame,
    entity_type::String,
    sample_size::Int,
    )
    refspace.width = 1
    refspace.entity_types[1] = entity_type
    dfh = size(df)[1]
    spacemaxh = size(refspace.space)[1]
    # TODO: warn about this next bit
    # note it in progress
    (sample_size > spacemaxh) && (sample_size = spacemaxh)
    if sample_size > dfh
        refspace.space[1:dfh, 1] = 1:dfh
        refspace.height = dfh
    else
        refspace.space[1:sample_size, 1] = StatsBase.sample(1:dfh, sample_size, replace=false)
        refspace.height = sample_size
    end
end


text:
export
    lastnb,
    latcyl,
    latrec,
    latsph,
    latsrf,
    lcase,
    ldpool,
    lgrind,
    limbpt,
    lmpool,
    lparse,
    lparsm,
    lparss,
    lspcn,
    lstlec,
    lstled,
    lstlei,
    lstle,
    lstltc,
    lstltd,
    lstlti,
    lstlt,
    ltime,
    lx4dec,
    lx4num,
    lx4sgn,
    lx4uns,
    lxqstr

function _lastnb(str)
    r = ccall((:lastnb_c, libcspice), SpiceInt, (Cstring,), str)
    handleerror()
    r + 1
end

@deprecate lastnb(str) findprev(!isspace, str, length(str))

"""
    lastnb(str)

!!! warning "Deprecated"
    Use `findprev(!isspace, str, length(str))` instead.
"""
lastnb

"""
    latcyl(radius, lon, lat)

Convert from latitudinal coordinates to cylindrical coordinates.

### Arguments ###

- `radius`: Distance of a point from the origin
- `lon`: Angle of the point from the XZ plane in radians
- `lat`: Angle of the point from the XY plane in radians

### Output ###

Return the tuple `(r, lonc, z)`.

- `r`: Distance of the point from the z axis
- `lonc`: Angle of the point from the XZ plane in radians. 'lonc' is set equal to 'lon'
- `z`: Height of the point above the XY plane

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/latcyl_c.html)
"""
function latcyl(radius, lon, lat)
    r = Ref{SpiceDouble}()
    lonc = Ref{SpiceDouble}()
    z = Ref{SpiceDouble}()
    ccall((:latcyl_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          radius, lon, lat, r, lonc, z)
    r[], lonc[], z[]
end

"""
    latrec(radius, lon, lat)

Convert from latitudinal coordinates to rectangular coordinates.

### Arguments ###

- `radius`: Distance of a point from the origin
- `lon`: Angle of the point from the XZ plane in radians
- `lat`: Angle of the point from the XY plane in radians

### Output ###

Return the rectangular coordinates vector of the point.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/latrec_c.html)
"""
function latrec(radius, lon, lat)
    rectan = Array{SpiceDouble}(undef, 3)
    ccall((:latrec_c, libcspice), Cvoid, (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          radius, lon, lat, rectan)
    rectan
end

"""
    latsph(radius, lon, lat)

Convert from latitudinal coordinates to rectangular coordinates.

### Arguments ###

- `radius`: Distance of a point from the origin
- `lon`: Angle of the point from the XZ plane in radians
- `lat`: Angle of the point from the XY plane in radians

### Output ###

Return the tuple `(rho, colat, lons)`.

- `rho`: Distance of the point from the origin
- `colat`: Angle of the point from positive z axis (radians)
- `lons`: Angle of the point from the XZ plane (radians)

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/latsph_c.html)
"""
function latsph(radius, lon, lat)
    rho = Ref{SpiceDouble}()
    colat = Ref{SpiceDouble}()
    lons = Ref{SpiceDouble}()
    ccall((:latsph_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble},  Ref{SpiceDouble},  Ref{SpiceDouble}),
          radius, lon, lat, rho, colat, lons)
    rho[], colat[], lons[]
end

"""
    latsrf(method, target, et, fixref, npts, lonlat)

Map array of planetocentric longitude/latitude coordinate pairs
to surface points on a specified target body.

The surface of the target body may be represented by a triaxial
ellipsoid or by topographic data provided by DSK files.

### Arguments ###

- `method`: Computation method
- `target`: Name of target body
- `et`: Epoch in TDB seconds past J2000 TDB
- `fixref`: Body-fixed, body-centered target body frame
- `lonlat`: Array of longitude/latitude coordinate pairs

### Output ###

Returns an array of surface points.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/latsrf_c.html)
"""
function latsrf(method, target, et, fixref, lonlat)
    npts = length(lonlat)
    lonlat_ = array_to_cmatrix(lonlat, n=2)
    srfpts = Array{SpiceDouble}(undef, 3, npts)
    ccall((:latsrf_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceDouble, Cstring, SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}),
          method, target, et, fixref, npts, lonlat_, srfpts)
    handleerror()
    cmatrix_to_array(srfpts)
end

_lcase(in) = _lcase(in,length(in)+1)

function _lcase(in,lenout)
    out = Array{SpiceChar}(undef, lenout)
    ccall((:lcase_c, libcspice), Cvoid, (Cstring, SpiceInt, Ref{SpiceChar}), in, lenout, out)
    handleerror()
    chararray_to_string(out)
end
"""
    lcase(in)

!!! warning "Deprecated"
    Use `lowercase(in)` instead.
"""
lcase

@deprecate lcase(in) lowercase(in)

"""
   ldpool(kernel)

Load the variables contained in a NAIF ASCII kernel file into the
kernel pool.

### Arguments ###

- `kernel`: Name of the kernel file

### Output ###

None

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ldpool_c.html)
"""
function ldpool(kernel)
    ccall((:ldpool_c, libcspice), Cvoid, (Cstring,), kernel)
    handleerror()
end

"""
    lgrind(xvals, yvals, x)

Evaluate a Lagrange interpolating polynomial for a specified
set of coordinate pairs, at a specified abscissa value.
Return the value of both polynomial and derivative.

### Arguments ###

- `xvals`: Abscissa values of coordinate pairs
- `yvals`: Ordinate values of coordinate pairs
- `x`: Point at which to interpolate the polynomial

### Output ###

- `p`: The value at x of the unique polynomial of
       degree n-1 that fits the points in the plane
       defined by xvals and yvals
- `dp`: The derivative at x of the interpolating
        polynomial described above

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lgrind_c.html)
"""
function lgrind(xvals, yvals, x)
    n = length(xvals)
    @checkdims n yvals
    p = Ref{SpiceDouble}()
    dp = Ref{SpiceDouble}()
    work = Matrix{SpiceDouble}(undef, 2, n)
    ccall((:lgrind_c, libcspice), Cvoid,
          (SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble},
           Ref{SpiceDouble}),
          n, xvals, yvals, work, x, p, dp)
    handleerror()
    p[], dp[]
end

"""
    limbpt(method, target, et, fixref, abcorr, corloc, obsrvr, refvec, rolstp, ncuts, schstp,
           soltol, maxn)

Find limb points on a target body. The limb is the set of points of tangency on the target of rays
emanating from the observer.  The caller specifies half-planes bounded by the observer-target
center vector in which to search for limb points.

The surface of the target body may be represented either by a triaxial ellipsoid or by
topographic data.

### Arguments ###

- `method`: Computation method
- `target`: Name of target body
- `et`: Epoch in ephemeris seconds past J2000 TDB
- `fixref`: Body-fixed, body-centered target body frame
- `abcorr`: Aberration correction
- `corloc`: Aberration correction locus
- `obsrvr`: Name of observing body
- `refvec`: Reference vector for cutting half-planes
- `rolstp`: Roll angular step for cutting half-planes
- `ncuts`: Number of cutting half-planes
- `schstp`: Angular step size for searching
- `soltol`: Solution convergence tolerance
- `maxn`: Maximum number of entries in output arrays

### Output ###

Returns the tuple `(npts, points, epochs, tangts)`.

- `npts`: Counts of limb points corresponding to cuts
- `points`: Limb points
- `epochs`: Times associated with limb points
- `tangts`: Tangent vectors emanating from the observer

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lgrind_c.html)
"""
function limbpt(method, target, et, fixref, abcorr, corloc, obsrvr, refvec,
                rolstp, ncuts, schstp, soltol, maxn)
    @checkdims 3 refvec
    npts = Array{SpiceInt}(undef, ncuts)
    points = Array{SpiceDouble}(undef, 3, maxn)
    epochs = Array{SpiceDouble}(undef, maxn)
    tangts = Array{SpiceDouble}(undef, 3, maxn)
    ccall((:limbpt_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceDouble, Cstring, Cstring, Cstring, Cstring,
          Ref{SpiceDouble}, SpiceDouble, SpiceInt, SpiceDouble, SpiceDouble,
          SpiceInt, Ref{SpiceInt}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          method, target, et, fixref, abcorr, corloc, obsrvr, refvec, rolstp,
          ncuts, schstp, soltol, maxn, npts, points, epochs, tangts)
    handleerror()
    npts, cmatrix_to_array(points), epochs, cmatrix_to_array(tangts)
end

@deprecate limb_pl02 limbpt
@deprecate llgrid_pl02 latsrf

"""
   lmpool(cvals)

Load the variables contained in an internal buffer into the
kernel pool.

### Arguments ###

- `cvals`: An array that contains a SPICE text kernel

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lmpool_c.html)
"""
function lmpool(cvals)
    cvals_, n, lenvals = chararray(cvals)
    ccall((:lmpool_c, libcspice), Cvoid, (Ref{SpiceChar}, SpiceInt, SpiceInt), cvals_, lenvals, n)
    handleerror()
end

function _lparse(list, delim, nmax)
    lenout = length(list)+1
    n = Ref{SpiceInt}()
    items = Array{SpiceChar}(undef, lenout, nmax)
    ccall((:lparse_c, libcspice), Cvoid, (Cstring, Cstring, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar}),
          list, delim, nmax, lenout, n , items)
    handleerror()
    chararray_to_string(items, n[])
end

@deprecate lparse(list, delim, nmax) split(list, delim, limit=nmax)

"""
   lparse(list, delim, nmax)

!!! warning "Deprecated"
    Use `split(list, delim, limit=nmax)` instead.
"""
lparse

function _lparsm(list, delims, nmax)
    lenout = length(list) + 1
    n = Ref{SpiceInt}()
    items = Array{SpiceChar}(undef, lenout, nmax)
    ccall((:lparsm_c, libcspice), Cvoid, (Cstring, Cstring, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar}),
          list, delims, nmax, lenout, n , items)
    handleerror()
    chararray_to_string(items, n[])
end

@deprecate lparsm(list, delim, nmax) split(list, delim, limit=nmax, keepempty=false)

"""
   lparsm(list, delims, nmax)

!!! warning "Deprecated"
    Use `split(list, delim, limit=nmax, keepempty=false)` instead.
"""
lparsm

function _lparss(list, delims)
    items = SpiceCharCell(length(list), length(list))
    ccall((:lparss_c, libcspice), Cvoid, (Cstring, Cstring, Ref{Cell{SpiceChar}}),
          list, delims, Ref(items.cell))
    handleerror()
    items
end

@deprecate lparss(list, delim) Set(split(list, collect(delim)))

"""
   lparss(list, delims)

!!! warning "Deprecated"
    Use `Set(split(list, collect(delim)))` instead.
"""
lparss

"""
   lspcn(body, et, abcorr)

Compute L_s, the planetocentric longitude of the sun, as seen
from a specified body.

### Arguments ###

- `body`: Name of the central body
- `et`: Epoch in seconds past J2000 TDB
- `abcorr`: Aberration correction

### Output ###

Returns the planetocentric longitude of the sun for the specified body
at the specified time in radians.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lspcn_c.html)
"""
function lspcn(body, et, abcorr)
    out = ccall((:lspcn_c, libcspice), SpiceDouble, (Cstring, SpiceDouble, Cstring),
         body, et, abcorr)
    handleerror()
    out
end

function _lstle(string::AbstractString, array)
    n = length(array)
    array_, n, lenvals = chararray(array)
    out = ccall((:lstlec_c, libcspice), SpiceInt, (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
                string, n, lenvals, array_)
    handleerror()
    out + 1
end

function _lstle(x::AbstractFloat, array)
    n = length(array)
    out = ccall((:lstled_c, libcspice), SpiceInt, (SpiceDouble, SpiceInt, Ref{SpiceDouble}),
                x, n, array)
    out + 1
end

function _lstle(x::Signed, array)
    n = length(array)
    array_ = SpiceInt.(array)
    out = ccall((:lstlei_c, libcspice), SpiceInt, (SpiceInt, SpiceInt, Ref{SpiceInt}),
                x, n, array_)
    out + 1
end

@deprecate lstlec(item, array) findfirst(item .<= array)
@deprecate lstled(item, array) searchsortedlast(array, item)
@deprecate lstlei(item, array) searchsortedlast(array, item)

"""
    lstlecd(x, array)

!!! warning "Deprecated"
    Use `findfirst(item .<= array)` instead.
"""
lstlec

"""
    lstle[di](x, array)

!!! warning "Deprecated"
    Use `searchsortedlast(array, x)` instead.
"""
lstled
lstlei

function _lstlt(string::AbstractString, array)
    n = length(array)
    array_, n, lenvals = chararray(array)
    out = ccall((:lstltc_c, libcspice), SpiceInt, (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
                string, n, lenvals, array_)
    handleerror()
    out + 1
end

function _lstlt(x::AbstractFloat, array)
    n = length(array)
    out = ccall((:lstltd_c, libcspice), SpiceInt, (SpiceDouble, SpiceInt, Ref{SpiceDouble}),
                x, n, array)
    out + 1
end

function _lstlt(x::Signed, array)
    n = length(array)
    array_ = SpiceInt.(array)
    out = ccall((:lstlti_c, libcspice), SpiceInt, (SpiceInt, SpiceInt, Ref{SpiceInt}),
                x, n, array_)
    out + 1
end

@deprecate lstltc(item, array) findfirst(item .< array)
@deprecate lstltd(item, array) searchsortedlast(array, item, lt=<=)
@deprecate lstlti(item, array) searchsortedlast(array, item, lt=<=)

"""
    lstltcd(x, array)

!!! warning "Deprecated"
    Use `findfirst(item .< array)` instead.
"""
lstltc

"""
    lstlt[di](x, array)

!!! warning "Deprecated"
    Use `searchsortedlast(array, x, lt=<=)` instead.
"""
lstltd
lstlti

"""
    ltime(etobs, obs, dir, targ)

This routine computes the transmit (or receive) time
of a signal at a specified target, given the receive
(or transmit) time at a specified observer. The elapsed
time between transmit and receive is also returned.

### Arguments ###

- `etobs`: Epoch of a signal at some observer
- `obs`: NAIF ID of some observer
- `dir`: Direction the signal travels ( "->" or "<-" )
- `targ`: Time between transmit and receipt of the signal

### Output ###

- `ettarg`: Epoch of the signal at the target
- `obs`: NAIF ID of some observer

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ltime_c.html)
"""
function ltime(etobs, obs, dir, targ)
    ettarg = Ref{SpiceDouble}()
    elapsd = Ref{SpiceDouble}()
    ccall((:ltime_c, libcspice), Cvoid,
          (SpiceDouble, SpiceInt, Cstring, SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}),
          etobs, obs, dir, targ, ettarg, elapsd)
    ettarg[], elapsd[]
end

"""
    lx4dec(string, first)

Scan a string from a specified starting position for the
end of a decimal number.

### Arguments ###

- `string`: Any character string
- `first`: First character to scan from in string

### Output ###

- `last`: Last character that is part of a decimal number. If there is no such
          character, last will be returned with the value first-1.
- `nchar`: Number of characters in the decimal number

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lx4dec_c.html)
"""
function lx4dec(string, first)
    last = Ref{SpiceInt}()
    nchar = Ref{SpiceInt}()
    ccall((:lx4dec_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
          string, first - 1, last, nchar)
    handleerror()
    last[] + 1, Int(nchar[])
end

"""
    lx4num(string, first)

Scan a string from a specified starting position for the
end of a number.

### Arguments ###

- `string`: Any character string
- `first`: First character to scan from in string

### Output ###

- `last`: Last character that is part of a number. If there is no such
          character, last will be returned with the value first-1.
- `nchar`: Number of characters in the number

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lx4num_c.html)
"""
function lx4num(string, first)
    last = Ref{SpiceInt}()
    nchar = Ref{SpiceInt}()
    ccall((:lx4num_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
          string, first - 1, last, nchar)
    handleerror()
    last[] + 1, Int(nchar[])
end

"""
    lx4sgn(string, first)

Scan a string from a specified starting position for the end of a signed integer.

### Arguments ###

- `string`: Any character string
- `first`: First character to scan from in string

### Output ###

- `last`: Last character that is part of a signed integer. If there is no such
          character, last will be returned with the value first-1.
- `nchar`: Number of characters in the signed integer

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lx4sgn_c.html)
"""
function lx4sgn(string, first)
    last = Ref{SpiceInt}()
    nchar = Ref{SpiceInt}()
    ccall((:lx4sgn_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
          string, first - 1, last, nchar)
    handleerror()
    last[] + 1, nchar[]
end

"""
    lx4uns(string, first)

Scan a string from a specified starting position for the
end of a unsigned integer.

### Arguments ###

- `string`: Any character string
- `first`: First character to scan from in string

### Output ###

- `last`: Last character that is part of an unsigned integer. If there is no such
          character, last will be returned with the value first-1.
- `nchar`: Number of characters in the unsigned integer

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lx4uns_c.html)
"""
function lx4uns(string, first)
    last = Ref{SpiceInt}()
    nchar = Ref{SpiceInt}()
    ccall((:lx4uns_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
          string, first - 1, last, nchar)
    handleerror()
    last[] + 1, nchar[]
end

"""
    lxqstr(string, qchar, first)

Lex (scan) a quoted string.

### Arguments ###

- `string`: String to be scanned
- `qchar`: Quote delimiter character
- `first`: Character position at which to start scanning

### Output ###

- `last`: Character position of end of token
- `nchar`: Number of characters in token

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/lxqstr_c.html)
"""
function lxqstr(string, qchar, first)
    last = Ref{SpiceInt}()
    nchar = Ref{SpiceInt}()
    ccall((:lxqstr_c, libcspice), Cvoid,
          (Cstring, SpiceChar, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
          string, qchar, first - 1, last, nchar)
    handleerror()
    last[] + 1, nchar[]
end


text:
@inline lassoonA!(args...) = lassoonA_glmnet!(args...)

################################################################## Helpers
function getβfromC(βs, C)
    p, K = size(βs)
    nprev = norm(βs[:, 1], 1) # = 0.
    for k in 2:K
        ncurr = norm(βs[:, k], 1)
        if C < ncurr
            # @debug format("\t{:6.4f}({}) < {} < {:6.4f}({})", nprev, k-1, C, ncurr, k)
            return (βs[:,k-1]*(ncurr-C) + βs[:,k]*(C-nprev))/(ncurr-nprev)
        end
        nprev = ncurr
    end
    # @debug format("\t{:6.4f}({}) < {}", nprev, K, C)
    return βs[:, end]
end

getdist(::Type{MSELoss}) = Normal()
getdist(::Type{LogitLoss}) = Binomial()

################################################################## GLMNet
using GLMNet

_gety(prob::Problem) = copy(prob.y)
_gety(prob::Problem{LogitLoss}) = [1-prob.y prob.y]

function lassoonA_glmnet!(ens::FWEnsemble, prob::Problem{L}, C) where L
  Z = getZ(ens)
  path = glmnet!(Z, _gety(prob), getdist(L), intercept=false, standardize=false)
  updateβ!(ens, getβfromC(path.betas, C))
end

################################################################## Lasso.jl
# using Lasso
#
# function lassoonA_lassojl(prob::Problem{L}, ens::Ensemble, C) where L
#     lfit = fit(LassoPath, Z, prob.y, getdist(L),
#                 intercept=false, standardize=false)
#     β = getβfromC(lfit.coefs, C)
# end


text:
#### Creating a small symbolic system for linear sums of the qnmfunctions
#### This is so we can write the derivative as a sum over multiple qnmfunctions

struct LinearCombinationOf{T <: CallableAtom} <: Callable
    dict::Dict{T,Complex{Float64}}
    constant::Complex{Float64}
end

LinearCombinationOf(x;c = Complex(0.0)) = LinearCombinationOf(x,c)

#(x::IdentityFunction)(s) = Complex(1.0)
#(x::HeunConfluentRadial)(s) = s^2 + 2s
(x::LinearCombinationOf{T})(s...) where T = sum(k(s...)*v for (k,v) ∈ x.dict) + x.constant

function postprocess(x::LinearCombinationOf{T}) where T
    if all(x == zero(x) for x ∈ values(x.dict))
        return LinearCombinationOf(Dict(first(x.dict)), first(values(x.dict)))
    end
    LinearCombinationOf(filter(a->a.second!=zero(a.second),x.dict),x.constant)
end

import Base.+, Base.-, Base.*
#const HCRNum = Union{HeunConfluentRadial, Number}
const LC = LinearCombinationOf{HeunConfluentRadial}
const HCR = HeunConfluentRadial
const LCHCR = Union{HeunConfluentRadial,LinearCombinationOf{HeunConfluentRadial}}

GetDict(x::Number) = Dict()
GetDict(x::CallableAtom; val = 1.0) = Dict(x => convert(Complex{Float64},val))
GetDict(x::LinearCombinationOf{T}) where T = x.dict

constant(x::Number) = x
constant(x::CallableAtom) = Complex(0.0)
constant(x::LinearCombinationOf{T}) where T = x.constant

## Overload Plus

function (+)(x::CallableAtom)
    LinearCombinationOf(GetDict(x))
end

function (+)(x::Callable,y::Callable)
    z = LinearCombinationOf(merge(+,GetDict(x),GetDict(y)),constant(x)+constant(y))
    postprocess(z)
end

function (+)(x::Number,y::Callable)
    LinearCombinationOf(GetDict(y),constant(x)+constant(y))
end

(+)(y::Callable, x::Number) = x+y

## Overload Minus

function (-)(x::CallableAtom)
    LinearCombinationOf(GetDict(x; val = -1.0))
end

function (-)(x::Callable,y::Callable)
    z = LinearCombinationOf(merge(+,GetDict(x),Dict((k,-v) for (k,v) ∈ GetDict(y))),constant(x)+constant(y))
    postprocess(z)
end

function (-)(x::Number,y::Callable)
    LinearCombinationOf(GetDict(y),constant(x)-constant(y))
end

function (-)(y::Callable,x::Number)
    LinearCombinationOf(GetDict(y),constant(y)-constant(x))
end

## Overload Times

function (*)(x::Number,y::Callable)
    ydict = GetDict(y);
    z = LinearCombinationOf(Dict((k,v*x) for (k,v) ∈ GetDict(y)),constant(y)*x)
    postprocess(z)
end

(*)(y::Callable,x::Number) = x*y

## Overload Show

function Base.show(io::IO, ψ::LinearCombinationOf{T}) where {T}
    N = length(ψ.dict); n = 1;
    print(io,"   ")
    for (k,v) in ψ.dict
        print(io,"(",v,")")
        print(io,k)
        if n < N
            println(io," + ")
        end
        n += 1
    end
end


text:
"""
    Transducer(o::OnlineStat)

Use an `OnlineStat` as a stateful transducer.

It is implemented as:

```julia
Scan(fit!, CopyInit(o)) |> Map(value)
```

# Examples
```jldoctest
julia> using Transducers
       using OnlineStats: Mean

julia> collect(Transducer(Mean()), 1:4)
4-element Array{Float64,1}:
 1.0
 1.5
 2.0
 2.5
```
"""
Transducer(o::OnlineStatsBase.OnlineStat) =
    Scan(OnlineStatsBase.fit!, CopyInit(o)) |> Map(OnlineStatsBase.value)
# TODO: implement `combine`


text:
using NCDatasets
using Interpolations
using Makie

#####
##### Read data
#####

ds = NCDataset("dry_convection.nc")

nx, ny, nz, nt = size(ds["ρe"])

ρe₀ = Array(ds["ρe₀"]);
ρe  = Array(ds["ρe"]);
ρe′ = similar(ρe);

for n in 1:length(ds["time"])
    @. ρe′[:, :, :, n] = ρe[:, :, :, n] - ρe₀;
end

#####
##### Interpolate to a finer grid with S::Int > 1
#####

function get_data(n, S=1)
    if S == 1
        return ρe′[:, :, :, n]
    else
        itp = interpolate(ρe′[:, :, :, n], BSpline(Cubic(Periodic(OnCell()))));
        n = size(ρe′[:, :, :, n])
        N = S .* n
        ℑρe′ = zeros(N)

        i_inds = range(1, n[1], length=N[1])
        j_inds = range(1, n[2], length=N[2])
        k_inds = range(1, n[3], length=N[3])

        i_inds = reshape(i_inds, (N[1], 1, 1))
        j_inds = reshape(j_inds, (1, N[2], 1))
        k_inds = reshape(k_inds, (1, 1, N[3]))

        @. ℑρe′ = itp(i_inds, j_inds, k_inds)

        return ℑρe′
    end
end

#####
##### Plot!
#####

α(ξ) = ξ  # Opacity/alpha along the cmap (0 <= ξ <= 1)

cmap_rgb = to_colormap(:ice)
A = α.(range(0, 1, length=length(cmap_rgb)))
cmap_rgba = RGBAf0.(cmap_rgb, A)

time_index = Node(1)
data = @lift get_data($time_index, 3)

function make_movie(scene)
    n_frames = nt
    θ = 2π / 360

    record(scene, "dry_convection_oceananigans_makie.mp4", 1:n_frames, framerate=30) do n
        @info "frame $n/$n_frames"
        time_index[] = n
        n == 1 && zoom!(scene, (0, 0, 0), -1.25, false)
        rotate_cam!(scene, (θ, 0, 0))
    end
end

scene = volume(0..nx, 0..ny, 0..nz, data, colorrange=(0, 400), colormap=cmap_rgba, algorithm=:absorption, absorption=10.0f0,
               backgroundcolor=cmap_rgb[2], show_axis=false, resolution = (1920, 1080))

make_movie(scene)


text:
module MDInclude

# Imports

using CommonMark

# Exports

export @mdinclude, mdinclude, markdown, isvalid, source, expression, setmodule, capture

# User Interface

"""
    @mdinclude(path, [config])

Include the given `path` in the current module. A configuration object may be
passed as well to control the behaviour of the file inclusion.
"""
macro mdinclude(args...)
    :(mdinclude(@__MODULE__, $(map(esc, args)...)))
end

"""
    mdinclude([module], path, [config])

Similar to `@mdinclude`, but allows an optional module to be given where the
`path` will be evaluated.
"""
function mdinclude(m::Module, path::AbstractString, config=nothing)
    path, prev = Base._include_dependency(m, path)
    ans = nothing
    for (node, enter) in markdown(config, path)
        if enter && isvalid(config, node)
            line = node.sourcepos[1][1]
            newm = setmodule(config, node, m)
            ans = capture(config) do
                _include_string(newm, source(config, node), path, prev) do expr
                    expression(config, _walk(x -> fixline(x, line), expr))
                end
            end
        end
    end
    return ans
end

# Configuration Interface

"""
    ast = markdown(config, path)

Returns the parsed AST for the given `path`. Can be used to customise the
markdown parser extensions.
"""
markdown(config, path) = open(Parser(), path)

"""
    valid = isvalid(config, node)

Checks whether the `node` should be `included` in the file's source code or
not. Returns either `true` or `false`.
"""
isvalid(config, node) = node.t isa CommonMark.CodeBlock && node.t.info == "julia"

"""
    src = source(config, src)

Returns the literal source that should be `included` for each valid AST node.
Can be used to perform textual transformations prior to Julia parsing the text.
"""
source(config, node) = node.literal

"""
    ex = expression(config, ex)

Returns the expression parsed by the Julia parser prior to being evaluated in
the module.
"""
expression(config, ex) = ex

"""
    mod = setmodule(config, node, default)

Select the module to evaluate the given node in. This allows for custom module
choices to be made based on the contents of the `node`. `default` is the
default module provided to `mdinclude`.
"""
setmodule(config, node, default) = default

"""
    value = capture(f, config)

Wraps the evaluation of each code block. Can be used to capture the `stdout`
and `stderr` streams during the evaluation as well as the value of the
evaluation. Should return the value of evaluating `f`.
"""
capture(f, config) = f()

# Utilities

# Traverse the expression and call `f` on each node.
function _walk(f, ex::Expr)
    for (nth, x) in enumerate(ex.args)
        ex.args[nth] = _walk(f, x)
    end
    return ex
end
_walk(f, @nospecialize(other)) = f(other)

# Based on `my_include_string` from NBInclude.jl.
function _include_string(mapexpr, m, txt, path, prev)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    try
        @static if VERSION < v"1.5"
            Core.eval(m, mapexpr(Base.parse_input_line(txt; filename=path)))
        else
            include_string(mapexpr, m, txt, path)
        end
    finally
        prev ≡ nothing ? delete!(tls, :SOURCE_PATH) : tls[:SOURCE_PATH] = prev
    end
end

# Adjust line numbers to match the source markdown file.
fixline(x, line) = x isa LineNumberNode ? LineNumberNode(x.line + line, x.file) : x

end # module


text:
@testset "rsoc problems" begin
    @testset "rsoc1_test" begin
    #=
    primal
        min 0a + 0b - 1x - 1y
    s.t.
        a    == 1/2 :w_2
        b    == 1   :w_1
        2a*b >= x^2+y^2  :w_3, w_4, w_5, w_6

   dual
        max w_2 + (1/2)w_1
    s.t.
        w_5 == -1
        w_6 == -1
        w_1 + w_3 == 0
        w_2 + w_4 == 0
        (w_3, w_4, w_5, w_6) \in RotatedSecondOrderCone
   =#
       primal_model = rsoc1_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 6
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.RotatedSecondOrderCone)
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.RotatedSecondOrderCone}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 4
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [0.5; 1.0]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == -1.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == -1.0
       eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       @test MOI.coefficient.(eq_con3_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con3_fun) == 0.0
       @test MOI.constant(eq_con3_set) == 0.0
       eq_con4_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(5))
       eq_con4_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(5))
       @test MOI.coefficient.(eq_con4_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con4_fun) == 0.0
       @test MOI.constant(eq_con4_set) == 0.0

       rsoc_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.RotatedSecondOrderCone}(1))
       @test rsoc_con.variables == VI.(1:4)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{SVF, MOI.EqualTo{Float64}}(3)] == [VI(5)]
       @test primal_con_dual_var[CI{SVF, MOI.EqualTo{Float64}}(4)] == [VI(6)]
       @test primal_con_dual_var[CI{VVF, MOI.RotatedSecondOrderCone}(1)] == VI.(1:4)

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
       @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
       @test primal_var_dual_con[VI(4)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(5)
   end

   @testset "rsoc2_test" begin
    #=
    primal
        min 0a + 0b - 1x - 1y
    s.t.
        a    == 1/2
        b    == 1
        2a*b >= x^2+y^2

    dual
        max w_2 + (1/2)w_1
    s.t.
        w_5 == -1
        w_6 == -1
        w_1 + w_3 == 0
        w_2 + w_4 == 0
        (w_3, w_4, w_5, w_6) \in RotatedSecondOrderCone
   =#
       primal_model = rsoc2_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.RotatedSecondOrderCone)
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.RotatedSecondOrderCone}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 2
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [-1.0; -0.5]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == -1.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == -1.0

       rsoc_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.RotatedSecondOrderCone}(1))
       @test rsoc_con.variables == VI.(1:4)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{VAF{Float64}, MOI.RotatedSecondOrderCone}(1)] == VI.(1:4)

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
   end
end