text:
@testset "Grids" begin

    # with a shift to the grid element, check if produce the correct floored index
    function check(grid, range, shift, idx_shift)
        for i in range
            @test(floor(grid, grid[i] + shift) == i + idx_shift)
            # if floor(grid, grid[i] + shift) != i + idx_shift
            #     return false
            # end
        end
        return true
    end

    @testset "UniformGrid" begin
        uniform = Grid.Uniform{Float64,4}(0.0, 1.0, (true, true))
        @test floor(uniform, 0.0) == 1
        @test floor(uniform, uniform[1]) == 1

        δ = 1.0e-12
        check(uniform, 2:uniform.size - 1, δ, 0)
        check(uniform, 2:uniform.size - 1, -δ, -1)

        @test floor(uniform, uniform[end]) == uniform.size - 1
        @test floor(uniform, 1.0) == uniform.size - 1

        uniform = Grid.Uniform{Float64,2}(0.0, 1.0, (true, true))
        @test floor(uniform, 0.0) == 1
        @test floor(uniform, uniform[1]) == 1

        δ = 1.0e-12
        check(uniform, 2:uniform.size - 1, δ, 0)
        check(uniform, 2:uniform.size - 1, -δ, -1)

        @test floor(uniform, uniform[end]) == uniform.size - 1
        @test floor(uniform, 1.0) == uniform.size - 1
    end

    @testset "Log Grid for Tau" begin
        β = 10.0
        tau = Grid.tau(β, 0.2β, 8)
        @test floor(tau, 0.0) == 1
        @test floor(tau, tau[1]) == 1 # tau[1]=0^+ is special

        δ = 1.0e-12
        check(tau, 2:tau.size - 1, δ, 0)
        check(tau, 2:tau.size - 1, -δ, -1)

        @test floor(tau, tau[end]) == tau.size - 1
        @test floor(tau, β) == tau.size - 1
    end

    @testset "Log Grid for fermiK" begin
        kF, maxK, halfLife = 1.0, 3.0, 0.5
        K = Grid.fermiK(kF, maxK, halfLife, 16)
        # println(K.grid)
        @test floor(K, 0.0) == 1
        @test floor(K, K[1]) == 1

        δ = 1.0e-12
        check(K, 2:K.size - 1, δ, 0)
        check(K, 2:K.size - 1, -δ, -1)

        @test floor(K, K[end] - δ) == K.size - 1
        @test floor(K, K[end]) == K.size - 1
        @test floor(K, maxK) == K.size - 1
    end

    @testset "Log Grid for boseK" begin
        kF, maxK, halfLife = 1.0, 3.0, 0.5
        K = Grid.boseK(kF, maxK, halfLife, 16)
        # println(K.grid)
        @test floor(K, 0.0) == 1
        @test floor(K, K[1]) == 1

        δ = 1.0e-12
        check(K, 2:K.size - 1, δ, 0)
        check(K, 2:K.size - 1, -δ, -1)

        @test floor(K, K[end] - δ) == K.size - 1
        @test floor(K, K[end]) == K.size - 1
        @test floor(K, maxK) == K.size - 1
    end

end

@testset "UniLog Grids" begin
    function check(grid, range, shift, idx_shift)
        for i in range
            @test(floor(grid, grid[i] + shift) == i + idx_shift)
            # if floor(grid, grid[i] + shift) != i + idx_shift
            #     return false
            # end
        end
        return true
    end

    @testset "UniLog" begin
        EPS = 1e-9
        bound = @SVector[1.0, 3.0]
        init = 4
        minterval = 0.01
        M = 2
        N = 3
        isopen = @SVector[false,false]
        d2s = true
        g = Grid.UniLog{Float64}(bound, init, minterval, M, N, d2s,isopen)

        val = Grid._grid(g,init)
        @test Grid._floor(g,val-EPS) == init
        @test Grid._floor(g,val+EPS) == init
        for i = init+1:(M+1)*N+init
            val = Grid._grid(g,i)
            # @printf("%10.6f , %10.6f\n", i, val)
            # @printf("\t%10.6f , %10.6f\n", val-EPS, Grid._floor(g,val-EPS))
            # @printf("\t%10.6f , %10.6f\n", val+EPS, Grid._floor(g,val+EPS))
            @test Grid._floor(g,val-EPS) == i-1
            @test Grid._floor(g,val+EPS) == i
        end
    end

    @testset "UniLogs" begin
        EPS = 1e-9
        seg = 4
        bounds = @SVector[0.0,1.0,2.0,3.0]
        M=2
        N=3
        minterval = 0.01
        isopen = @SVector[false,false]
        g = Grid.UniLogs{Float64,(M+1)*N*seg+1,seg}(bounds,minterval,M,N,isopen,[true,true])

        val = g.grid[1]
        @test floor(g,val-EPS) == 1
        @test floor(g,val+EPS) == 1
        for i = 2:g.size-1
            val = g.grid[i]
            # @printf("%10.6f , %10.6f\n", i, val)
            # @printf("\t%10.6f , %10.6f\n", val-EPS, floor(g,val-EPS))
            # @printf("\t%10.6f , %10.6f\n", val+EPS, floor(g,val+EPS))
            @test floor(g,val-EPS) == i-1
            @test floor(g,val+EPS) == i
        end
        val = g.grid[g.size]
        @test floor(g,val-EPS) == g.size-1
        @test floor(g,val+EPS) == g.size-1

        for i =1:seg-1
            val = g.segment[i]
            # @printf("%10.6f , %10.6f\n", g.segindex[i], val)
            # @printf("\t%10.6f , %10.6f\n", val-EPS, floor(g,val-EPS))
            # @printf("\t%10.6f , %10.6f\n", val+EPS, floor(g,val+EPS))
            @test floor(g,val-EPS) == g.segindex[i]-1
            @test floor(g,val+EPS) == g.segindex[i]
        end
    end

    @testset "UniLog Grid for Tau" begin
        β = 10.0
        tau = Grid.tauUL(β, 0.001, 2,3)
        @test floor(tau, 0.0) == 1
        @test floor(tau, tau[1]) == 1 # tau[1]=0^+ is special

        δ = 1.0e-12
        check(tau, 2:tau.size - 1, δ, 0)
        check(tau, 2:tau.size - 1, -δ, -1)

        @test floor(tau, tau[end]) == tau.size - 1
        @test floor(tau, β) == tau.size - 1
    end

    @testset "UniLog Grid for fermiK" begin
        kF, maxK, minterval = 1.0, 3.0, 0.001
        K = Grid.fermiKUL(kF, maxK, minterval, 2,3)
        # println(K.grid)
        @test floor(K, 0.0) == 1
        @test floor(K, K[1]) == 1

        δ = 1.0e-12
        check(K, 2:K.size - 1, δ, 0)
        check(K, 2:K.size - 1, -δ, -1)

        @test floor(K, K[end] - δ) == K.size - 1
        @test floor(K, K[end]) == K.size - 1
        @test floor(K, maxK) == K.size - 1
    end

    @testset "UniLog Grid for boseK" begin
        kF, maxK, minterval = 1.0, 3.0, 0.001
        K = Grid.boseKUL(kF, maxK, minterval,2,3)
        # println(K.grid)
        @test floor(K, 0.0) == 1
        @test floor(K, K[1]) == 1

        δ = 1.0e-12
        check(K, 2:K.size - 1, δ, 0)
        check(K, 2:K.size - 1, -δ, -1)

        @test floor(K, K[end] - δ) == K.size - 1
        @test floor(K, K[end]) == K.size - 1
        @test floor(K, maxK) == K.size - 1
    end

end

@testset "Interpolate" begin
    function kernelFermiT(t, ω, β=1)
        return exp(-ω * t) / (1 + exp(-ω * β))
    end

    @testset "Linear1D" begin
        β = π
        N,M= 3, 4
        tgrid = Grid.tauUL(β, 0.01β, M, N)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = t
        data = zeros(tgrid.size)

        for (ti, t) in enumerate(tgrid.grid)
            data[ti] = f(t)
        end

        for ti = 1:tgrid.size - 1
            t = tgrid[ti] + 1.e-6
            fbar = Grid.linear1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) < fbar
            @test f(tgrid[ti + 1]) > fbar
        end
        for ti = 2:tgrid.size
            t = tgrid[ti] - 1.e-6
            fbar = Grid.linear1D(data, tgrid, t)
            @test abs(f(tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            @test f(tgrid[ti]) > fbar
            @test f(tgrid[ti - 1]) < fbar
        end

        t = tgrid[1] + eps(Float64)*1e3
        fbar = Grid.linear1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        t = tgrid[tgrid.size] - eps(Float64)*1e3
        fbar = Grid.linear1D(data, tgrid, t)
        @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt

        tlist = rand(10) * β
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            fbar = Grid.linear1D(data, tgrid, t)
            # println("$k, $t, $fbar, ", f(k, t))
            @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
        end
    end

    @testset "Linear2D" begin
        β, kF, maxK = 10.0, 1.0, 3.0
        Nt, Nk = 8, 7
        tgrid = Grid.tau(β, 0.2β, Nt)
        kgrid = Grid.boseK(kF, maxK, 0.2kF, Nk)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(k, t) = k + t * 1.1
        data = zeros((Nk, Nt))

        for (ti, t) in enumerate(tgrid.grid)
            for (ki, k) in enumerate(kgrid.grid)
                data[ki, ti] = f(k, t)
            end
        end

        for ti = 1:tgrid.size - 1
            for ki = 1:kgrid.size - 1
                t = tgrid[ti] + 1.e-6
                k = kgrid[ki] + 1.e-6
                fbar = Grid.linear2D(data, kgrid, tgrid, k, t)
                @test abs(f(kgrid[ki], tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
                @test f(kgrid[ki], tgrid[ti]) < fbar
                @test f(kgrid[ki], tgrid[ti + 1]) > fbar
                @test f(kgrid[ki + 1], tgrid[ti]) > fbar
                @test f(kgrid[ki + 1], tgrid[ti + 1]) > fbar
            end
        end

        for ti = 2:tgrid.size
            for ki = 2:kgrid.size
                t = tgrid[ti] - 1.e-6
                k = kgrid[ki] - 1.e-6
                fbar = Grid.linear2D(data, kgrid, tgrid, k, t)
                @test abs(f(kgrid[ki], tgrid[ti]) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
                @test f(kgrid[ki], tgrid[ti]) > fbar
                @test f(kgrid[ki], tgrid[ti - 1]) < fbar
                @test f(kgrid[ki - 1], tgrid[ti]) < fbar
                @test f(kgrid[ki - 1], tgrid[ti - 1]) < fbar
            end
        end

        tlist = rand(10) * β
        klist = rand(10) * maxK
        # println(tlist)

        for (ti, t) in enumerate(tlist)
            for (ki, k) in enumerate(klist)
                fbar = Grid.linear2D(data, kgrid, tgrid, k, t)
                # println("$k, $t, $fbar, ", f(k, t))
                @test abs(f(k, t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
            end
        end
    end

    @testset "TestInterpolation1D for tauUL" begin
        β, minterval = 1.0, 0.0001
        M, N = 11, 4
        tgrid1 = Grid.tauUL(β, minterval, M, N)
        tgrid2 = Grid.tauUL(β, minterval, 2M, 2N)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(t) = kernelFermiT(t/β,0.01)

        d_max, std = Grid.testInterpolation1D(f, tgrid1, tgrid2, true)
        println("Testing interpolation for tauUL grid")
        println("d_max=",d_max,"\t std=", std)
    end

    @testset "TestInterpolation1D for fermiKUL" begin
        M, N= 15, 8
        kF, maxK, minterval = 1.0, 3.0, 0.00001
        kgrid1 = Grid.fermiKUL(kF, maxK, minterval, M, N)
        kgrid2 = Grid.fermiKUL(kF, maxK, minterval, 2M, 2N)
        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
        f(k) = 1.0/(0.0004π^2+(k^2-kF^2)^2)

        d_max, std = Grid.testInterpolation1D(f, kgrid1, kgrid2, true)
        println("Testing interpolation for fermiKUL grid")
        println("d_max=",d_max,"\t std=", std)
    end

    @testset "Optimize UniLog grid" begin
        struct Para
            β::Float64
            τ_min::Float64
            k_min::Float64
            k_max::Float64
            kF::Float64

            function Para()
                return new(10.0, 0.001, 0.0001, 3.0, 1.0)
            end
        end

        para = Para()
        MN = 64

        println("Testing optimization of tauUL")
        f(t) = kernelFermiT(t/para.β,1e0*para.kF^2*para.β)
        M, N, d_max = Grid.optimizeUniLog(Grid.tauUL, para, MN, f)
        println(MN,"\t", M,"\t", N,"\t", d_max)

        tgrid1 = Grid.tauUL(para, M, N)
        tgrid2 = Grid.tauUL(para, 2M, 2N)
        d_max, std = Grid.testInterpolation1D(f, tgrid1, tgrid2, true)
        println("Testing interpolation for tauUL grid")
        println("d_max=",d_max,"\t std=", std)

        MN = 64
        println("Testing optimization of fermiKUL")
        f2(k) = 1.0/((2π/para.β)^2+(k^2-para.kF^2)^2)
        M, N, d_max= Grid.optimizeUniLog(Grid.fermiKUL, para, MN, f2)
        println(MN,"\t", M,"\t", N,"\t", d_max)

        kgrid1 = Grid.fermiKUL(para, M, N)
        kgrid2 = Grid.fermiKUL(para, 2M, 2N)
        d_max, std, index = Grid.testInterpolation1D(f2, kgrid1, kgrid2, true, true)
        println("Testing interpolation for fermiKUL grid")
        println("d_max=",d_max,"\t std=", std)
        # println("index=", index, "\t x=", kgrid2[index])
        # index1=(floor(kgrid1, kgrid2[index]))
        # println(kgrid1[index1],"\t",kgrid1[index1+1])
        # println(kgrid1.grid)

    end

end


text:
"""
$(TYPEDSIGNATURES)

Call the method matching `node.name` from [`tagmap`](@ref) if that mapping exists,
otherwise parse as [`unclaimed_label`](@ref) wrapped in a [`PnmlLabel`](@ref).
"""
function parse_node end
parse_node(node::XMLNode; kw...) = parse_node(node, PnmlCore(); kw...)
function parse_node(node::XMLNode, pntd::PNTD; kw...) where {PNTD <: PnmlType}
    if haskey(tagmap, nodename(node))
        return tagmap[nodename(node)](node, pntd; kw...) # Various types returned here.
    else
        return PnmlLabel(unclaimed_label(node, pntd; kw...), node)
    end
end

#TODO test pnml_namespace
"""
$(TYPEDSIGNATURES)

Return namespace of `node. When `node` does not have a namespace return default value [`pnml_ns`](@ref)."
"""
function pnml_namespace(node::XMLNode; missing_ns_fatal=false, default_ns=pnml_ns)
    if EzXML.hasnamespace(node)
         return EzXML.namespace(node)
    else
        emsg = "$(nodename(node)) missing namespace"
        missing_ns_fatal===false ? @warn(emsg) : error(emsg)
        return default_ns
    end
end

"""
$(TYPEDSIGNATURES)

Start parse from the pnml root `node` of a well formed XML document.
Return a [`PnmlModel`](@ref)..
"""
function parse_pnml(node; kw...)
    nn = nodename(node)
    nn == "pnml" || error("element name wrong: $nn" )

    @assert haskey(kw, :reg)
    # Do not yet have a PNTD defined, so call parse_net directly.
    PnmlModel(parse_net.(allchildren("net", node), nothing; kw...),
              pnml_namespace(node),
              kw[:reg],
              node)
end

"""
$(TYPEDSIGNATURES)
Return a dictonary of the pnml net with keys matching their XML tag names.
"""
function parse_net(node, pntd=nothing; kw...)::PnmlNet
    nn = nodename(node)
    nn == "net" || error("element name wrong: $nn")
    EzXML.haskey(node, "id")   || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "type") || throw(MalformedException("$nn missing type", node))
    @assert haskey(kw, :reg)

    # Missing the page level in the pnml heirarchy causes nodes to be placed in :labels.
    # May result in undefined behavior and/or require ideosyncratic parsing.
    isempty(allchildren("page", node)) &&
         throw(MalformedException("$nn does not have any pages"))

    # Although the petri net type definition (pntd) must be attached to the <net> element,
    # it is allowed by this package to override that value.
    pntypedef = pnmltype(node["type"])
    if isnothing(pntd)
        pntd = pntypedef
    else
        @info """
        parse_net pntd set to $pntd
                    should be $pntypedef
        """
    end
    # Create a PnmlDict with keys for possible child tags.
    # Some keys have known/required values.
    # Optional key values are nothing for single object or empty vector when multiples.
    # Keys that have pural names usually have a vector value.
    # The 'graphics' key is one exception and has a single value.
    d = pnml_node_defaults(node, :tag => Symbol(nn),
                           :id => register_id!(kw[:reg], node["id"]),
                           :pages => Page[],
                           :declaration => Declaration())

    #!##pntd = d[:pntd] # Pass the PNTD down the parse tree with keyword arguments.

    # Go through children looking for expected tags, delegating common tags and labels.
    foreach(elements(node)) do child
        @match nodename(child) begin
            "page"         => push!(d[:pages], parse_page(child, pntd; kw...))

            # For nets and pages the <declaration> tag is optional
            # <declaration> ia a High-Level Annotation with a <structure> holding
            # zero or more <declarations>. Is complicated. You have been warned!
            "declaration"  => (d[:declaration] = parse_declaration(child, pntd; kw...))
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    PnmlNet(pntd, d[:id], d[:pages], d[:declaration], d[:name], ObjectCommon(d), node)
end
# Expected XML structure:
#    <declaration> <structure> <declarations> <namedsort id="weight" name="Weight"> ...

"""
$(TYPEDSIGNATURES)

PNML requires at least one page.
"""
function parse_page(node, pntd; kw...)
    nn = nodename(node)
    nn == "page" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag => Symbol(nn),
                           :id => register_id!(kw[:reg], node["id"]),
                           :places => Place[],
                           :trans => Transition[],
                           :arcs => Arc[],
                           :refP => RefPlace[],
                           :refT => RefTransition[],
                           :declaration => Declaration(),
                           :pages => Page[])

    foreach(elements(node)) do child
        @match nodename(child) begin
            "place"       => push!(d[:places], parse_place(child, pntd; kw...))
            "transition"  => push!(d[:trans], parse_transition(child, pntd; kw...))
            "arc"         => push!(d[:arcs], parse_arc(child, pntd; kw...))
            "referencePlace"      => push!(d[:refP], parse_refPlace(child, pntd; kw...))
            "referenceTransition" => push!(d[:refT], parse_refTransition(child, pntd; kw...))
            # See note above about declarations vs. declaration.
            "declaration" => (d[:declaration] = parse_declaration(child, pntd; kw...))
            "page"        => push!(d[:pages], parse_page(child, pntd; kw...))
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    Page(pntd, d[:id],
        d[:places], d[:refP],
        d[:trans], d[:refT],
        d[:arcs],
        d[:declaration], d[:pages], d[:name], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_place(node, pntd; kw...)
    nn = nodename(node)
    nn == "place" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag => Symbol(nn),
                           :id => register_id!(kw[:reg], node["id"]),
                           :marking => default_marking(pntd),
                           :type => default_sort(pntd)) # Different from net's.
    foreach(elements(node)) do child
        @match nodename(child) begin
            # Tags initialMarking and hlinitialMarking are mutually exclusive.
            "initialMarking"   => (d[:marking] = parse_initialMarking(child, pntd; kw...))
            "hlinitialMarking" => (d[:marking] = parse_hlinitialMarking(child, pntd; kw...))
            # Here type means `sort`. Re: Many-sorted algebra.
            "type"             => (d[:type] = parse_type(child, pntd; kw...))
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end

    Place(pntd, d[:id], d[:marking], d[:type], d[:name], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_transition(node, pntd; kw...)
    nn = nodename(node)
    nn == "transition" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag=>Symbol(nn),
                           :id=>register_id!(kw[:reg], node["id"]),
                           :condition=>default_condition(pntd))
    foreach(elements(node)) do child
        @match nodename(child) begin
            "condition"    => (d[:condition] = parse_condition(child, pntd; kw...))
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    Transition(pntd, d[:id], d[:condition], d[:name], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_arc(node, pntd; kw...)
    nn = nodename(node)
    nn == "arc" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    @assert haskey(node, "source")
    @assert haskey(node, "target")
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag=>Symbol(nn),
                           :id=>register_id!(kw[:reg], node["id"]),
                           :source=>Symbol(node["source"]),
                           :target=>Symbol(node["target"]),
                           :inscription=>default_inscription(pntd))
    foreach(elements(node)) do child
        @match nodename(child) begin
            # Mutually exclusive tags: inscription, hlinscription
            "inscription"    => (d[:inscription] = parse_inscription(child, pntd; kw...))
            "hlinscription"  => (d[:inscription] = parse_hlinscription(child, pntd; kw...))
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    Arc(pntd, d[:id], d[:source], d[:target], d[:inscription], d[:name], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_refPlace(node, pntd; kw...)
    nn = nodename(node)
    nn == "referencePlace" || error("element name wrong: $nn")
    EzXML.haskey(node, "id")  || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "ref") || throw(MalformedException("$(nn) missing ref attribute", node))
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag=>Symbol(nn),
                           :id=>register_id!(kw[:reg], node["id"]),
                           :ref=>Symbol(node["ref"]))
    foreach(elements(node)) do child
        @match nodename(child) begin
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    RefPlace(pntd, d[:id], d[:ref], d[:name], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_refTransition(node, pntd; kw...)
    nn = nodename(node)
    nn == "referenceTransition" || error("element name wrong: $nn")
    EzXML.haskey(node, "id") || throw(MissingIDException(nn, node))
    EzXML.haskey(node, "ref") || throw(MalformedException("$(nn) missing ref attribute", node))
    @assert haskey(kw, :reg)

    d = pnml_node_defaults(node, :tag=>Symbol(nn),
                           :id=>register_id!(kw[:reg], node["id"]),
                           :ref=>Symbol(node["ref"]))
    foreach(elements(node)) do child
        @match nodename(child) begin
            _ => parse_pnml_node_common!(d, child, pntd; kw...)
        end
    end
    RefTransition(pntd, d[:id], d[:ref], d[:name], ObjectCommon(d))
end

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

"""
$(TYPEDSIGNATURES)

Return the stripped string of nodecontent.
"""
function parse_text(node, pntd; kw...)
    nn = nodename(node)
    nn == "text" || error("$nn nodename wrong")
    string(strip(nodecontent(node)))
end

"""
$(TYPEDSIGNATURES)

Return [`Name`](@ref) holding text value and optional tool & GUI information.
"""
function parse_name(node, pntd; kw...)
    nn = nodename(node)
    nn == "name" || error("element nodename wrong")

    # Using firstchild or allchildren can cause parse_node to be passed nothing
    # for optional or missing child nodes.

    textnode = firstchild("text", node)
    if isnothing(textnode)
        @warn "$(nn) missing <text> element"
        # There are pnml files that break the rules & do not have a text element here.
        # Ex: PetriNetPlans-PNP/parallel.jl
        # Attempt to harvest content of <name> element instead of the child <text> element.
        # Assumes there are no other children elements.
        text = string(strip(nodecontent(node)))
    else
        text = string(strip(nodecontent(textnode)))
    end

    graphicsnode = firstchild("graphics", node)
    graphics = isnothing(graphicsnode) ? nothing : parse_graphics(graphicsnode, pntd; kw..., verbose=false)

    toolspecific = allchildren("toolspecific", node)
    tools = isempty(toolspecific) ? nothing : parse_toolspecific.(toolspecific, Ref(pntd); kw..., verbose=false)

    Name(; text, graphics, tools)
end

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

"""
$(TYPEDSIGNATURES)

Return [`Structure`](@ref) wrapping a `PnmlDict` holding a <structure>.
Should be inside of an label.
A "claimed" label usually elids the <structure> level (does not call this method).
"""
function parse_structure(node, pntd; kw...)
    nn = nodename(node)
    nn == "structure" || error("element name wrong: $nn")
    Structure(unclaimed_label(node, pntd; kw...), node)
end

#----------------------------------------------------------
#
# PNML annotation-label XML element parsers.
#
#----------------------------------------------------------

# Place Transition nets (PT-Nets) use only the text tag of a label for
# the meaning of marking and inscriptions.

"""
$(TYPEDSIGNATURES)
"""
function parse_initialMarking(node, pntd; kw...)
    nn = nodename(node)
    nn == "initialMarking" || error("element name wrong: $nn")
    d = pnml_label_defaults(node, :tag=>Symbol(nn),
        :value=>isempty(nodecontent(node)) ? _evaluate(default_marking(pntd)) :
                    number_value(strip(nodecontent(node)))
                    )

    foreach(elements(node)) do child
        @match nodename(child) begin
            # We extend to real numbers.
            "text" => (d[:value] = number_value(string(strip(nodecontent(child)))))
            _ => parse_pnml_label_common!(d, child, pntd; kw...)
        end
    end
    PTMarking(d[:value], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_inscription(node, pntd; kw...)
    nn = nodename(node)
    nn == "inscription" || error("element name wrong: $nn'")
    d = pnml_label_defaults(node, :tag=>Symbol(nn),
        :value=>isempty(nodecontent(node)) ? _evaluate(default_inscription(pntd)) :
                                    number_value(strip(nodecontent(node))))

    foreach(elements(node)) do child
        @match nodename(child) begin
            "text" => (d[:value] = number_value(string(strip(nodecontent(child)))))
            # Should not have a sturucture.
            _ => parse_pnml_label_common!(d, child, pntd; kw...)
        end
    end
    PTInscription(d[:value], ObjectCommon(d))
end

# High-Level Nets, includeing PT-HLPNG, are expected to use the structure child node to
# define the semantics of marking and inscriptions.

"""
$(TYPEDSIGNATURES)
"""
function parse_hlinitialMarking(node, pntd; kw...)
    nn = nodename(node)
    nn == "hlinitialMarking" || error("element name wrong: $nn")
    d = pnml_label_defaults(node, :tag=>Symbol(nn),
                            :text=>nothing,
                            :structure=>nothing)
    foreach(elements(node)) do child
        @match nodename(child) begin
            "structure" => (d[:structure] =
                haselement(child) ? parse_term(firstelement(child), pntd; kw...) :
                !isempty(nodecontent(child)) ? number_value(strip(nodecontent(child))) :
                default_tmarking(pntd)())
            _ => parse_pnml_label_common!(d, child, pntd; kw...)
        end
    end
    HLMarking(d[:text], d[:structure], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)
"""
function parse_hlinscription(node, pntd; kw...)
    @debug node
    nn = nodename(node)
    nn == "hlinscription" || error("element name wrong: $nn'")
    d = pnml_label_defaults(node, :tag=>Symbol(nn))
    foreach(elements(node)) do child
        @match nodename(child) begin
        "structure" => (d[:structure] =
                haselement(child) ? parse_term(firstelement(child), pntd; kw...) :
                !isempty(nodecontent(child)) ? number_value(strip(nodecontent(child))) :
                default_inscription(pntd)())
     _ => parse_pnml_label_common!(d, child, pntd; kw...)
        end
    end
    HLInscription(d[:text], d[:structure], ObjectCommon(d))
end

"""
$(TYPEDSIGNATURES)

Label of transition nodes.

# Details

Condition is defined by the ISO Specification as a High-level Annotation,
meaning it has <text> and <structure> elements. With all meaning in the element
that the <structure> holds.

We extend Condition by allowing <structure> to have _context_ instead of of a child element.
This content is treated as either an Integer or Float64.

A Condition should evaluate to a boolean. We defer that evaluation to a higher level.
See [`AbstractTerm`](@ref).
"""
function parse_condition(node, pntd; kw...)
    @debug node
    nn = nodename(node)
    nn == "condition" || error("element name wrong: $nn")
    d = pnml_label_defaults(node, :tag=>Symbol(nn))
    foreach(elements(node)) do child
        @match nodename(child) begin
            "structure" => (d[:structure] =
                    haselement(child) ? parse_term(firstelement(child), pntd; kw...) :
                    !isempty(nodecontent(child)) ? number_value(strip(nodecontent(child))) :
                    default_condition(pntd)())
            _ => parse_pnml_label_common!(d, child, pntd; kw...)
        end
    end
    Condition(d[:text], d[:structure], ObjectCommon(d))
end

#---------------------------------------------------------------------
#TODO Will unclaimed_node handle this?
"""
$(TYPEDSIGNATURES)

Should not often have a '<label>' tag, this will bark if one is found.
Return minimal PnmlDict holding (tag,node), to defer parsing the xml.
"""
function parse_label(node, pntd; kw...)
    nn = nodename(node)
    nn == "label" || error("element name wrong: $nn")
    @warn "parse_label '$(node !== nothing && nn)'"
    PnmlDict(:tag=>Symbol(nn), :xml=>node) # Always add xml because this is unexpected.
end

ext: jl


text:
tol = 1e-11

# Check if besselj(nu, besselroots(nu, n) ) is small:

for nu = 0.:0.1:5.
    n = 10
    @test norm( besselj(nu, besselroots(nu, n) ), Inf ) < tol

end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

function Geometries(T::DataType, geometries::Vector, inputunit::Unitful.Units)::Vector{AbstractGeometry}
    return [Geometry(T, Val{Symbol(geometry["type"])}(), Dict{Any,Any}(geometry)::Union{Dict{String,Any},Dict{Any,Any}}, inputunit) for geometry in geometries]
end


text:
vtrait(X) = X |> trait |> Val

const REQUIRE = "(requires MLJBase to be loaded)"

errlight(s) = throw(InterfaceError("Only `MLJModelInterface` is loaded. " *
                                    "Import `MLJBase` in order to use `$s`."))

# ------------------------------------------------------------------------
# categorical, note: not exported to avoid clashes; this is fine because
# MLJBase loads CategoricalArrays and MLJ interfaces should use qualified
# statements.

categorical(a...; kw...) = categorical(get_interface_mode(), a...; kw...)

categorical(::LightInterface, a...; kw...) = errlight("categorical")

# ------------------------------------------------------------------------
# matrix

"""
    matrix(X; transpose=false)

If `X <: AbstractMatrix`, return `X` or `permutedims(X)` if `transpose=true`.
If `X` is a Tables.jl compatible table source, convert `X` into a `Matrix`.

"""
matrix(X; kw...) = matrix(get_interface_mode(), vtrait(X), X; kw...)

matrix(::Mode, ::Val{:other}, X::AbstractMatrix; transpose=false) =
    transpose ? permutedims(X) : X

matrix(::Mode, ::Val{:other}, X; kw...) =
    throw(ArgumentError("Function `matrix` only supports AbstractMatrix or " *
                        "containers implementing the Tables interface."))

matrix(::LightInterface, ::Val{:table}, X; kw...) = errlight("matrix")

# ------------------------------------------------------------------------
# int

"""
   int(x; type=nothing)

The positional integer of the `CategoricalString` or `CategoricalValue` `x`, in
the ordering defined by the pool of `x`. The type of `int(x)` is the reference
type of `x`.

Not to be confused with `x.ref`, which is unchanged by reordering of the pool
of `x`, but has the same type.

    int(X::CategoricalArray)
    int(W::Array{<:CategoricalString})
    int(W::Array{<:CategoricalValue})

Broadcasted versions of `int`.

    julia> v = categorical([:c, :b, :c, :a])
    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001

See also: [`decoder`](@ref).
"""
function int(x; type::Union{Nothing,Type{T}}=nothing) where T <: Real
    type === nothing && return int(get_interface_mode(), x)
    return convert.(T, int(get_interface_mode(), x))
end

int(::LightInterface, x) = errlight("int")

# ------------------------------------------------------------------------
# classes

"""
    classes(x)

All the categorical elements with the same pool as `x` (including
`x`), returned as a list, with an ordering consistent with the pool.
Here `x` has `CategoricalValue` or `CategoricalString` type, and
`classes(x)` is a vector of the same eltype. Note that `x in
classes(x)` is always true.

Not to be confused with `levels(x.pool)`. See the example below.

    julia>  v = categorical([:c, :b, :c, :a])
    4-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :c
     :b
     :c
     :a

    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

    julia> x = v[4]
    CategoricalArrays.CategoricalValue{Symbol,UInt32} :a

    julia> classes(x)
    3-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :a
     :b
     :c

    julia> levels(x.pool)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

"""
classes(x) = classes(get_interface_mode(), x)

classes(::LightInterface, x) = errlight("classes")

# ------------------------------------------------------------------------
# schema

"""
    schema(X)

If `X` is tabular, inspect the column types and scitypes, otherwise return
`nothing`.
"""
schema(X; kw...) = schema(get_interface_mode(), vtrait(X), X; kw...)

schema(::Mode, ::Val{:other}, X; kw...) = nothing

schema(::LightInterface, ::Val{:other}, X; kw...) = errlight("schema")

schema(::LightInterface, ::Val{:table}, X; kw...) = errlight("schema")

# ------------------------------------------------------------------------
# istable

"""
    istable(X)

Return true if `X` is tabular.
"""
istable(X) = istable(get_interface_mode(), vtrait(X))

istable(::Mode, ::Val{:other}) = false

istable(::Mode, ::Val{:table}) = true

# ------------------------------------------------------------------------
# decoder

"""
    d = decoder(x)

A callable object for decoding the integer representation of a
`CategoricalString` or `CategoricalValue` sharing the same pool as
`x`. (Here `x` is of one of these two types.) Specifically, one has
`d(int(y)) == y` for all `y in classes(x)`. One can also call `d` on
integer arrays, in which case `d` is broadcast over all elements.

    julia> v = categorical([:c, :b, :c, :a])
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001
    julia> d = decoder(v[3])
    julia> d(int(v)) == v
    true

*Warning:* It is *not* true that `int(d(u)) == u` always holds.

See also: [`int`](@ref), [`classes`](@ref).
"""
decoder(x) = decoder(get_interface_mode(), x)

decoder(::LightInterface, x) = errlight("decoder")

# ------------------------------------------------------------------------
# table

"""
    table(columntable; prototype=nothing)

Convert a named tuple of vectors or tuples `columntable`, into a table
of the "preferred sink type" of `prototype`. This is often the type of
`prototype` itself, when `prototype` is a sink; see the Tables.jl
documentation. If `prototype` is not specified, then a named tuple of
vectors is returned.

    table(A::AbstractMatrix; names=nothing, prototype=nothing)

Wrap an abstract matrix `A` as a Tables.jl compatible table with the
specified column `names` (a tuple of symbols). If `names` are not
specified, `names=(:x1, :x2, ..., :xn)` is used, where `n=size(A, 2)`.

If a `prototype` is specified, then the matrix is materialized as a table of
the preferred sink type of `prototype`, rather than wrapped. Note that if
`prototype` is *not* specified, then `matrix(table(A))` is essentially a no-op.
"""
table(X; kw...) = table(get_interface_mode(), X; kw...)

table(::LightInterface, X; kw...) = errlight("table")

# ------------------------------------------------------------------------
# nrows, select, selectrows, selectcols

"""
    nrows(X)

Return the number of rows for a table, abstract vector or matrix `X`.
"""
nrows(X) = nrows(get_interface_mode(), vtrait(X), X)

nrows(::Mode, ::Val{:other}, X::AbstractVecOrMat) = size(X, 1)
nrows(::Mode, ::Val{:other}, X::Nothing) = 0

nrows(::Mode, ::Val{:other}, X) =
    throw(ArgumentError("Function `nrows` only supports AbstractVector or " *
                        "AbstractMatrix or containers implementing the " *
                        "Tables interface."))

nrows(::LightInterface, ::Val{:table}, X) = errlight("table")

"""
    selectrows(X, r)

Select single or multiple rows from a table, abstract vector or matrix
`X`. If `X` is tabular, the object returned is a table of the
preferred sink type of `typeof(X)`, even if only a single row is
selected.

If the object is neither a table, abstract vector or matrix, `X` is
returned and `r` is ignored.

"""
selectrows(X, r) = selectrows(get_interface_mode(), vtrait(X), X, r)

# fall-back is to return object, ignoring vector of row indices, `r`:
selectrows(::Mode, ::Val{:other}, X::Any, r) = X

selectrows(::Mode, ::Val{:other}, X::AbstractVector, r)          = X[r]
selectrows(::Mode, ::Val{:other}, X::AbstractVector, r::Integer) = X[r:r]
selectrows(::Mode, ::Val{:other}, X::AbstractVector, ::Colon)    = X

selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, r)          = X[r, :]
selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, r::Integer) = X[r:r, :]
selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, ::Colon)    = X

selectrows(::LightInterface, ::Val{:table}, X, r; kw...) =
    errlight("selectrows")

"""
    selectcols(X, c)

Select single or multiple columns from a matrix or table `X`. If `c`
is an abstract vector of integers or symbols, then the object returned
is a table of the preferred sink type of `typeof(X)`. If `c` is a
*single* integer or column, then an `AbstractVector` is returned.

"""
selectcols(X, c) = selectcols(get_interface_mode(), vtrait(X), X, c)

selectcols(::Mode, ::Val{:other}, ::Nothing, c) = nothing

selectcols(::Mode, ::Val{:other}, X::AbstractMatrix, r)       = X[:, r]
selectcols(::Mode, ::Val{:other}, X::AbstractMatrix, ::Colon) = X

selectcols(::Mode, ::Val{:other}, X, r) =
    throw(ArgumentError("Function `selectcols` only supports AbstractMatrix " *
                        "or containers implementing the Tables interface."))

selectcols(::LightInterface, ::Val{:table}, X, c; kw...) =
    errlight("selectcols")

"""
    select(X, r, c)

Select element(s) of a table or matrix at row(s) r and column(s) c. An object
of the sink type of X (or a matrix) is returned unless c is a single integer or
symbol. In that case a vector is returned, unless r is a single integer, in
which case a single element is returned.

See also: [`selectrows`](@ref), [`selectcols`](@ref).
"""
select(X, r, c) = select(get_interface_mode(), vtrait(X), X, r, c)

# only used here to denote "group of indices"
const MIdx = Union{AbstractArray,Colon}

select(::Mode, ::Val, X, r::MIdx, c)       = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r, c::MIdx)       = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r::MIdx, c::MIdx) = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r, c) = _squeeze(selectcols(selectrows(X, r), c))

_squeeze(::Nothing) = nothing
_squeeze(v) = first(v)

# ------------------------------------------------------------------------
# UnivariateFinite

const UNIVARIATE_FINITE_DOCSTRING =
"""
    UnivariateFinite(support,
                     probs;
                     pool=nothing,
                     augmented=false,
                     ordered=false)

Construct a discrete univariate distribution whose finite support is
the elements of the vector `support`, and whose corresponding
probabilities are elements of the vector `probs`. Alternatively,
construct an abstract *array* of `UnivariateFinite` distributions by
choosing `probs` to be an array of one higher dimension than the array
generated.

Here the word "probabilities" is an abuse of terminology as there is
no requirement that probabilities actually sum to one, only that they
be non-negative. So `UnivariateFinite` objects actually implement
arbitrary non-negative measures over finite sets of labelled points. A
`UnivariateDistribution` will be a bona fide probability measure when
constructed using the `augment=true` option (see below) or when
`fit` to data.

Unless `pool` is specified, `support` should have type
 `AbstractVector{<:CategoricalValue}` and all elements are assumed to
 share the same categorical pool, which may be larger than `support`.

*Important.* All levels of the common pool have associated
probabilities, not just those in the specified `support`. However,
these probabilities are always zero (see example below).

If `probs` is a matrix, it should have a column for each class in
`support` (or one less, if `augment=true`). More generally, `probs`
will be an array whose size is of the form `(n1, n2, ..., nk, c)`,
where `c = length(support)` (or one less, if `augment=true`) and the
constructor then returns an array of `UnivariateFinite` distributions
of size `(n1, n2, ..., nk)`.

```
using CategoricalArrays
v = categorical([:x, :x, :y, :x, :z])

julia> UnivariateFinite(classes(v), [0.2, 0.3, 0.5])
UnivariateFinite{Multiclass{3}}(x=>0.2, y=>0.3, z=>0.5)

julia> d = UnivariateFinite([v[1], v[end]], [0.1, 0.9])
UnivariateFinite{Multiclass{3}(x=>0.1, z=>0.9)

julia> rand(d, 3)
3-element Array{Any,1}:
 CategoricalArrays.CategoricalValue{Symbol,UInt32} :z
 CategoricalArrays.CategoricalValue{Symbol,UInt32} :z
 CategoricalArrays.CategoricalValue{Symbol,UInt32} :z

julia> levels(d)
3-element Array{Symbol,1}:
 :x
 :y
 :z

julia> pdf(d, :y)
0.0
```

### Specifying a pool

Alternatively, `support` may be a list of raw (non-categorical)
elements if `pool` is:

- some `CategoricalArray`, `CategoricalValue` or `CategoricalPool`,
  such that `support` is a subset of `levels(pool)`

- `missing`, in which case a new categorical pool is created which has
  `support` as its only levels.

In the last case, specify `ordered=true` if the pool is to be
considered ordered.

```
julia> UnivariateFinite([:x, :z], [0.1, 0.9], pool=missing, ordered=true)
UnivariateFinite{OrderedFactor{2}}(x=>0.1, z=>0.9)

julia> d = UnivariateFinite([:x, :z], [0.1, 0.9], pool=v) # v defined above
UnivariateFinite(x=>0.1, z=>0.9) (Multiclass{3} samples)

julia> pdf(d, :y) # allowed as `:y in levels(v)`
0.0

v = categorical([:x, :x, :y, :x, :z, :w])
probs = rand(100, 3)
probs = probs ./ sum(probs, dims=2)
julia> UnivariateFinite([:x, :y, :z], probs, pool=v)
100-element UnivariateFiniteVector{Multiclass{4},Symbol,UInt32,Float64}:
 UnivariateFinite{Multiclass{4}}(x=>0.194, y=>0.3, z=>0.505)
 UnivariateFinite{Multiclass{4}}(x=>0.727, y=>0.234, z=>0.0391)
 UnivariateFinite{Multiclass{4}}(x=>0.674, y=>0.00535, z=>0.321)
   ⋮
 UnivariateFinite{Multiclass{4}}(x=>0.292, y=>0.339, z=>0.369)
```

### Probability augmentation

If `augment=true` the provided array is augmented by inserting
appropriate elements *ahead* of those provided, along the last
dimension of the array. This means the user only provides probabilities
for the classes `c2, c3, ..., cn`. The class `c1` probabilities are
chosen so that each `UnivariateFinite` distribution in the returned
array is a bona fide probability distribution.

---

    UnivariateFinite(prob_given_class; pool=nothing, ordered=false)

Construct a discrete univariate distribution whose finite support is
the set of keys of the provided dictionary, `prob_given_class`, and
whose values specify the corresponding probabilities.

The type requirements on the keys of the dictionary are the same as
the elements of `support` given above with this exception: if
non-categorical elements (raw labels) are used as keys, then
`pool=...` must be specified and cannot be `missing`.

If the values (probabilities) are arrays instead of scalars, then an
abstract array of `UnivariateFinite` elements is created, with the
same size as the array.

"""
UNIVARIATE_FINITE_DOCSTRING
UnivariateFinite(d::AbstractDict; kwargs...) =
    UnivariateFinite(get_interface_mode(), d; kwargs...)
UnivariateFinite(support::AbstractVector, probs; kwargs...) =
    UnivariateFinite(get_interface_mode(), support, probs; kwargs...)
UnivariateFinite(probs; kwargs...) =
    UnivariateFinite(get_interface_mode(), probs; kwargs...)
UnivariateFinite(::LightInterface, a...; kwargs...) =
    errlight("UnivariateFinite")


text:
#
# Exponential fit
#

struct SingleExponential
  a :: Float64
  b :: Float64
  c :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

struct MultipleExponential
  n :: Int64
  a :: Vector{Float64}
  b :: Vector{Float64}
  c :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

function sum_of_exps(x :: Real, p :: AbstractVector)
  n = div(length(p)-1,2)
  f = p[length(p)]
  for i in 1:n
    f = f + p[i]*exp(-x/p[n+i])
  end
  return f
end

function exp_model(X :: AbstractVector, p :: AbstractVector)
  f = Vector{Float64}(undef,length(X))
  for i in 1:length(X)
    f[i] = sum_of_exps(X[i],p)
  end
  return f
end

function sum_of_exps_const(x :: Real, p :: AbstractVector, c)
  n = div(length(p),2)
  f = c
  for i in 1:n
    f = f + p[i]*exp(-x/p[n+i])
  end
  return f
end

function exp_model_const(X :: AbstractArray{<:Real}, p :: AbstractVector, c)
  f = Vector{Float64}(undef,length(X))
  for i in 1:length(X)
    f[i] = sum_of_exps_const(X[i],p,c)
  end
  return f
end

"""
`fitexponential(x,y;n :: Int = 1)`  or  `fitexp(x,y;n :: Int = 1)`

Obtains single or multiexponential fits: ``y = a*exp(-x/b) + c`` or  ``y = sum(a[i]*exp(-x/b[i]) for i in 1:N) + c``

Lower and upper bounds can be optionall set, and the intercept `c` can set to be constant:

For single exponentials all parameters are scalars:

```
fitexp(x,y,l=lower(b=0.),u=upper(b=10.),c=5.)
```

For multiple exponentials, `a` and `b` bounds must be vectors of dimension `N`.

```
fitexp(x,y,n=2,l=lower(a=[0.,0.]),u=upper(b=[-100.,-5.]))
```

# Examples
```jldoctest
julia> x = sort(rand(10)); y = rand()*exp.(sort(rand(10)));

julia> fit = fitexp(x,y,l=lower(b=[0.,0.]),n=2)

 -------- Multiple-exponential fit -------------

 Equation: y = sum(a[i] exp(-x/b[i]) for i in 1:2) + c

 With: a = [6.60693727987886e-13, 0.6249999999993409]
       b = [0.02688289803014393, 0.5000000000002596]
       c = 0.37499999999999856

 Pearson correlation coefficient, R = 1.0
 Average square residue = 1.1639900380979497e-29

 Predicted Y: ypred = [1.0000000000000002, 0.4595845520228801...
 residues = [2.220446049250313e-16, -2.831068712794149e-15...

 -----------------------------------------------

```
"""
function fitexponential(X :: AbstractArray{<:Real}, Y :: AbstractArray{<:Real};
                        n :: Int = 1, c = nothing,
                        l :: lower = lower(), u :: upper = upper(),
                        options :: Options = Options())

  # Check data
  X, Y = checkdata(X,Y,options)
  if c == nothing
    # Set model
    model(x,p) = exp_model(x,p)
    # Number of exponentials
    # Set bounds
    if n == 1
      vars = [ VarType(:a,Number,1),
               VarType(:b,Number,1),
               VarType(:c,Nothing,1) ]
      lower, upper = setbounds(vars,l,u)
    else
      vars = [ VarType(:a,Vector,n),
               VarType(:b,Vector,n),
               VarType(:c,Nothing,1) ]
      lower, upper = setbounds(vars,l,u)
    end
    # Fit
    fit = find_best_fit(model,X,Y,2*n+1,options,lower,upper)
    # Analyze and return
    R = pearson(X,Y,model,fit)
    x, y, ypred = finexy(X,options.fine,model,fit)
    if n == 1
      return SingleExponential(fit.param[1],fit.param[2],fit.param[3],
                               R,x,y,ypred,fit.resid)
    else
      ind = collect(1:n)
      sort!(ind,by=i->fit.param[n+i])
      a = fit.param[1:n][ind]
      b = fit.param[n+1:2*n][ind]
      return MultipleExponential(n,a,b,fit.param[2*n+1],
                                 R,x,y,ypred,fit.resid)
    end
  else
    # Set model
    model_const(x,p) = exp_model_const(x,p,c)
    # Number of exponentials
    # Set bounds
    if n == 1
      vars = [ VarType(:a,Number,1),
               VarType(:b,Number,1),
               VarType(:c,Nothing,1) ]
      lower, upper = setbounds(vars,l,u)
    else
      vars = [ VarType(:a,Vector,n),
               VarType(:b,Vector,n),
               VarType(:c,Nothing,1) ]
      lower, upper = setbounds(vars,l,u)
    end
    lower = lower[1:2*n]
    upper = upper[1:2*n]
    # Fit
    fit = find_best_fit(model_const,X,Y,2*n,options,lower,upper)
    # Analyze and return
    R = pearson(X,Y,model_const,fit)
    x, y, ypred = finexy(X,options.fine,model_const,fit)
    if n == 1
      return SingleExponential(fit.param[1],fit.param[2],c,R,x,y,ypred,fit.resid)
    else
      ind = collect(1:n)
      sort!(ind,by=i->fit.param[n+i])
      a = fit.param[1:n][ind]
      b = fit.param[n+1:2*n][ind]
      return MultipleExponential(n,a,b,c,R,x,y,ypred,fit.resid)
    end
  end
end
fitexp = fitexponential

function Base.show( io :: IO, fit :: SingleExponential )
  println("")
  println(" ------------ Single Exponential fit ----------- ")
  println("")
  println(" Equation: y = a exp(-x/b) + c")
  println("")
  println(" With: a = ", fit.a)
  println("       b = ", fit.b)
  println("       c = ", fit.c)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println(" Average square residue = ", mean(fit.residues.^2))
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

function Base.show( io :: IO, fit :: MultipleExponential )
  println("")
  println(" -------- Multiple-exponential fit ------------- ")
  println("")
  println(" Equation: y = sum(a[i] exp(-x/b[i]) for i in 1:$(fit.n)) + c ")
  println("")
  println(" With: a = ", fit.a)
  println("       b = ", fit.b)
  println("       c = ", fit.c)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println(" Average square residue = ", mean(fit.residues.^2))
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitexp, fitexponential


text:
__precompile__()

module ImmutableNestedSetModel

export Node, AnyNode, root, prepare, parents, parentscount, NodeIterator,
  Tree, iter, StructItems, KeyedItems, Table, KeyedTable, StructTable

import DataStructures: OrderedDict

include("types/Node.jl")
include("types/Tree.jl")

end # module


text:
struct ConstrainedTimeInvariantLQR{T <: Number}
    """
    The predicted system
    """
    sys::AbstractStateSpace

    """
    The horizon length
    """
    N::Integer

    """
    The state weight matrix
    """
    Q::AbstractMatrix{T}

    """
    The Q weighting matrix taking into account the prestabilizing controller
    """
    Qₖ::AbstractMatrix{T}

    """
    The input weight matrix
    """
    R::AbstractMatrix{T}

    """
    The terminal state weight matrix
    """
    P::AbstractMatrix{T}

    """
    The cross-term weight matrix
    """
    S::AbstractMatrix{T}

    """
    A prestabilizing controller (if any) that forms A-BK
    """
    K::AbstractMatrix{T}

    """
    The stage state constraint coefficient matrix
    """
    E::AbstractMatrix{T}

    """
    The stage input constraint coefficient matrix
    """
    F::AbstractMatrix{T}

    """
    The stage right hand side of the constraints
    """
    g::AbstractVector{T}

    """
    The lower bounds on the state variables
    """
    xₗ::AbstractVector{Union{T, Missing}}

    """
    The upper bounds on the state variables
    """
    xᵤ::AbstractVector{Union{T, Missing}}

    """
    The lower bounds on the input variables
    """
    uₗ::AbstractVector{Union{T, Missing}}

    """
    The upper bounds on the input variables
    """
    uᵤ::AbstractVector{Union{T, Missing}}
end

function ConstrainedTimeInvariantLQR( sys::AbstractStateSpace, N::Integer, Q::AbstractNumOrUniform, R::AbstractNumOrUniform, P::AbstractNumOrUniformOrSymbol = :dare;
                                      S::AbstractNumOrMatrix = 0, K::AbstractNumOrMatrixOrSymbol = 0, E::AbstractNumOrMatrix = 0, F::AbstractNumOrMatrix = 0, g::AbstractNumOrVector = 0,
                                      xₗ::AbstractNumOrVector = -Inf, xᵤ::AbstractNumOrVector = Inf, uₗ::AbstractNumOrVector = -Inf, uᵤ::AbstractNumOrVector = Inf )
    T = promote_type( eltype( Q ), eltype( R ), eltype( F ), eltype( E ), eltype( g ), eltype( S ) )
    nₓ = nstates( sys )
    nᵤ = ninputs( sys )

    # Handle the UniformScaling cases
    if( typeof( Q ) <: UniformScaling )
        Q = to_matrix( T, Q, nₓ )
    else
        Q = to_matrix( T, Q )
    end

    if( typeof( R ) <: UniformScaling )
        R = to_matrix( T, R, nᵤ )
    else
        R = to_matrix( T, R )
    end

    # Verify arguments
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    if( !isposdef( R ) || !issymmetric( R ) )
        throw( DomainError( R, "R must be symmetric positive definite") )
    end

    if( !issymmetric( Q ) )
        throw( DomainError( Q, "Q must be symmetric" ) )
    end

    if( size( Q, 1 ) != nₓ )
        throw( DimensionMismatch( "Q must be square with dimension the same as the number of states" ) )
    end

    if( size( R, 1 ) != nᵤ )
        throw( DimensionMismatch( "R must be square with dimension the same as the number of rows" ) )
    end

    if( S != 0 )
        if( size( S, 1 ) != nₓ || size( S, 2 ) != nᵤ )
            throw( DimensionMismatch( "S must have $nₓ rows and $nᵤ columns" ) )
        end
        S = to_matrix( T, S )
    else
        S = zeros( T, nₓ, nᵤ )
    end

    if( typeof( K ) <: Symbol )
        if( K == :dlqr )
            _, _, K = ared( sys.A, sys.B, R, Q )
        else
            throw(  ArgumentError( "Unknown value $K for K. Allowed values are :dlqr or a $nₓ by $nᵤ matrix." ) )
        end

    elseif( K != 0 )
        if( size( K, 1 ) != nᵤ || size( K, 2 ) != nₓ )
            throw( DimensionMismatch( "K must have $nₓ rows and $nᵤ columns" ) )
        end
        K = to_matrix( T, K )

    else
        K = zeros( T, nᵤ, nₓ )
    end

    # Create the controlled state weighting matrix (this is just the normal weighting matrix if there is no prestabilization)
    Qₖ = Q + K'*R*K

    # Handle the final state weighting matrix
    if( typeof( P ) <: Symbol )
        if( P == :Q )
            P = Qₖ

        elseif( P == :dare )
            # This computation is done using the original weighting matrices
            P, = ared( sys.A, sys.B, R, Q )

        elseif( P == :dlyap )
            # This computation is done using the original weighting matrices
            P = lyapd( sys.A', Q )

        else
            throw( ArgumentError( "Unknown value $P for P. Allowed values are :Q, :dare, :dlyap or a $nₓ by $nₓ matrix." ) )
        end

    elseif( typeof( P ) <: UniformScaling )
        P = to_matrix( T, P, nₓ )

    else
        if( size( P, 1 ) != nₓ )
            throw( DimensionMismatch( "P must be square with dimension the same as the number of states" ) )
        elseif( !issymmetric(P) )
            throw( DomainError( P, "P must be symmetric" ) )
        end

        P = to_matrix( T, P )
    end

    # Verify the constraints when there are some
    if( g != 0 )
        ( size( E, 1 ) != size( F, 1 ) ) && throw( DimensionMismatch( "E and F must have the same number of rows" ) )
        ( size( g, 1 ) != size( F, 1 ) ) && throw( DimensionMismatch( "g must have $size( F, 1 ) rows" ) )
        ( size( E, 2 ) != nₓ ) && throw( DimensionMismatch( "E must have $nₓ columns" ) )
        ( size( F, 2 ) != nᵤ ) && throw( DimensionMismatch( "F must have $nᵤ columns" ) )

        E = to_matrix( T, E )
        F = to_matrix( T, F )
    else
        # Create an empty matrix for the constraints
        E = zeros( T, 1, nₓ )
        F = zeros( T, 1, nᵤ )
        g = zeros( T, 1 )
    end

    # Create the bounds for the state and inputs
    ( size( xₗ, 1 ) != 1 ) && ( size( xₗ, 1 ) != nₓ ) && throw( DimensionMismatch( "xₗ must have $nₓ rows" ) )
    ( size( xᵤ, 1 ) != 1 ) && ( size( xᵤ, 1 ) != nₓ ) && throw( DimensionMismatch( "xᵤ must have $nₓ rows" ) )
    ( size( uₗ, 1 ) != 1 ) && ( size( uₗ, 1 ) != nᵤ ) && throw( DimensionMismatch( "uₗ must have $nᵤ rows" ) )
    ( size( uᵤ, 1 ) != 1 ) && ( size( uᵤ, 1 ) != nᵤ ) && throw( DimensionMismatch( "uᵤ must have $nᵤ rows" ) )

    xₗ = convertboundstomissing( T, xₗ, nₓ )
    xᵤ = convertboundstomissing( T, xᵤ, nₓ )
    uₗ = convertboundstomissing( T, uₗ, nᵤ )
    uᵤ = convertboundstomissing( T, uᵤ, nᵤ )

    ConstrainedTimeInvariantLQR{T}( sys, N, Q, Qₖ, R, P, S, K, E, F, g, xₗ, xᵤ, uₗ, uᵤ )
end

# Helper function to perform the Inf -> missing conversion on the bounds
function convertboundstomissing( T, v, n )
    if isa( v, AbstractVector )
        # Convert infinite bounds to missing
        temp = Vector{Union{T, Missing}}(undef, n)
        temp[isinf.( v )]  .= missing
        temp[.!isinf.( v )] .= v[.!isinf.( v )]

        return temp
    else
        return fill( isinf( v ) ? missing : v, (n, ) )
    end
end

"""
    getsystem(clqr::ConstrainedTimeInvariantLQR; prestabilized::Bool = true)

Constructs the potentially prestabilized system for a specific constrained LQR controller.

If a prestabilizing controller `K` is given in the LQR controller, that controller will be applied
to the system forming `A - BK` when `prestabilized` is true and the controlled system will be returned. If `prestabilized`
is false, the raw system is returned.
"""
function getsystem( clqr::ConstrainedTimeInvariantLQR; prestabilized::Bool = true )
    sys = clqr.sys

    if( prestabilized )
        return StateSpace( sys.A - sys.B*clqr.K, sys.B, sys.C, sys.D, sys.Ts )
    else
        return sys
    end
end

function isprestabilized( clqr::ConstrainedTimeInvariantLQR )
    return !iszero( clqr.K )
end

################################################################################
# Pretty print the CLQR struct
################################################################################
function _string_mat_with_headers(X::Matrix)
    p = (io, m) -> Base.print_matrix(io, m)
    return replace(sprint(p, X), "\"" => " ")
end

function print_constraints( io, clqr )
    rowsE = UnitRange( axes( clqr.E, 1 ) )
    colsE = UnitRange( axes( clqr.E, 2 ) )
    alignE = Base.alignment( io, clqr.E, rowsE, colsE, typemax( Int ), typemax( Int ), 2 )

    rowsF = UnitRange( axes( clqr.F, 1 ) )
    colsF = UnitRange( axes( clqr.F, 2 ) )
    alignF = Base.alignment( io, clqr.F, rowsF, colsF, typemax( Int ), typemax( Int ), 2 )

    rowsg = UnitRange( axes( clqr.g, 1 ) )
    colsg = UnitRange( axes( clqr.g, 2 ) )
    aligng = Base.alignment( io, clqr.g, rowsg, colsg, typemax( Int ), typemax( Int ), 2 )

    numConstraints = length( rowsF )

    # Figure out if a space should be added before the first elements
    v = view( clqr.E, :, 1 )
    spaceE = ( minimum( v ) >= 0 )

    v = view( clqr.F, :, 1 )
    spaceF = ( minimum( v ) >= 0 )

    v = view( clqr.g, :, 1 )
    spaceg = ( minimum( v ) >= 0 )

    # rowsF = rowsA - otherwise constraints are invalid
    for i in rowsF
        if( i == first( rowsF ) )
            premat  = "⎡"
            postmat = " ⎤"
        elseif( i == last( rowsF ) )
            premat  = "⎣"
            postmat = " ⎦"
        else
            premat  = "⎢"
            postmat = " ⎥"
        end

        print( io,  spaceE ? premat * " " : premat )
        Base.print_matrix_row( io, clqr.E, alignE, i, colsE, "  " )
        print( io, postmat )

        if( i == ceil( numConstraints / 2 ) )
            print( io, " x + " )
        else
            print( io, "     " )
        end

        print( io,  spaceF ? premat * " " : premat )
        Base.print_matrix_row( io, clqr.F, alignF, i, colsF, "  " )
        print( io, postmat )

        if( i == ceil( numConstraints / 2 ) )
            print( io, " u ≤ " )
        else
            print( io, "     " )
        end

        print( io,  spaceg ? premat * " " : premat )
        Base.print_matrix_row( io, clqr.g, aligng, i, colsg, "  " )
        print( io, postmat )

        println( io )
    end
end

function convertboundsfrommissing( T, v, n, neg )
    temp = Vector{Union{T, Float64}}( undef, n )
    temp[ismissing.( v )]  .= ( neg ? -Inf : Inf )
    temp[.!ismissing.( v )] .= v[.!ismissing.( v )]

    return temp
end

function print_bounds( io, clqr )
    nₓ = nstates( getsystem( clqr ) )
    nᵤ = ninputs( getsystem( clqr ) )
    T  = eltype( clqr.Q )

    xᵤ = convertboundsfrommissing( T, clqr.xᵤ, nₓ, false )
    xₗ = convertboundsfrommissing( T, clqr.xₗ, nₓ, true )
    uᵤ = convertboundsfrommissing( T, clqr.uᵤ, nᵤ, false )
    uₗ = convertboundsfrommissing( T, clqr.uₗ, nᵤ, true )

    rowsxᵤ  = UnitRange( axes( xᵤ, 1 ) )
    colsxᵤ  = UnitRange( axes( xᵤ, 2 ) )
    alignxᵤ = Base.alignment( io, xᵤ, rowsxᵤ, colsxᵤ, typemax( Int ), typemax( Int ), 2 )

    rowsxₗ  = UnitRange( axes( xₗ, 1 ) )
    colsxₗ  = UnitRange( axes( xₗ, 2 ) )
    alignxₗ = Base.alignment( io, xₗ, rowsxₗ, colsxₗ, typemax( Int ), typemax( Int ), 2 )

    rowsuᵤ  = UnitRange( axes( uᵤ, 1 ) )
    colsuᵤ  = UnitRange( axes( uᵤ, 2 ) )
    alignuᵤ = Base.alignment( io, uᵤ, rowsuᵤ, colsuᵤ, typemax( Int ), typemax( Int ), 2 )

    rowsuₗ  = UnitRange( axes( uₗ, 1 ) )
    colsuₗ  = UnitRange( axes( uₗ, 2 ) )
    alignuₗ = Base.alignment( io, uₗ, rowsuₗ, colsuₗ, typemax( Int ), typemax( Int ), 2 )

    # Figure out if a space should be added before the first elements
    spacexᵤ = any( xᵤ .< 0 )
    spacexₗ = any( xₗ .< 0 )
    spaceuᵤ = any( uᵤ .< 0 )
    spaceuₗ = any( uₗ .< 0 )

    anyinfxᵤ = any( isinf.( xᵤ ) )
    anyinfxₗ = any( isinf.( xₗ ) )
    anyinfuᵤ = any( isinf.( uᵤ ) )
    anyinfuₗ = any( isinf.( uₗ ) )

    havex = !all( ismissing.( clqr.xᵤ ) ) || !all( ismissing.( clqr.xₗ ) )
    haveu = !all( ismissing.( clqr.uᵤ ) ) || !all( ismissing.( clqr.uₗ ) )

    printx  = true
    printxₙ = true
    printu  = true
    printuₙ = true

    for i in max( rowsxᵤ, rowsuᵤ )
        if havex
            printx = printxₙ

            if ( i == first( rowsxᵤ ) ) && ( i == last( rowsxᵤ ) )
                printxₙ = false
                prematx  = "["
                postmatx = " ]"
            elseif i == first( rowsxᵤ )
                printxₙ = true
                prematx  = "⎡"
                postmatx = " ⎤"
            elseif  i == last( rowsxᵤ )
                printxₙ = false
                prematx  = "⎣"
                postmatx = " ⎦"
            elseif printx == true
                prematx  = "⎢"
                postmatx = " ⎥"
            else
                prematx  = " "
                postmatx = "  "
            end

            if printx == true
                print( io, anyinfxₗ ? ( isinf( xₗ[i] ) ? prematx * " " : prematx ) : prematx * " " )
                Base.print_matrix_row( io, xₗ, alignxₗ, i, colsxₗ, "  " )
                print( io, postmatx )

                if( i == ceil(  nₓ/ 2 ) )
                    print( io, " ≤ x ≤ " )
                else
                    print( io, "       " )
                end

                if spacexᵤ
                    # This case happens if there are any negative numbers in the array
                    # In that case we need to add 2 spaces before any infinities, and only
                    # 1 space in front of the negative numbers
                    print( io, isinf( xᵤ[i] ) ? prematx * "  " : prematx * " " )
                elseif anyinfxᵤ
                    # This case happens if there are infinities and only positive numbers.
                    # In that case we only add the space before the infinities.
                    print( io, isinf( xᵤ[i] ) ? prematx * " " : prematx )
                else
                    # In the case of all positive numbers, always add a space before every number
                    print( io, prematx * " " )
                end
                Base.print_matrix_row( io, xᵤ, alignxᵤ, i, colsxᵤ, "  " )
                print( io, postmatx )
            end
        end

        if havex && haveu
            print( io, "   " )
        end

        if haveu
            printu = printuₙ

            if ( i == first( rowsuᵤ ) ) && ( i == last( rowsuᵤ ) )
                printuₙ = false
                prematu  = "["
                postmatu = " ]"
            elseif i == first( rowsuᵤ )
                printuₙ = true
                prematu  = "⎡"
                postmatu = " ⎤"
            elseif i == last( rowsuᵤ )
                printuₙ = false
                prematu  = "⎣"
                postmatu = " ⎦"
            elseif printu == true
                prematu  = "⎢"
                postmatu = " ⎥"
            else
                prematu  = " "
                postmatu = "  "
            end

            if printu == true
                print( io, anyinfuₗ ? ( isinf( uₗ[i] ) ? prematu * " " : prematu ) : prematu * " " )
                Base.print_matrix_row( io, uₗ, alignuₗ, i, colsuₗ, "  " )
                print( io, postmatu )

                if( i == ceil(  nᵤ/ 2 ) )
                    print( io, " ≤ u ≤ " )
                else
                    print( io, "       " )
                end

                if spaceuᵤ
                    # This case happens if there are any negative numbers in the array
                    # In that case we need to add 2 spaces before any infinities, and only
                    # 1 space in front of the negative numbers
                    print( io, isinf( uᵤ[i] ) ? prematu * "  " : prematu * " " )
                elseif anyinfuᵤ
                    # This case happens if there are infinities and only positive numbers.
                    # In that case we only add the space before the infinities.
                    print( io, isinf( uᵤ[i] ) ? prematu * " " : prematu )
                else
                    # In the case of all positive numbers, always add a space before every number
                    print( io, prematu * " " )
                end
                Base.print_matrix_row( io, uᵤ, alignuᵤ, i, colsuᵤ, "  " )
                print( io, postmatu )
            end
        end

        println( io )
    end
end

Base.print(io::IO, clqr::ConstrainedTimeInvariantLQR) = show(io, clqr)

function Base.show(io::IO, clqr::ConstrainedTimeInvariantLQR)
    println( io, "Cost function:" )
    println( io, "Q = \n", _string_mat_with_headers( clqr.Q ) )
    println( io, "R = \n", _string_mat_with_headers( clqr.R ) )
    println( io, "P = \n", _string_mat_with_headers( clqr.P ) )

    if !iszero( clqr.S )
        println( io, "S = \n", _string_mat_with_headers( clqr.S ) )
    end
    println( io )

    println( io, "Predicted system:" )
    print( io, clqr.sys )
    println( io )
    println( io )

    if isprestabilized( clqr )
        println( io, "Prestabilizing controller:" )
        println( io, "K = \n", _string_mat_with_headers( clqr.K ) )
    else
        println( io, "No prestabilizing controller" )
    end

    println( io )

    if !all( ismissing.( clqr.xᵤ ) ) || !all( ismissing.( clqr.xₗ ) ) || !all( ismissing.( clqr.uᵤ ) ) || !all( ismissing.( clqr.uₗ ) )
        println( io, "Bounds:" )
        print_bounds( io, clqr )
    end

    if !iszero( clqr.E ) || ! iszero( clqr.F )
        println( io, "Constraints:" )
        println( io )
        print_constraints( io, clqr )
    end
end


text:
#Parameter extraction of a parametric sinusoidal "simulation"
#-------------------------------------------------------------------------------
module CMDimData_SampleUsage

using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot
CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
WIDTH_LEGEND = 250
WIDTH_DATA = 600 #Approx
WIDTH_PLOT = WIDTH_DATA+WIDTH_LEGEND
HEIGHT_PLOT = round(Int, WIDTH_DATA/3) #Total is roughly square
HEIGHT_CANVAS = round(Int, 3*HEIGHT_PLOT*1.2) #Add extra for overhead (titles, etc)

#==Read in results of parametric "simulation"
===============================================================================#
signal = evalfile("parametric_sin_data.jl")

#==Parameter extraction
===============================================================================#
sdim1 = 4e3 #Sample dimension 1 @ 𝑓=4kHz
sdim2 = .002 #Sample dimension 2 @ A=.002

println()
@info("Information on `signal`:")
@show ndims(signal), paramlist(signal)

#Convenience: x-value resulting from paramter reduction operations:
x_paramreduce = reverse(paramlist(signal))

#Generate DataRS with value of amplitude (A) in each parameter combination:
ampvalue = parameter(signal, "A")

#Normalize signal amplitudes for ALL signals SIMULTANEOUSLY:
signal_norm = signal / ampvalue

#Compute continuous-time signal rate for ALL signals SIMULTANEOUSLY:
#   signal(𝜙, A, 𝑓, t) -> rate(𝜙, A, 𝑓):
rate = deriv(signal)

println()
@info("Reduction 0: Compute `fallx`: First fall-crossing point for ALL signals SIMULTANEOUSLY:")
#   signal(𝜙, A, 𝑓, t) -> fallx(𝜙, A, 𝑓):
fallx = xcross1(signal, xstart=0, allow=CrossType(:fall))
@show ndims(fallx), paramlist(fallx)
@show xred1 = x_paramreduce[1] #x-value for 1st fallx reduction

println()
@info("Reduction 1: Evaluate `fallx` @ $xred1=$sdim1 (=>`fallx_red1`):")
#   fallx(𝜙, A, 𝑓) -> fallx_red1(𝜙, A):
fallx_red1 = value(fallx, x=sdim1)
@show ndims(fallx_red1), paramlist(fallx_red1)
@show xred2 = x_paramreduce[2] #x-value for 2nd fallx reduction

println()
@info("Reduction 2: Evaluate `fallx_red1` @ $xred2=$sdim2 (=>`fallx_red2`):")
#   fallx_red1(𝜙, A) -> fallx_red2(𝜙):
fallx_red2 = value(fallx_red1, x=sdim2)
@show ndims(fallx_red2), paramlist(fallx_red2)
@show xred3 = x_paramreduce[3] #x-value for 3rd fallx reduction

data = (
	signal = signal, signal_norm = signal_norm,
	rate = rate, fallx = fallx,
	fallx_red1 = fallx_red1, fallx_red2 = fallx_red2,
	xred1 = xred1, xred2 = xred2,	xred3 = xred3,
	sdim1 = sdim1,	sdim2 = sdim2,
)

#==Generate plots
===============================================================================#
PB = EasyPlot.load_plotbuilders(@__DIR__,
	initial = "bld_parametric_sin_initial.jl",
	explore = "bld_parametric_sin_explore.jl",
)

function adjust_legend(gtkplot)
	mplot = gtkplot.src
	for sp in mplot.subplots
		sp.layout[:halloc_legend] = WIDTH_LEGEND
	end
end
plotdisplay = EasyPlot.GUIDisplay(:InspectDR, postproc=adjust_legend)
	opt = EasyPlot.ShowOptions(dim=set(w=WIDTH_PLOT, h=HEIGHT_CANVAS))

plotset1 = EasyPlot.build(PB[:initial], data)
	plotgui1 = display(plotdisplay, plotset1)
	EasyPlot._write(:png, "parametric_sin_initial.png", opt, plotgui1)
plotset2 = EasyPlot.build(PB[:explore], data)
	plotgui2 = display(plotdisplay, plotset2)
	EasyPlot._write(:png, "parametric_sin_explore.png", opt, plotgui2)

end #module
:SampleCode_Executed


text:
using ForwardDiff
using DelimitedFiles
using MAT

nn_width = 20
nn_depth = 3
σY = 100.0
dσ = 0.1
activation = "tanh"
exp_id = 1

if length(ARGS)==6
    nn_width = parse(Int64, ARGS[1])
    nn_depth = parse(Int64, ARGS[2])
    σY = parse(Float64, ARGS[3])
    dσ = parse(Float64, ARGS[4])
    activation = ARGS[5]
    exp_id = parse(Int64, ARGS[6])
end
FILEID = "nn_width$(nn_width)nn_depth$(nn_depth)sigmaY$(σY)dsigma$(dσ)activation$(activation)exp_id$exp_id"
FILEID2 = "nn_width$(nn_width)nn_depth$(nn_depth)sigmaY$(σY)dsigma$(dσ)activation$(activation)"

@info FILEID

config = [[nn_width for i=1:nn_depth]...,1]

stress_scale = 1e2
strain_scale = 1e-3

E0 = 200.0e3/(stress_scale/strain_scale)

function nn(ε, ε0, σ0)
    @show ε, ε0, σ0
    x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
    H = ae(x, config, "default", activation = activation)^2
    s = σ0^2
    i = sigmoid((s - σY^2)/(dσ*σY^2))

    y = ( H.* i + E0 * (1-i) ) .* (ε-ε0)*stress_scale/strain_scale + σ0
    return y

end

function sigmoid_(z)
    return 1.0 / (1.0 + exp(-z))
end

function post_nn(ε::Float64, ε0::Float64, σ0::Float64, Δt::Float64)
    f = x -> begin
        H = nndefault(reshape([x/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]^2

        # @show H
        s = σ0^2
        i = sigmoid_((s - σY^2)/(dσ*σY^2))

        (H * i + E0 * (1.0-i) )*stress_scale/strain_scale * (x-ε0) + σ0

    end
    df = ForwardDiff.derivative(f, ε)
    return f(ε), df
end


text:
const DATA_START = 632 # Start byte for a SAC file's data section
const ASCII_DATA_START = 1672 # Start character for ASCII SAC file's data section

abstract AbstractSACData
abstract AbstractTimeSeries <: AbstractSACData
abstract AbstractSpectrum <: AbstractSACData

"""
An evenly spaced time series.

Fields
======

- `hdr::SACDatHeader` - The file's header.
- `data::Vector{Float32}` - The dependent variable.

See Also
========

- `readsac`
"""
type EvenTimeSeries <: AbstractTimeSeries
    hdr::Header
    data::Vector{Float32}
end

"""
An unevenly spaced time series.

Fields
======

- `hdr::Header` - The file's header.
- `ddata::Vector{Float32}` - The dependent variable.
- `idata::Vector{Float32}` - The independent variable.

See Also
========

- `readsac`
"""
type UnevenTimeSeries <: AbstractTimeSeries
    hdr::Header
    ddata::Vector{Float32} # Dependent variable
    idata::Vector{Float32} # Independent variable
end

"""
Spectral data in complex form.

Fields
======

- `hdr::Header` - The file's header.
- `data::Vector{Complex{Float32}}` - The spectral data.

See Also
========

- `readsac`
"""
type ComplexSpectrum <: AbstractSpectrum
    hdr::Header
    data::Vector{Complex{Float32}}
end

"""
Spectral data in amplitude and phase form.

Fields
======

- `hdr::Header` - The file's header.
- `ampdata::Vector{Float32}` - The amplitude data.
- `phasedata::Vector{Float32}` - The phase data.

See Also
========

- `readsac`
"""
type AmplitudeSpectrum <: AbstractSpectrum
    hdr::Header
    ampdata::Vector{Float32}
    phasedata::Vector{Float32}
end

"""
General x vs. y data.

Fields
======

- `hdr::Header` - The file's header.
- `y::Vector{Float32}` - The y variable's data.
- `x::Vector{Float32}` - The x variable's data.

See Also
========

- `readsac`
"""
type GeneralXY <: AbstractSACData
    hdr::Header
    y::Vector{Float32}
    x::Vector{Float32}
end

# Maps types onto their Header enums.
const FILE_TYPE_ENUMS = Dict{Type,HeaderEnum}(EvenTimeSeries    => itime,
                                              UnevenTimeSeries  => itime,
                                              AmplitudeSpectrum => iamph,
                                              ComplexSpectrum   => irlim,
                                              GeneralXY         => ixy)


text:
using ImageSegmentationEvaluation
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end
@testset "Boundary extraction" begin
    img = [0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 1 1 1 1 1 0 0;
           0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0]

    expected = [0 0 0 0 0 0 0 0 0;
                0 0 1 1 1 1 1 0 0;
                0 1 1 1 1 1 1 1 0;
                0 1 1 0 0 0 1 1 0;
                0 1 1 0 0 0 1 1 0;
                0 1 1 0 0 0 1 1 0;
                0 1 1 0 0 0 1 1 0;
                0 1 1 1 1 1 1 1 0;
                0 0 1 1 1 1 1 0 0;
                0 0 0 0 0 0 0 0 0]
    @test expected == boundary_map(BoundaryGradient(), img)

    expected = [ 0  0  0  0  0  0  0  0  0;
                 0  1  1  1  1  1  1  0  0;
                 0  1  0  0  0  0  1  0  0;
                 0  1  0  0  0  0  1  0  0;
                 0  1  0  0  0  0  1  0  0;
                 0  1  0  0  0  0  1  0  0;
                 0  1  0  0  0  0  1  0  0;
                 0  1  1  1  1  1  1  0  0;
                 0  0  0  0  0  0  0  0  0;
                 0  0  0  0  0  0  0  0  0]
    @test expected == boundary_map(BoundaryShift(), img)
    display(boundary_map(BoundaryMooreTracing(), img)  )
    @test expected == boundary_map(BoundaryMooreTracing(), img)

end
@testset "Boundary Displacement Error" begin
    function distances(B1, B2)
        d=[]
        for j = 1:size(B1,1)
            p = B1[j,:]
            push!(d, minimum([norm(B2[k,:]-p) for k = 1:size(B2,1)]))
        end
        return d
    end
    img = [0 0 0 0 0;
           0 1 1 1 0;
           0 1 1 1 0;
           0 1 1 1 0;
           0 0 0 0 0]
    gt = [0 0 0 0 0;
          0 0 1 1 0;
          0 0 1 1 0;
          0 0 1 1 0;
          0 0 0 0 0]

    img_boundary =  [0 1 1 1 0;
                     1 1 1 1 1;
                     1 1 0 1 1;
                     1 1 1 1 1;
                     0 1 1 1 0]
    gt_boundary = []
    display(Integer.(boundary_map(BoundaryGradient(), gt)))
    B1 = ind2sub(size(img),find(img.==1))
    B1 = hcat(B1...)
    B2 = ind2sub(size(gt), find(gt.==1))
    B2 = hcat(B2...)

    d1 = mean(distances(B1, B2))
    d2 = mean(distances(B2, B1))
    println(d1)
    println(d2)
    @test evaluate(BoundaryDisplacementError(), img, gt) == (d1+d2)/2
end

# write your own tests here
@test 1 == 2


text:
# a simple julia script to collect all our images into arrays and be happy there
# basically we collect all our images, read them into a massive array and then save them as numpy files, and possibly view them. it shouldn't be that difficult to be honest, but it's just what we do. we've got to figure out the walkdir file so let's just figure this out to be honest, I really don't know

using Base.Filesystem
using Images
using ImageView
using NPZ
using FileIO

const extensions = ("jpg","png","jpeg")

function getCurrentDirectory()
	return pwd()
end

function iterateThroughDirs(rootdir, basename="", save::Bool=true,verbose::Bool=true)
	for dirs in walkdir(rootdir)
		if verbose==true
			println("walking through dirs")
		end
		for dir in dirs
			if verbose==true
				println("Processing dir: " * string(dir))
			end
			arr = getImgArray(dir)
			if save==true
				fname = basename * "_" * string(dir) * "_images"
				saveNpzArray(fname, arr)
			end
		end
	end
end

function getImgArray(dir)
	#init array
	img_names = readdir(dir)
	counter = 0
	img_array = 0
	for img_name in img_names
	#for the moment I will just define it wrong so it's okay but then die afterwards  idon't know
		if counter ==0
			img_array = load(img_name)
		end
		if counter >0
			cat(3, img_array, load(img_name))
		end
		counter +=1
	end
	return img_array
end

function saveNpzArray(fname, arr)
	npz.write(fname*".npy", arr)
end

max_stars_repo_name: Bmillidgework/Misc-Maths

max_issues_repo_name: Bmillidgework/Misc-Maths

max_forks_repo_name: Bmillidgework/Misc-Maths


text:
function dnsol(n, L::TronSparseMatrixCSC, r)
    # Solve L*x = r and store the result in r.

    @inbounds for j=1:n
        temp = r[j]/L.diag_vals[j]
        for k=L.colptr[j]:L.colptr[j+1]-1
            r[L.rowval[k]] = r[L.rowval[k]] - L.tril_vals[k]*temp
        end
        r[j] = temp
    end

    return
end

function dnsol(n, L::TronDenseMatrix, r)
    # Solve L*x = r and store the result in r.

    @inbounds for j=1:n
        temp = r[j]/L.vals[j,j]
        for k=j+1:n
            r[k] = r[k] - L.vals[k,j]*temp
        end
        r[j] = temp
    end
end

@inline function dnsol(n::Int, L::CuDeviceArray{Float64,2},
                       r::CuDeviceArray{Float64,1})
    # Solve L*x = r and store the result in r.

    tx = threadIdx().x
    ty = threadIdx().y

    @inbounds for j=1:n
        if tx == 1 && ty == 1
            r[j] = r[j] / L[j,j]
        end
        CUDA.sync_threads()

        if tx > j && tx <= n && ty == 1
            r[tx] = r[tx] - L[tx,j]*r[j]
        end
        CUDA.sync_threads()
    end

    return
end


text:
# Run tests on functions in src/by_reference_doc/genres.jl

@testset verbose = true "GET-request endpoints for genres" begin

    @test ~isempty(Spotify.genres_get()[1])

end


text:
export set!
for (specC, tJ, tC) in (("",    BigFloat,   Ref{BigFloat}),
                        ("_ui", CulongMax,  Culong),
                        ("_si", ClongMax,   Clong),
                        ("_d",  CdoubleMax, Cdouble),
                        ("_z",  BigInt,     Ref{BigInt}))

    @eval begin
        function set!(rop::BigFloat, op::$tJ, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_set", specC)), "libmpfr"), Int32,
                  (Ref{BigFloat}, $tC, MPFR.MPFRRoundingMode),
                  rop, op, rnd)
            return rop
        end
    end
end

@doc """
    set!(rop, op, rnd)

Set the value of `rop` from `op`, rounded toward the given direction `rnd`.
""" set!


text:
"""
    BoundingBox(xmin, ymin, zmin, xmax, ymax, zmax)
    BoundingBox(corner_min, corner_max)

Construct a three-dimensional `BoundingBox` representing a closed, axis-aligned rectangular
prism with the given coordinates.

Note that users will usually want to call the `boundingbox` function, which constructs a
`BoundingBox` guaranteed to contain all of the inputed geometry.

See also `AbstractRegion`
"""
struct BoundingBox{T <: Real} <: AbstractRegion{T}
    xmin::T
    ymin::T
    zmin::T
    xmax::T
    ymax::T
    zmax::T
end

@inline function BoundingBox(xmin, ymin, zmin, xmax, ymax, zmax)
    T = promote_type(typeof(xmin), promote_type(typeof(ymin), promote_type(typeof(zmin), promote_type(typeof(xmax), promote_type(typeof(ymax), typeof(zmax))))))
    BoundingBox(convert(T, xmin), convert(T, ymin), convert(T, zmin), convert(T, xmax), convert(T, ymax), convert(T, zmax))
end

@inline BoundingBox(v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}) = boundingbox(v1, v2)

@inline function Base.:(==)(bb1::BoundingBox, bb2::BoundingBox)
    bb1.xmin == bb2.xmin &&
        bb1.xmax == bb2.xmax &&
        bb1.ymin == bb2.ymin &&
        bb1.ymax == bb2.ymax &&
        bb1.zmin == bb2.zmin &&
        bb1.zmax == bb2.zmax
end

@inline function in(p::StaticVector{3, <:Real}, bb::BoundingBox)
    @inbounds return (p[1] >= bb.xmin) & (p[2] >= bb.ymin) & (p[3] >= bb.zmin) &
        (p[1] <= bb.xmax) & (p[2] <= bb.ymax) & (p[3] <= bb.zmax)
end

"""
    intersects(bb1, bb2)

Returns true if the bounding boxes `bb1` and `bb2` intersect.
"""
function intersects(bb1::BoundingBox, bb2::BoundingBox)
    ! ( bb1.xmax < bb2.xmin || bb1.xmin > bb2.xmax ||
        bb1.ymax < bb2.ymin || bb1.ymin > bb2.ymax ||
        bb1.zmax < bb2.zmin || bb1.zmin > bb2.zmax )
end

@inline might_intersect(bb1::BoundingBox, bb2::BoundingBox) = intersects(bb1, bb2)

"""
    isempty(bb)

Returns true if the bounding box bb is empty.
"""
function isempty(bb::BoundingBox)
    bb.xmin > bb.xmax || bb.ymin > bb.ymax || bb.zmin > bb.zmax
end

# contains (issubset, ⊑), intersect (∩)

"""
    pad(bb::BoundingBox, dist)

Pad out the bounding box `bb` by distance `dist`.
"""
pad(bb::BoundingBox, dist) = BoundingBox(bb.xmin-dist, bb.ymin-dist, bb.zmin-dist, bb.xmax+dist, bb.ymax+dist, bb.zmax+dist)

"""
    wireframe(bb::BoundingBox)

Returns a vector of 12 `Line`s representing the wireframe of the bounding box
(primarily for plotting purposes).
"""
function wireframe(bb::BoundingBox)
    if isempty(bb)
        return Vector{Line{3, eltype(eltype(bb))}}()
    else
        @static if VERSION < v"0.7"
            lines = Vector{Line{3, eltype(eltype(bb))}}(12)
        else
            lines = Vector{Line{3, eltype(eltype(bb))}}(undef, 12)
        end
        @inbounds begin
            lines[1]  = Line(SVector(bb.xmin, bb.ymin, bb.zmin), SVector(bb.xmax, bb.ymin, bb.zmin))
            lines[2]  = Line(SVector(bb.xmin, bb.ymax, bb.zmin), SVector(bb.xmax, bb.ymax, bb.zmin))
            lines[3]  = Line(SVector(bb.xmin, bb.ymin, bb.zmax), SVector(bb.xmax, bb.ymin, bb.zmax))
            lines[4]  = Line(SVector(bb.xmin, bb.ymax, bb.zmax), SVector(bb.xmax, bb.ymax, bb.zmax))

            lines[5]  = Line(SVector(bb.xmin, bb.ymin, bb.zmin), SVector(bb.xmin, bb.ymax, bb.zmin))
            lines[6]  = Line(SVector(bb.xmax, bb.ymin, bb.zmin), SVector(bb.xmax, bb.ymax, bb.zmin))
            lines[7]  = Line(SVector(bb.xmin, bb.ymin, bb.zmax), SVector(bb.xmin, bb.ymax, bb.zmax))
            lines[8]  = Line(SVector(bb.xmax, bb.ymin, bb.zmax), SVector(bb.xmax, bb.ymax, bb.zmax))

            lines[9]  = Line(SVector(bb.xmin, bb.ymin, bb.zmin), SVector(bb.xmin, bb.ymin, bb.zmax))
            lines[10] = Line(SVector(bb.xmax, bb.ymin, bb.zmin), SVector(bb.xmax, bb.ymin, bb.zmax))
            lines[11] = Line(SVector(bb.xmin, bb.ymax, bb.zmin), SVector(bb.xmin, bb.ymax, bb.zmax))
            lines[12] = Line(SVector(bb.xmax, bb.ymax, bb.zmin), SVector(bb.xmax, bb.ymax, bb.zmax))
        end
        return lines
    end
end

"""
    boundingbox(geometry...)

Construct a `BoundingBox` which encapsulates all the `geometry` given. The geometry might
including objects like points (length-3 vectors), other bounding boxes, and other Roames
geometry types.
"""
boundingbox(bb::BoundingBox) = bb
boundingbox(x) = error("BoundingBox for type $(typeof(x)) not implemented.")
@inline boundingbox(geom1, geoms...) = boundingbox(boundingbox(geom1), boundingbox(geoms...))
@inline function boundingbox(bb::BoundingBox, geoms...)
    bb2 = boundingbox(geoms...)
    T = promote_type(eltype(eltype(bb)), eltype(eltype(bb2)))
    if isempty(bb2)
        BoundingBox(convert(T, bb.xmin), convert(T, bb.ymin), convert(T, bb.zmin),
                    convert(T, bb.xmax), convert(T, bb.ymax), convert(T, bb.zmax))
    elseif isempty(bb)
        BoundingBox(convert(T, bb2.xmin), convert(T, bb2.ymin), convert(T, bb2.zmin),
                    convert(T, bb2.xmax), convert(T, bb2.ymax), convert(T, bb2.zmax))
    else
        BoundingBox(min(bb.xmin, bb2.xmin), min(bb.ymin, bb2.ymin), min(bb.zmin, bb2.zmin),
                    max(bb.xmax, bb2.xmax), max(bb.ymax, bb2.ymax), max(bb.zmax, bb2.zmax))
    end
end

"""
    boundingbox(geometries::AbstractVector)

Construct a `BoundingBox` from `geometries`, which is a vector containing an arbitrary
number of points, lines, bounding boxes or any other Roames geometry type.
"""
boundingbox(geometries::AbstractVector) = mapreduce(boundingbox, boundingbox, geometries)

# micro-optimization for collections of points (3x faster)
function boundingbox(points::AbstractVector{<:StaticVector{2, T}}) where T <: Real
    xmin = typemax(T)
    xmax = typemin(T)
    ymin = typemax(T)
    ymax = typemin(T)

    @inbounds for p ∈ points
        x = p[1]
        y = p[2]

        if x < xmin
            xmin = x
        end
        if y < ymin
            ymin = y
        end
        if x > xmax
            xmax = x
        end
        if y > ymax
            ymax = y
        end
    end

    return BoundingBox(xmin, ymin, typemin(T), xmax, ymax, typemax(T))
end

function boundingbox(points::AbstractVector{<:StaticVector{3, T}}) where T <: Real
    xmin = typemax(T)
    xmax = typemin(T)
    ymin = typemax(T)
    ymax = typemin(T)
    zmin = typemax(T)
    zmax = typemin(T)

    @inbounds for p ∈ points
        x = p[1]
        y = p[2]
        z = p[3]

        if x < xmin
            xmin = x
        end
        if y < ymin
            ymin = y
        end
        if z < zmin
            zmin = z
        end
        if x > xmax
            xmax = x
        end
        if y > ymax
            ymax = y
        end
        if z > zmax
            zmax = z
        end
    end

    return BoundingBox(xmin, ymin, zmin, xmax, ymax, zmax)
end

# Point
function boundingbox(point::StaticVector{2, T}) where {T <: Real}
    @inbounds return BoundingBox(point[1], point[2], typemin(T), point[1], point[2], typemax(T))
end

function boundingbox(point::StaticVector{3, T}) where {T <: Real}
    @inbounds return BoundingBox(point[1], point[2], point[3], point[1], point[2], point[3])
end

# Line
boundingbox(line::Line) = boundingbox(line.p1, line.p2)

# LineString
boundingbox(ls::LineString) = boundingbox(ls.points)

# Polygon
boundingbox(p::Polygon) = boundingbox(p.exterior)

# Catenary
function boundingbox(cat::Catenary)
    if cat.lmin < 0 && cat.lmax > 0
        boundingbox(cat[cat.lmin], cat[0], cat[cat.lmax])
    else
        boundingbox(cat[cat.lmin], cat[cat.lmax])
    end
end

# Quadratic (TODO)
function boundingbox(q::Quadratic)
    error("BoundingBox for quadratic is not yet implemented")
end