text:
abstract type ObjWrapper end

"""
    ScalarObjective

"""
struct ScalarObjective{Tf, Tg, Tfg, Tfgh, Th, Thv, Tbf, P}
    f::Tf
    g::Tg
    fg::Tfg
    fgh::Tfgh
    h::Th
    hv::Thv
    batched_f::Tbf
    param::P
end
ScalarObjective(;f=nothing, g=nothing, fg=nothing, fgh=nothing, h=nothing) = ScalarObjective(f, g, fg, fgh, h, nothing, nothing, nothing)
has_param(so::ScalarObjective) = so.param === nothing ? false : true
function value(so::ScalarObjective, x)
    if has_param(so)
        return so.f(x, so.param)
    else
        return so.f(x)
    end
end
# need fall back for the case where fgh is not there
function upto_gradient(so::ScalarObjective, ∇f, x)
    if has_param(so)
        return so.fg(∇f, x, so.param)
    else
        return so.fg(∇f, x)
    end
end
function upto_hessian(so::ScalarObjective, ∇f, ∇²f, x)
    if has_param(so)
        return so.fgh(∇f, ∇²f, x, so.param)
    else
        return so.fgh(∇f, ∇²f, x)
    end
end
has_batched_f(so::ScalarObjective) = !(so.batched_f === nothing)
"""
    batched_value(obj, X)

Return the objective evaluated at all elements of X. If obj contains
a batched_f it will have X passed collectively, else f will be broadcasted
across the elements of X.
"""
function batched_value(so::ScalarObjective, X)
    if has_batched_f(so)
        if has_param(so)
            return so.batched_f(X, so.param)
        else
            return so.batched_f(X)
        end
    else
        if has_param(so)
            return so.f.(X, Ref(so.param))
        else
            return so.f.(X)
        end
    end
end
function batched_value(so::ScalarObjective, F, X)

    if has_batched_f(so)# add
        if has_param(so)
            return F = so.batched_f(F, X, so.param)
        else
            return F = so.batched_f(F, X)
        end
    else
        if has_param(so)
            F .= so.f.(X, Ref(so.param))
            return F
        else
            F .= so.f.(X)
            return F
        end
    end
end

## If prob is a NEqProblem, then we can just dispatch to least squares MeritObjective
# if fast JacVec exists then maybe even line searches that updates the gradient can be used???
struct LineObjective!{TP, T1, T2, T3}
    prob::TP
    ∇fz::T1
    z::T2
    x::T2
    d::T2
    φ0::T3
    dφ0::T3
end
function (le::LineObjective!)(λ)
    z = retract!(_manifold(le.prob), le.z, le.x, le.d, λ)
    ϕ = value(le.prob, z)
    (ϕ=ϕ, z=z)
end
function (le::LineObjective!)(λ, calc_grad::Bool)
    f, g = upto_gradient(le.prob, le.∇fz, retract!(_manifold(le.prob), le.z, le.x, le.d, λ))
    (ϕ=f, dϕ=real(dot(g, le.d))) # because complex dot might not have exactly zero im part and it's the wrong type
end
struct LineObjective{TP, T1, T2, T3}
    prob::TP
    ∇fz::T1
    z::T2
    x::T2
    d::T2
    φ0::T3
    dφ0::T3
end
function (le::LineObjective)(λ)
    z = retract(_manifold(le.prob), le.x, le.d, λ)
    _value = value(le.prob, z)
    if le.prob.objective isa MeritObjective
        return (ϕ = _value.ϕ, _value.Fx)
    end
    (ϕ = _value,)
end
function (le::LineObjective)(λ, calc_grad::Bool)
    f, g = upto_gradient(le.prob, le.∇fz, retract(_manifold(le.prob), le.x, le.d, λ))
    (ϕ=f, dϕ=real(dot(g, le.d))) # because complex dot might not have exactly zero im part and it's the wrong type
end

# We call real on dφ0 because x and df might be complex
_lineobjective(mstyle::InPlace, prob::AbstractProblem, ∇fz, z, x, d, φ0, dφ0) = LineObjective!(prob, ∇fz, z, x, d, φ0, real(dφ0))
_lineobjective(mstyle::OutOfPlace, prob::AbstractProblem, ∇fz, z, x, d, φ0, dφ0) = LineObjective(prob, ∇fz, z, x, d, φ0, real(dφ0))

struct MeritObjective{TP, T1, T2, T3, T4, T5}
  prob::TP
  F::T1
  FJ::T2
  Fx::T3
  Jx::T4
  d::T5
end
function value(mo::MeritObjective, x)
  Fx = mo.F(mo.Fx, x)
  (ϕ=(norm(Fx)^2)/2, Fx = Fx)
end

struct LsqWrapper{Tobj, TF, TJ} <: ObjWrapper
  R::Tobj
  F::TF
  J::TJ
end
function (lw::LsqWrapper)(x)
  F = lw.R(lw.F, x)
  sum(abs2, F)/2
end
function (lw::LsqWrapper)(∇f, x)
  _F, _J = lw.R(lw.F, lw.J, x)
  copyto!(∇f, sum(_J; dims=1))
  sum(abs2, _F), ∇f
end

struct LeastSquaresObjective{TFx, TJx, Tf, Tfj, Td}
    Fx::TFx
    Jx::TJx
    F::Tf
    FJ::Tfj
    ydata::Td
end
has_ydata(lsq::LeastSquaresObjective) = !(lsq.ydata === nothing)
function value(lsq::LeastSquaresObjective, x)
    # Evaluate the residual system or the "predicted" value for LeastSquaresObjective
    Fx = lsq.F(lsq.Fx, x)

    # If this comes from a LeastSquaresProblem there will be a lhs to subtract
    if has_ydata(lsq)
        Fx .= Fx .- lsq.ydata
    end
    lsq.Fx .= Fx

    # Least Squares
    f = (norm(Fx)^2)/2
    return f
end
function batched_value(lsq::LeastSquaresObjective, F, X)
   F .= value.(Ref(lsq), X)
end
function upto_gradient(lsq::LeastSquaresObjective, Fx, x)
    # Evaluate the residual system or the "predicted" value for LeastSquares
    # and the Jacobian of either one
    Fx_sq, Jx_sq = lsq.FJ(lsq.Fx, lsq.Fx*x', x)

    # If this comes from a LeastSquaresProblem there will be a lhs to subtract
    if has_ydata(lsq)
        Fx_sq .= Fx_sq .- lsq.ydata
    end
    lsq.Fx .= Fx_sq

    # Least Squares
    f = (norm(Fx_sq)^2)/2
    Fx .= Jx_sq'*Fx_sq
   return f, Fx
end
function upto_hessian(lsq::LeastSquaresObjective, Fx, Jx, x)  #Fx is the gradient and Jx is the Hessian
    Fx_sq, Jx_sq = lsq.FJ(lsq.Fx, lsq.Jx, x)

    lsq.Fx .= Fx_sq
    f = (norm(Fx)^2)/2
    # this is the gradient
    Fx .= Jx_sq'*Fx_sq
    # As you may notice, this can be expensive... Because the gradient
    # is going to be very simple. May want to create a
    # special type or way to hook into trust regions here. We can exploit
    # that we only need the cauchy and the newton steps, not any shifted
    # systems. There is no need to get the full hessian. because these two
    # steps are don't need these multiplies
    # This is the Hessian
    Jx .= Jx_sq'*Jx_sq
    return f, Fx, Jx
end


text:
module TestSets

using Test
using VecMathOptInterface
const MOI = VecMathOptInterface

include("dummy.jl")

"""
    MutLessThan{T<:Real} <: MOI.AbstractScalarSet

A mutable `LessThan`-like set to test `copy` of indicator set
"""
mutable struct MutLessThan{T<:Real} <: MOI.AbstractScalarSet
    upper::T
    MutLessThan(v::T) where {T<:Real} = new{T}(v)
end

Base.copy(mlt::MutLessThan) = MutLessThan(Base.copy(mlt.upper))

function test_sets_equals()
    # By default, `==` redirects to `===`, it works for bits type
    # but not for `BigInt`s. We define functions creating different
    # instances so that `a() !== a()`.
    a() = big(1)
    b() = big(2)
    @test a() !== a()
    @test b() !== b()
    for S in [
        MOI.LessThan,
        MOI.GreaterThan,
        MOI.EqualTo,
        MOI.PowerCone,
        MOI.DualPowerCone,
    ]
        @test S(a()) == S(a())
        @test S(a()) != S(b())
        @test S(b()) == S(b())
        @test S(b()) != S(a())
    end
    for S in [MOI.Interval, MOI.Semicontinuous, MOI.Semiinteger]
        @test S(a(), b()) == S(a(), b())
        @test S(a(), b()) != S(b(), a())
        @test S(a(), b()) != S(b(), b())
        @test S(a(), b()) != S(a(), a())
        @test S(a(), a()) != S(b(), b())
        @test S(a(), a()) == S(a(), a())
    end
    S = MOI.Indicator
    A() = MOI.LessThan(a())
    B() = MOI.LessThan(b())
    @test S{MOI.ACTIVATE_ON_ZERO}(A()) == S{MOI.ACTIVATE_ON_ZERO}(A())
    @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ONE}(A())
    @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ZERO}(B())
    @test S{MOI.ACTIVATE_ON_ONE}(A()) != S{MOI.ACTIVATE_ON_ONE}(B())
end

function test_sets_sos1_copy()
    s = MOI.SOS1([1.0])
    s_copy = copy(s)
    s_copy.weights[1] = 2.0
    @test s.weights[1] == 1.0
end

function test_sets_sos2_copy()
    s = MOI.SOS2([1.0])
    s_copy = copy(s)
    s_copy.weights[1] = 2.0
    @test s.weights[1] == 1.0
end

function test_sets_indicator_copy()
    s1 = MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(4.0))
    s2 = MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(MOI.GreaterThan(4.0))
    s1_copy = copy(s1)
    s2_copy = copy(s2)
    @test s1_copy isa MOI.Indicator{MOI.ACTIVATE_ON_ONE}
    @test s1 == s1_copy
    @test s2_copy isa MOI.Indicator{MOI.ACTIVATE_ON_ZERO}
    @test s2 == s2_copy
    s3 = MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(MutLessThan(4.0))
    s3_copy = copy(s3)
    @test s3.set.upper ≈ 4.0
    s3_copy.set.upper = 5.0
    @test s3.set.upper ≈ 4.0
    @test s3_copy.set.upper ≈ 5.0
end

function test_sets_broadcast()
    model = DummyModelWithAdd()
    x = MOI.add_variables(model, 3)
    cis = MOI.add_constraint.(model, x, MOI.EqualTo(0.0))
    @test cis isa
          Vector{MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}}
    @test length(cis) == 3
end

function test_sets_dimension()
    @test MOI.dimension(MOI.ExponentialCone()) == 3
    @test MOI.dimension(MOI.DualExponentialCone()) == 3
    @test MOI.dimension(MOI.PowerCone(1 / 2)) == 3
    @test MOI.dimension(MOI.DualPowerCone(1 / 2)) == 3
    @test MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(5)) == 15
    @test MOI.dimension(MOI.PositiveSemidefiniteConeSquare(5)) == 25
    @test MOI.dimension(MOI.LogDetConeTriangle(5)) == 17
    @test MOI.dimension(MOI.LogDetConeSquare(5)) == 27
    @test MOI.dimension(MOI.RootDetConeTriangle(5)) == 16
    @test MOI.dimension(MOI.RootDetConeSquare(5)) == 26
    @test MOI.dimension(MOI.SOS1([1.0, 2.0])) == 2
    @test MOI.dimension(MOI.SOS2([1.0, 2.0])) == 2
    @test MOI.dimension(
        MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(1.0)),
    ) == 2
    @test MOI.dimension(MOI.Complements(10)) == 10
end

function test_sets_DimensionMismatch()
    for (S, min_dimension) in (
        (MOI.Reals, 0),
        (MOI.Zeros, 0),
        (MOI.Nonnegatives, 0),
        (MOI.Nonpositives, 0),
        (MOI.NormInfinityCone, 1),
        (MOI.NormOneCone, 1),
        (MOI.SecondOrderCone, 1),
        (MOI.RotatedSecondOrderCone, 2),
        (MOI.GeometricMeanCone, 2),
        (MOI.Complements, 0),
        (MOI.RelativeEntropyCone, 1),
        (MOI.PositiveSemidefiniteConeTriangle, 0),
        (MOI.PositiveSemidefiniteConeSquare, 0),
        (MOI.LogDetConeTriangle, 0),
        (MOI.LogDetConeSquare, 0),
        (MOI.RootDetConeTriangle, 0),
        (MOI.RootDetConeSquare, 0),
    )
        @test_throws DimensionMismatch S(min_dimension - 1)
        @test S(min_dimension) isa S
    end
    @test_throws DimensionMismatch MOI.NormSpectralCone(-1, 0)
    @test_throws DimensionMismatch MOI.NormSpectralCone(0, -1)
    @test MOI.NormSpectralCone(0, 0) isa MOI.NormSpectralCone
    @test_throws DimensionMismatch MOI.NormNuclearCone(-1, 0)
    @test_throws DimensionMismatch MOI.NormNuclearCone(0, -1)
    @test MOI.NormNuclearCone(0, 0) isa MOI.NormNuclearCone
    # Other dimension checks
    @test_throws DimensionMismatch MOI.RelativeEntropyCone(2)
    @test_throws DimensionMismatch MOI.Complements(-3)
    @test_throws DimensionMismatch MOI.Complements(3)
    return
end

function _dual_set_test(set1, set2)
    @test MOI.dual_set(set1) == set2
    @test MOI.dual_set_type(typeof(set1)) == typeof(set2)
    @test MOI.dual_set(set2) == set1
    @test MOI.dual_set_type(typeof(set2)) == typeof(set1)
end

function _self_dual_set_test(set)
    @test MOI.dual_set(set) == set
    @test MOI.dual_set_type(typeof(set)) == typeof(set)
end

function test_sets_dual_nonpositives()
    nonpositives3 = MOI.Nonpositives(3)
    nonpositives4 = MOI.Nonpositives(4)
    _self_dual_set_test(nonpositives3)
    @test MOI.dual_set(nonpositives3) != nonpositives4
    _self_dual_set_test(nonpositives4)
    return
end

function test_sets_dual_nonnegatives()
    nonnegatives3 = MOI.Nonnegatives(3)
    nonnegatives4 = MOI.Nonnegatives(4)
    _self_dual_set_test(nonnegatives3)
    @test MOI.dual_set(nonnegatives3) != nonnegatives4
    _self_dual_set_test(nonnegatives4)
    return
end

function test_sets_dual_zeroreal()
    zeros3 = MOI.Zeros(3)
    zeros4 = MOI.Zeros(4)
    reals3 = MOI.Reals(3)
    reals4 = MOI.Reals(4)
    _dual_set_test(zeros3, reals3)
    @test MOI.dual_set(reals3) != zeros4
    _dual_set_test(zeros4, reals4)
    @test MOI.dual_set(zeros4) != reals3
    return
end

function test_sets_dual_norm()
    norminf2 = MOI.NormInfinityCone(2)
    norminf3 = MOI.NormInfinityCone(3)
    normone2 = MOI.NormOneCone(2)
    normone3 = MOI.NormOneCone(3)
    _dual_set_test(norminf2, normone2)
    _dual_set_test(norminf3, normone3)
    @test MOI.dual_set(norminf2) != normone3
    @test MOI.dual_set(normone2) != norminf3
    return
end

function test_sets_dual_soc()
    soc2 = MOI.SecondOrderCone(2)
    soc3 = MOI.SecondOrderCone(3)
    _self_dual_set_test(soc2)
    @test MOI.dual_set(soc2) != soc3
    _self_dual_set_test(soc3)
    return
end

function test_sets_dual_rsoc()
    rsoc2 = MOI.RotatedSecondOrderCone(2)
    rsoc3 = MOI.RotatedSecondOrderCone(3)
    _self_dual_set_test(rsoc2)
    @test MOI.dual_set(rsoc2) != rsoc3
    _self_dual_set_test(rsoc3)
    return
end

function test_sets_dual_normspectral()
    normspec22 = MOI.NormSpectralCone(2, 2)
    normspec23 = MOI.NormSpectralCone(2, 3)
    normnuc22 = MOI.NormNuclearCone(2, 2)
    normnuc23 = MOI.NormNuclearCone(2, 3)
    _dual_set_test(normspec23, normnuc23)
    _dual_set_test(normspec22, normnuc22)
    @test MOI.dual_set(normspec22) != normnuc23
    @test MOI.dual_set(normnuc22) != normspec23
    return
end

function test_sets_dual_psdtriangle()
    psd2 = MOI.PositiveSemidefiniteConeTriangle(2)
    psd3 = MOI.PositiveSemidefiniteConeTriangle(3)
    _self_dual_set_test(psd2)
    @test MOI.dual_set(psd2) != psd3
    _self_dual_set_test(psd3)
    return
end

function test_sets_dual_exponential()
    exp = MOI.ExponentialCone()
    dual_exp = MOI.DualExponentialCone()
    _dual_set_test(exp, dual_exp)
    @test MOI.dual_set(exp) != exp
    _dual_set_test(dual_exp, exp)
    @test MOI.dual_set(dual_exp) != dual_exp
    return
end

function test_sets_dual_power()
    pow03 = MOI.PowerCone(0.3)
    pow04 = MOI.PowerCone(0.4)
    dual_pow03 = MOI.DualPowerCone(0.3)
    _dual_set_test(pow03, dual_pow03)
    @test MOI.dual_set(pow03) != pow03
    _dual_set_test(dual_pow03, pow03)
    @test MOI.dual_set(dual_pow03) != pow04
    @test MOI.dual_set(dual_pow03) != dual_pow03
    return
end

function test_sets_dual_psdsquare()
    s = MOI.PositiveSemidefiniteConeSquare(4)
    err = ErrorException(
        """Dual of `PositiveSemidefiniteConeSquare` is not defined in VecMathOptInterface.
            For more details see the comments in `src/Bridges/Constraint/square.jl`.""",
    )
    @test_throws err MOI.dual_set(MOI.PositiveSemidefiniteConeSquare(4))
    @test_throws err MOI.dual_set_type(MOI.PositiveSemidefiniteConeSquare)
end

function test_sets_dual_nonimplemented()
    s = MOI.LogDetConeTriangle(4)
    err = ErrorException("Dual of $s is not implemented.")
    @test_throws err MOI.dual_set(MOI.LogDetConeTriangle(4))
    err = ErrorException("Dual type of $(typeof(s)) is not implemented.")
    @test_throws err MOI.dual_set_type(MOI.LogDetConeTriangle)
end

function test_sets_Interval()
    @test MOI.Interval(MOI.GreaterThan(1.0)) === MOI.Interval(1.0, Inf)
    @test MOI.Interval(MOI.LessThan(2.0)) === MOI.Interval(-Inf, 2.0)
    @test MOI.Interval(MOI.EqualTo(3.0)) === MOI.Interval(3.0, 3.0)

    @test MOI.Interval(MOI.GreaterThan(1.0f0)) === MOI.Interval(1.0f0, Inf32)
    @test MOI.Interval(MOI.LessThan(2.0f0)) === MOI.Interval(-Inf32, 2.0f0)
    @test MOI.Interval(MOI.EqualTo(3.0f0)) === MOI.Interval(3.0f0, 3.0f0)

    @test_throws MethodError MOI.Interval(MOI.GreaterThan(false))
    @test_throws MethodError MOI.Interval(MOI.LessThan(true))
    @test MOI.Interval(MOI.EqualTo(true)) === MOI.Interval(true, true)

    @test_throws MethodError MOI.Interval(MOI.GreaterThan(1))
    @test_throws MethodError MOI.Interval(MOI.LessThan(2))
    @test MOI.Interval(MOI.EqualTo(3)) === MOI.Interval(3, 3)

    @test MOI.Interval(MOI.Interval(1, 2)) == MOI.Interval(1, 2)
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestSets.runtests()

ext: jl


text:
push!(LOAD_PATH, "../src")

using Tree
using Base.Test

# Create a manual testing tree
root = Root("luca")
children = Dict{String, Child}()
children["a"] = Child("a", root, 1<<0)
children["b"] = Child("b", root, 1<<1)
children["c"] = Child("c", children["a"], 1<<2)
children["d"] = Child("d", children["b"], 1<<3)
children["e"] = Child("e", children["c"], 1<<4)
children["f"] = Child("f", children["d"], 1<<5)
children["g"] = Child("g", children["d"], 1<<6)
children["h"] = Child("h", children["f"], 1<<7)

@testset "Instantiation" begin
    # Test Root
    @test root.children == Child[children["a"], children["b"]]
    @test root.name == "luca"

    # Test Child
    @test children["e"].children == Child[]
    @test children["a"].children == Child[children["c"]]
    @test children["e"].name == "e"
    @test children["a"].parent == root
    @test children["c"].parent == children["a"]
    @test children["a"].length == 1
end

@testset "Parsing/Lexing Newick" begin
    newickstring = "(Bovine:0.69395,(Gibbon:0.36079,(:0.33636,(Gorilla:0.17147,(Chimp, Human:0.11927):0.08386):0.06124):0.15057):0.54939, Mouse:1.21460);"
    parsedstring = "(Bovine:0.69395,(Gibbon:0.36079,(3:0.33636,(Gorilla:0.17147,(Chimp:0.0,Human:0.11927)5:0.08386)4:0.06124)2:0.15057)1:0.54939,Mouse:1.2146)root;"
    luca = newick(newickstring)

    @test luca.name == "root"
    @test [child.name for child in luca.children] == ["Bovine", "1", "Mouse"]
    @test [child.length for child in luca.children] == [0.69395, 0.54939, 1.2146]
    @test length(descendantsof(luca.children[2])) == 8
    @test newick(luca) == parsedstring
end

@testset "Parsing/Lexing JSON" begin
    jsonstring = """{"name":"root","children":[{"length":0.69395,"name":"Bovine","parent":"root","children":[]},{"length":0.54939,"name":"1","parent":"root","children":[{"length":0.36079,"name":"Gibbon","parent":"1","children":[]},{"length":0.15057,"name":"2","parent":"1","children":[{"length":0.33636,"name":"3","parent":"2","children":[]},{"length":0.06124,"name":"4","parent":"2","children":[{"length":0.17147,"name":"Gorilla","parent":"4","children":[]},{"length":0.08386,"name":"5","parent":"4","children":[{"length":0.0,"name":"Chimp","parent":"5","children":[]},{"length":0.11927,"name":"Human","parent":"5","children":[]}]}]}]}]},{"length":1.2146,"name":"Mouse","parent":"root","children":[]}]}"""
    luca = json(jsonstring)

    @test luca.name == "root"
    @test [child.name for child in luca.children] == ["Bovine", "1", "Mouse"]
    @test [child.length for child in luca.children] == [0.69395, 0.54939, 1.2146]
    @test length(descendantsof(luca.children[2])) == 8
    @test json(luca) == jsonstring
end

@testset "namemapof" begin
    # On root
    namemap = namemapof(root)

    @test length(namemap) == 9
    @test namemap["a"] == root.children[1]

    vals = sort([i for i in values(namemap)], by=x -> x.name)[1:8]
    @test vals == sort([i for i in values(children)], by=x -> x.name)

    # On child
    namemap = namemapof(children["b"])
    @test length(namemap) == 5
    @test namemap["d"] == children["d"]

    # Temporarily add a new child with an existing name
    Child("d", children["c"], 9.2)

    namemap = namemapof(root)
    @test length(namemap) == 9

    # With unique
    @test_throws ErrorException namemap = namemapof(root, unique=true)
end

@testset "delete!" begin
    # Begin by deleting the extra child we added before
    @test length(descendantsof(root)) == 9

    delete!(children["c"].children[2])

    @test length(descendantsof(root)) == 8

    # Test that children are properly added and length preserved
    new_f_length = children["d"].length + children["f"].length
    new_g_length = children["d"].length + children["g"].length

    delete!(children["d"])

    @test children["b"].children == [children["f"], children["g"]]
    @test children["f"].length == new_f_length
    @test children["g"].length == new_g_length
end

@testset "insert" begin
    # Insert the "d" we deleted before
    newchild = insert("d", children["g"], children["g"].length - children["d"].length)

    @test children["g"].length == 1 << 6
    @test newchild.length == children["d"].length
    @test newchild.children == [children["g"]]
    children["d"] = newchild
end

@testset "move!" begin
    num_children = length(descendantsof(root))
    move!(children["f"], children["d"])
    children["f"].length -= 1 << 3

    @test length(descendantsof(root)) == num_children
    @test children["d"].children == [children["g"], children["f"]]
    @test descendantsof(children["d"]) == [children["g"], children["f"], children["h"]]
end

@testset "isleaf" begin
    @test isleaf(root) == false
    @test isleaf(children["a"]) == false
    @test isleaf(children["e"]) == true
end

@testset "isredundant" begin
    nonredundant = Root("1")
    child1 = Child("1", nonredundant, 1)
    child2 = Child("1", nonredundant, 2)
    child3 = Child("1", child2, 4)
    child4 = Child("1", child2, 8)

    @test isredundant(nonredundant) == false

    for child in descendantsof(nonredundant)
        @test isredundant(child) == false
    end

    @test isredundant(root) == true
    @test isredundant(children["d"]) == true
    @test isredundant(children["e"]) == false
end

@testset "ispolytomic" begin
    polytomic = Root("1")
    child1 = Child("1", polytomic, 1)
    child2 = Child("1", polytomic, 2)
    child3 = Child("1", polytomic, 4)
    child4 = Child("1", child1, 8)
    child5 = Child("1", child1, 16)
    child6 = Child("1", child2, 32)

    @test ispolytomic(polytomic) == true

    for child in descendantsof(polytomic)
        @test ispolytomic(child) == false
    end

    @test ispolytomic(root) == false
end

@testset "lineageof/childlineageof" begin
    @test lineageof(root) == Node[root]
    @test lineageof(children["h"]) == Node[children["h"], children["f"], children["d"], children["b"], root]

    @test_throws MethodError childlineageof(root)
    @test childlineageof(children["h"]) == Child[children["h"], children["f"], children["d"], children["b"]]
end

# Test descendantsof
@test descendantsof(root) == sort(collect(values(children)), by=x -> x.name)
@test descendantsof(children["e"]) == Child[]

# Test mrcaof
@test mrcaof([root]) == root
@test mrcaof([children["h"], children["h"]]) == children["h"]
@test mrcaof([children["a"], children["e"]]) == children["a"]
@test mrcaof([children["h"], children["g"]]) == children["d"]

# Test distance
@test distance(root, root) == 0.0
@test distance(children["b"], children["b"]) == 0.0
@test distance(root, children["b"]) == 2.0
@test distance(children["a"], children["b"]) == 3.0
@test distance(children["g"], children["h"]) == 224.0

@test distance(root, root, root) == 0.0
@test distance(children["b"], children["b"],  children["b"]) == 0.0
@test distance(root, children["b"], root) == 2.0
@test distance(children["a"], children["b"], root) == 3.0
@test distance(children["g"], children["h"], children["d"]) == 224.0

# Test simplify
root2 = deepcopy(root)
simplify!(root2)
@test [child.name for child in root2.children] == ["e", "d"]
@test length(descendantsof(root2)) == 4
@test length(descendantsof(root2.children[2])) == 2
@test distance(root2.children[1], root2.children[2].children[1]) == 191.0

# Test subtree
# Test furthestleaf - lots of methods!
# Test diameter
# Test shufflenames!

max_stars_repo_name: jakobnissen/tree

max_issues_repo_name: jakobnissen/tree

max_forks_repo_name: jakobnissen/tree


text:
using Empirikos

prior = Normal()

@test Empirikos._support(Normal()) == Interval(nothing, nothing)

lfsr_0 = Empirikos.PosteriorProbability(StandardNormalSample(0.0), Interval(0,nothing))

@test lfsr_0(-1.0) == 0.0
@test lfsr_0(0.1) == 1.0

lfsr_numerator = numerator(lfsr_0)
lfsr_denominator = denominator(lfsr_0)

@test Empirikos._support(lfsr_numerator) == Interval(0, nothing)

_num = lfsr_numerator(prior)
_denom = lfsr_denominator(prior)
@test _num / _denom == 0.5

@test Empirikos.default_target_computation(lfsr_numerator, location(lfsr_numerator), prior) == Empirikos.NumeratorOfConjugate()
@test Empirikos.compute_target(Empirikos.NumeratorOfConjugate(), lfsr_numerator, prior) == _num

@test Empirikos.compute_target(Empirikos.QuadgkQuadrature(), lfsr_numerator, prior) ≈ _num
@test Empirikos.compute_target(Empirikos.QuadgkQuadrature(), lfsr_denominator, prior) ≈ _denom

@test numerator(lfsr_0)(Uniform(-1.0,0.0)) == 0.0

postmean_target = PosteriorMean(StandardNormalSample(2.0))
postmean_target_binomial = PosteriorMean(BinomialSample(3,10))

@test postmean_target(MixtureModel([prior; prior], [0.3;0.7])) ≈ postmean_target(prior)

@test postmean_target(2.1) == 2.1
@test postmean_target_binomial(0.6) == 0.6

num_target_gaussian = Empirikos.PosteriorTargetNumerator(postmean_target)
num_target_gaussian(1.0)

num_target_binomial = Empirikos.PosteriorTargetNumerator(postmean_target_binomial)
num_target_binomial(Beta(2,1))
num_target_binomial(0.4)

#Empirikos.PosteriorTargetNumerator{PosteriorMean{BinomialSample{Int64,Int64}}})

@test Empirikos.default_target_computation(num_target_gaussian, location(num_target_gaussian), prior) == Empirikos.NumeratorOfConjugate()
@test Empirikos.compute_target(Empirikos.NumeratorOfConjugate(), num_target_gaussian, prior) == num_target_gaussian(prior)

Empirikos.compute_target(Empirikos.QuadgkQuadrature(), lfsr_numerator, prior)

Z = location(num_target_gaussian)

@test quadgk(μ-> likelihood(Z, μ)*μ*pdf(prior,μ), -Inf, +Inf)[1] ≈  num_target_gaussian(prior)
@test quadgk(μ-> likelihood(Z, μ)*μ*pdf(prior,μ), -Inf, +Inf)[1] ==
     Empirikos.compute_target(Empirikos.QuadgkQuadrature(), numerator(postmean_target), prior)

_f_num_postmean_gaussian_quad(μ) = Empirikos.compute_target(Empirikos.QuadgkQuadrature(), numerator(postmean_target), Normal(μ))
@test ForwardDiff.derivative(_f_num_postmean_gaussian_quad, 0.0) ≈ central_fdm(5,1)(_f_num_postmean_gaussian_quad, 0.0)
@test ForwardDiff.derivative(_f_num_postmean_gaussian_quad, -2.0) ≈ central_fdm(5,1)(_f_num_postmean_gaussian_quad, -2.0)

_f_postmean_gaussian(μ) = postmean_target(Normal(μ))
@test ForwardDiff.derivative(_f_postmean_gaussian, 0.0) ≈ central_fdm(5,1)(_f_postmean_gaussian, 0.0)
@test ForwardDiff.derivative(_f_postmean_gaussian, 2.0) ≈ central_fdm(5,1)(_f_postmean_gaussian, 2.0)

_f_lfsr(μ) = lfsr_0(Normal(μ))
@test ForwardDiff.derivative(_f_lfsr, 0.0) ≈ central_fdm(5,1)(_f_lfsr, 0.0)
@test ForwardDiff.derivative(_f_lfsr, 2.0) ≈ central_fdm(5,1)(_f_lfsr, 2.0)


text:
using ArgParse
using DataFrames
using Statistics

import Colors
import CSV
import Pkg.TOML

function parse_commandline(args::Union{Nothing, Array{String, 1}}=nothing) # TODO: add verbosity level
    s = ArgParseSettings(prog="baysor run")
    @add_arg_table! s begin
        "--config", "-c"
            help = "TOML file with config"
        "--x-column", "-x"
            help = "Name of x column. Overrides the config value."
        "--y-column", "-y"
            help = "Name of y column. Overrides the config value."
        "--z-column", "-z"
            help = "Name of z column. Overrides the config value."
        "--gene-column", "-g"
            help = "Name of gene column. Overrides the config value."

        "--iters", "-i"
            help = "Number of iterations"
            arg_type = Int
            default = 500
        "--min-molecules-per-cell", "-m"
            help = "Minimal number of molecules for a cell to be considered as real. It's an important parameter, as it's used to infer several other parameters. Overrides the config value."
            arg_type = Int
        "--n-clusters"
            help = "Number of molecule clusters, i.e. major cell types. Depends on protocol resolution, but should not be too high. In most cases something between 3 and 15 should work well."
            arg_type = Int
            default=4
        "--num-cells-init"
            help = "Initial number of cells."
            arg_type = Int
        "--output", "-o"
            help = "Name of the output file or path to the output directory"
            default = "segmentation.csv"
        "--plot", "-p"
            help = "Save pdf with plot of the segmentation"
            action = :store_true
        "--save-polygons"
            help = "Save estimated cell boundary polygons to a file with a specified FORMAT. Only 'GeoJSON' format is currently supported. The option requires setting '-p' to work."
            range_tester = (x -> in(lowercase(x), ["geojson"]))
            metavar = "FORMAT"
        "--force-2d"
            help = "Ignores z-column in the data if it is provided"
            action = :store_true
        "--exclude-genes"
            help = "Comma-separated list of genes or regular expressions to ignore during segmentation. Example: --exclude-genes='Blank*,MALAT1'"
        "--nuclei-genes"
            help = "Comma-separated list of nuclei-specific genes. If provided, `cyto-genes` has to be set, as well."
        "--cyto-genes"
            help = "Comma-separated list of cytoplasm-specific genes. If provided, `nuclei-genes` has to be set, as well."
        "--scale-std"
            help = "Standard deviation of scale across cells. Can be either number, which means absolute value of the std, or string ended with '%' to set it relative to scale. Default: 25%"
            default = "25%"

        "--scale", "-s"
            help = "Scale parameter, which suggest approximate cell radius for the algorithm. Must be in the same units as 'x' and 'y' molecule coordinates. Overrides the config value. Sets 'estimate-scale-from-centers' to false."
            arg_type = Float64
        "--prior-segmentation-confidence"
            help = "Confidence of the `prior_segmentation` results. Value in [0; 1]. If you want the final segmentation not contradicting to prior_segmentation, set it to 1. Otherwise, if you assume errors in prior_segmentation, values in [0.2-0.7] allow flexibility for the algorithm."
            arg_type = Float64
            default = 0.2

        "coordinates"
            help = "CSV file with coordinates of molecules and gene type"
            required = true
        "prior_segmentation"
            help = "Image or a MAT file with segmentation mask (either boolean or component indexing) or CSV column with integer segmentation labels. If it's the column name, it should be preceded ':' symbol (e.g. :cell)"
    end

    return (args === nothing) ? parse_args(s) : parse_args(args, s)
end

function parse_configs(args::Union{Nothing, Array{String, 1}}=nothing)
    r = parse_commandline(args)
    if r === nothing
        return nothing
    end

    if r["scale"] !== nothing
        r["estimate-scale-from-centers"] = false
    end

    if r["config"] !== nothing
        extend_params_with_config!(r, parse_toml_config(r["config"]))
    else
        @warn "No config file provided. Using default parameters."
        extend_params_with_config!(r, get_default_config())
    end

    for k in ["gene-column", "x-column", "y-column"]
        if !(k in keys(r)) || (r[k] === nothing) # should never be the case as we have defaults
            @warn "$k must be specified"
            return nothing
        end
        r[k] = Symbol(r[k])
    end

    r["z-column"] = Symbol(r["z-column"])

    if r["prior_segmentation"] === nothing && r["scale"] === nothing
        @warn "Either `prior_segmentation` or `scale` must be provided."
        return nothing
    end

    if r["min-molecules-per-segment"] === nothing
        r["min-molecules-per-segment"] = default_param_value(:min_molecules_per_segment, r["min-molecules-per-cell"])
    end

    if isdir(r["output"]) || isdirpath(r["output"])
        r["output"] = joinpath(r["output"], "segmentation.csv")
    end

    if (r["save-polygons"] !== nothing) && (!r["plot"])
        @warn "--plot option is required for saving polygons (--save-polygons). The polygons will not be saved."
    end

    if xor(r["nuclei-genes"] === nothing, r["cyto-genes"] === nothing)
        @warn "Only one of `nuclei-genes` and `cyto-genes` is provided. It has to be either both or none."
        return nothing
    end

    if (r["nuclei-genes"] !== nothing) && (r["n-clusters"] > 1)
        @warn "Setting n-clusters > 1 is not recommended with compartment-specific expression patterns (nuclei- and cyto-genes parameters)."
    end

    return r
end

function plot_diagnostics_panel(df_res::DataFrame, assignment::Array{Int, 1}, tracer::Dict, args::Dict;
        clust_res::Union{NamedTuple, Nothing}=nothing, comp_segs::Union{NamedTuple, Nothing}=nothing)
    @info "Plot diagnostics"
    open(append_suffix(args["output"], "diagnostics.html"), "w") do io
        vega_plots = Dict{String, VegaLite.VLSpec}()
        println(io, "")
        println(io, vega_header("Report"))
        println(io, vega_style())
        println(io, "")
        # Molecule clustering convergence
        if clust_res !== nothing
            println(io, "

")
            vega_plots["vg_clust_conv"] = plot_clustering_convergence(clust_res, "Molecule clustering convergence")
        end

        if comp_segs !== nothing
            println(io, "

")
            vega_plots["vg_compart_conv"] = plot_clustering_convergence(comp_segs, "Compartment segmentation convergence")
        end

        println(io, "

")
        # Main algorithm convergence
        if (:n_components in keys(tracer)) && length(tracer[:n_components]) != 0
            p_conv = plot_num_of_cells_per_iterarion(tracer);
            show(io, MIME("text/html"), p_conv)
        end

        println(io, "

")

        # Confidence per molecule
        if :confidence in propertynames(df_res)
            println(io, "

")
            vega_plots["vg_mol_confidence"] = plot_confidence_distribution(df_res.confidence, assignment, size=(500, 250))
        end

        # Assignment confidence
        if :assignment_confidence in propertynames(df_res)
            println(io, "

")
            vega_plots["vg_assignment_confidence"] = estimate_hist(df_res.assignment_confidence[assignment .> 0]; nbins=30) |>
                @vlplot() +
                @vlplot(:bar, x={:s, title="Assignment confidence"}, x2=:e, y={:h, title="Num. molecules"},
                    title="Assignment confidence", width=500, height=250) +
                @vlplot(:rule, x={datum=0.95})
        end

        println(io, "

")

        # Num. of molecules per cell
        n_mols_per_cell = count_array(assignment, drop_zero=true)
        n_mols_per_cell = n_mols_per_cell[(n_mols_per_cell .> 1) .& (n_mols_per_cell .< quantile(n_mols_per_cell, 0.99) / 0.99)]

        println(io, "

")
        vega_plots["vg_num_transcripts"] = estimate_hist(n_mols_per_cell, nbins=50, center=false) |>
        @vlplot(:bar, x={:s, scale={domain=[1, maximum(n_mols_per_cell) + 1]}, title="Num. molecules per cell"}, x2=:e, y={:h, title="Num. cells"},
            width=500, height=250, title="Num. molecules per cell")

        println(io, "")
        if length(vega_plots) > 0
            println(io, vega_plot_html(vega_plots))
        end
        println(io, "")
    end
end

function plot_transcript_assignment_panel(df_res::DataFrame, assignment::Vector{Int}, args::Dict;
        clusters::Vector{Int}, prior_polygons::Array{Matrix{Float64}, 1}, gene_colors::Union{Vector{<:Colors.ColorTypes.Color}, Nothing}=nothing)
    if gene_colors === nothing
        @info "Estimating local colors"
        gene_colors = gene_composition_colors(df_res, args["gene-composition-neigborhood"])
    end

    @info "Estimating boundary polygons" # For some parameters, this step takes a lot of time and memory
    grid_step = args["scale"] / args["min-pixels-per-cell"];
    polygons = poly_joint = boundary_polygons(df_res, assignment; grid_step=grid_step, bandwidth=args["scale"]/10);

    if (args["save-polygons"] !== nothing) && (:z in propertynames(df_res))
        if length(unique(df_res.z)) > (size(df_res, 1) / args["min-molecules-per-cell"])
            @warn "To many values of z. Using 2D polygons"
        else
            z_vals = sort(unique(df_res.z))
            mask_per_z = [(df_res.z .≈ z) for z in z_vals]
            poly_per_z = progress_map(mask -> boundary_polygons(df_res[mask,:], assignment[mask]; grid_step=grid_step, bandwidth=args["scale"]/10), mask_per_z);
            poly_per_z = Dict("$k" => p for (k,p) in zip(z_vals, poly_per_z))
            poly_per_z["joint"] = polygons
            polygons = poly_per_z
        end
    end

    if gene_colors !== nothing
        gene_colors = Colors.alphacolor.(gene_colors, 0.5)
    end

    @info "Plot transcript assignment"
    gc_plot = plot_dataset_colors(df_res, gene_colors; polygons=poly_joint, prior_polygons=prior_polygons, min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Local expression similarity")

    clust_plot = nothing
    if !isempty(clusters)
        clust_plot = plot_dataset_colors(df_res, gene_colors; polygons=poly_joint, prior_polygons=prior_polygons, annotation=clusters, min_molecules_per_cell=args["min-molecules-per-cell"],
            min_pixels_per_cell=args["min-pixels-per-cell"], title="Molecule clustering")
    end

    open(append_suffix(args["output"], "borders.html"), "w") do io
        show(io, MIME("text/html"), gc_plot)

        if clust_plot !== nothing
            show(io, MIME("text/html"), clust_plot)
        end
    end

    return polygons
end

function parse_prior_assignment(df_spatial::DataFrame, args::Dict{String, Any})
    prior_col = Symbol(args["prior_segmentation"][2:end])
    prior_segmentation = df_spatial[!, prior_col]
    try
        prior_segmentation = Int.(prior_segmentation);
    catch
        error("The prior segmentation column '$prior_col' must be of integer type")
    end
    if minimum(prior_segmentation) < 0
        error("The prior segmentation column '$prior_col' must not contain negative numbers")
    end

    filter_segmentation_labels!(prior_segmentation, min_molecules_per_segment=args["min-molecules-per-segment"])
    if args["scale"] === nothing
        pos_data = position_data(df_spatial)
        scale, scale_std = estimate_scale_from_assignment(pos_data, prior_segmentation;
            min_mols_per_cell=args["min-molecules-per-cell"])
    else
        scale, scale_std = args["scale"], args["scale-std"]
    end

    return prior_segmentation, Matrix{Float64}[], scale, scale_std
end

function load_prior_segmentation(df_spatial::DataFrame, args::Dict{String, Any})
    if args["prior_segmentation"][1] == ':'
        return parse_prior_assignment(df_spatial, args)
    end

    @info "Loading segmentation mask..."

    prior_seg_labels = load_segmentation_mask(args["prior_segmentation"])
    prior_segmentation = Int.(staining_value_per_transcript(df_spatial, prior_seg_labels));
    filter_segmentation_labels!(prior_seg_labels, prior_segmentation; min_molecules_per_segment=args["min-molecules-per-segment"])
    GC.gc()

    scale, scale_std = args["scale"], args["scale-std"]
    if args["estimate-scale-from-centers"]
        scale, scale_std = estimate_scale_from_centers(prior_seg_labels)
    end
    @info "Done"

    prior_polygons = nothing
    if args["plot"]
        @info "Estimating prior segmentation polygons..."
        prior_polygons = extract_polygons_from_label_grid(Matrix{UInt32}(prior_seg_labels[1:5:end, 1:5:end]); grid_step=5.0) # subset to save memory and time
        @info "Done"
    end

    return prior_segmentation, prior_polygons, scale, scale_std
end

function dump_parameters(args::Dict{String, Any}, args_str::String)
    if args["config"] !== nothing
        dump_dst = abspath(append_suffix(args["output"], "config.toml"))
        if abspath(args["config"]) != dump_dst
            cp(args["config"], dump_dst, force=true)
        end
    end

    open(append_suffix(args["output"], "params.dump"), "w") do f
        println(f, "# CLI params: `$args_str`")
        TOML.print(f, Dict(k => (v !== nothing) ? ((typeof(v) === Symbol) ? String(v) : v) : "" for (k,v) in args))
    end
end

function load_and_preprocess_data!(args::Dict{String, Any})
    @info get_baysor_run_str()
    @info "Loading data..."
    df_spatial, gene_names = load_df(args, filter_cols=false)
    df_spatial[!, :molecule_id] = 1:size(df_spatial, 1)

    @info "Loaded $(size(df_spatial, 1)) transcripts"
    if size(df_spatial, 1) != size(unique(df_spatial), 1)
        @warn "$(size(df_spatial, 1) - size(unique(df_spatial), 1)) records are duplicates. You may need to filter them beforehand."
    end

    if args["gene-composition-neigborhood"] === nothing
        args["gene-composition-neigborhood"] = default_param_value(:composition_neighborhood, args["min-molecules-per-cell"], n_genes=length(gene_names))
    end

    prior_polygons = Matrix{Float64}[]
    if args["prior_segmentation"] !== nothing
        df_spatial[!, :prior_segmentation], prior_polygons, args["scale"], args["scale-std"] = load_prior_segmentation(df_spatial, args)
    end
    GC.gc()

    confidence_nn_id = default_param_value(:confidence_nn_id, args["min-molecules-per-cell"])

    @info "Estimating noise level"
    prior_seg = (args["prior_segmentation"]===nothing) ? nothing : df_spatial.prior_segmentation
    append_confidence!(df_spatial, prior_seg, nn_id=confidence_nn_id, prior_confidence=args["prior-segmentation-confidence"])
    @info "Done"

    return df_spatial, gene_names, prior_polygons
end

function estimate_molecule_clusters(df_spatial::DataFrame, n_clusters::Int)
    @info "Clustering molecules..."
    # , adjacency_type=:both, k_adj=fmax(1, div(args["min-molecules-per-cell"], 2))
    adjacent_points, adjacent_weights = build_molecule_graph_normalized(df_spatial, :confidence, filter=false);

    mol_clusts = cluster_molecules_on_mrf(df_spatial, adjacent_points, adjacent_weights; n_clusters=n_clusters, weights_pre_adjusted=false)

    @info "Done"
    return mol_clusts
end

function estimate_molecule_compartments(df_spatial::DataFrame, gene_names::Vector{String}, args::Dict{String, Any})
    @info "Estimating compartment regions..."
    comp_genes = Dict{String, Vector{String}}()

    # Parse genes from CLI
    for k in ["nuclei-genes", "cyto-genes"]
        c_genes = String.(strip.(Base.split(args[k], ",")))
        all(length.(c_genes) .> 0) || @warn "Empty genes were provided in $k"

        missing_genes = c_genes[.!in.(c_genes, Ref(Set(gene_names)))]
        (length(missing_genes) == 0) || @warn "Genes $(join(missing_genes, ',')) are missing from the data"
        c_genes = intersect(c_genes, gene_names)
        length(c_genes) > 0 || error("No genes left in $k after filtration")
        comp_genes[k] = c_genes
    end

    # Run segmentation
    adjacent_points, adjacent_weights = build_molecule_graph_normalized(df_spatial, :confidence, filter=false);

    init_probs, is_locked = init_nuclei_cyto_compartments(position_data(df_spatial), df_spatial.gene;
        nuclei_genes=comp_genes["nuclei-genes"], cyto_genes=comp_genes["cyto-genes"], gene_names=gene_names, scale=args["scale"]);

    comp_segs = segment_molecule_compartments(init_probs, is_locked, adjacent_points, adjacent_weights, df_spatial.confidence);

    @info "Done"

    id_per_gene = Dict(g => i for (i,g) in enumerate(gene_names))
    comp_genes = [[id_per_gene[g] for g in gs] for gs in values(comp_genes)]

    return comp_segs, comp_genes
end

function save_segmentation_results(bm_data::BmmData, gene_names::Vector{String}, args::Dict{String, Any};
        mol_clusts::Union{NamedTuple, Nothing}, comp_segs::Union{NamedTuple, Nothing}, prior_polygons::Union{Vector{Matrix{Float64}}, Nothing})
    segmentated_df = get_segmentation_df(bm_data, gene_names)
    cell_stat_df = get_cell_stat_df(bm_data, segmentated_df; add_qc=true)

    @info "Estimating local colors"
    gene_colors = gene_composition_colors(bm_data.x, args["gene-composition-neigborhood"])
    segmentated_df[!, :ncv_color] = "#" .* Colors.hex.(gene_colors)

    @info "Saving results to $(args["output"])"
    CSV.write(args["output"], segmentated_df[sortperm(segmentated_df.molecule_id), :]);
    CSV.write(append_suffix(args["output"], "cell_stats.csv"), cell_stat_df);

    cm = convert_segmentation_to_counts(bm_data.x.gene, bm_data.assignment; gene_names=gene_names)
    count_str = join(names(cm), "\t") * "\n" * join([join(["$v" for v in r], '\t') for r in eachrow(cm)], '\n');
    open(append_suffix(args["output"], "counts.tsv"), "w") do f; print(f, count_str) end

    if args["plot"]
        plot_diagnostics_panel(segmentated_df, segmentated_df.cell, bm_data.tracer, args; clust_res=mol_clusts, comp_segs=comp_segs)
        polygons = plot_transcript_assignment_panel(bm_data.x, bm_data.assignment, args; clusters=bm_data.cluster_per_molecule, prior_polygons=prior_polygons,
            gene_colors=gene_colors)

        if args["save-polygons"] !== nothing
            if lowercase(args["save-polygons"]) == "geojson"
                save_polygons_to_geojson(polygons, append_suffix(args["output"], "polygons.json"))
            else
                @warn "Unknown polygon format: $(args["save-polygons"])"
            end
        end
    end
end

function run_cli_main(args::Union{Nothing, Array{String, 1}}=nothing)
    args_str = join(args, " ")
    args = parse_configs(args)
    (args !== nothing) || return 1

    dump_parameters(args, args_str)

    # Set up logger

    log_file = open(append_suffix(args["output"], "log.log"), "w")
    Base.CoreLogging.global_logger(DoubleLogger(log_file, stdout; force_flush=true))

    # Load data

    df_spatial, gene_names, prior_polygons = load_and_preprocess_data!(args)

    # Region-based segmentations

    comp_segs, comp_genes = nothing, Vector{Int}[]
    adjacent_points, adjacent_weights = build_molecule_graph(df_spatial; use_local_gene_similarities=false, adjacency_type=:triangulation)[1:2]
    if args["nuclei-genes"] !== nothing
        comp_segs, comp_genes = estimate_molecule_compartments(df_spatial, gene_names, args)
        df_spatial[!, :compartment] = ["Nuclei", "Cyto", "Unknown"][comp_segs.assignment];
        df_spatial[!, :nuclei_probs] = 1 .- comp_segs.assignment_probs[2,:];

        adjacent_points, adjacent_weights = adjust_mrf_with_compartments(adjacent_points, adjacent_weights,
            comp_segs.assignment_probs[1,:], comp_segs.assignment_probs[2,:]);
    end

    mol_clusts = nothing
    if args["n-clusters"] > 1
        mol_clusts = estimate_molecule_clusters(df_spatial, args["n-clusters"])
        df_spatial[!, :cluster] = mol_clusts.assignment;
    end

    # Cell segmentation

    bm_data = initialize_bmm_data(df_spatial; scale=args["scale"], scale_std=args["scale-std"],
            n_cells_init=args["num-cells-init"], prior_seg_confidence=args["prior-segmentation-confidence"],
            min_molecules_per_cell=args["min-molecules-per-cell"], confidence_nn_id=0,
            adjacent_points=adjacent_points, adjacent_weights=adjacent_weights, na_genes=Vector{Int}(vcat(comp_genes...)));

    @info "Using $(size(position_data(bm_data), 1))D coordinates"

    history_depth = round(Int, args["iters"] * 0.1)
    bm_data = bmm!(bm_data; n_iters=args["iters"], new_component_frac=args["new-component-fraction"], new_component_weight=args["new-component-weight"],
            min_molecules_per_cell=args["min-molecules-per-cell"], assignment_history_depth=history_depth);

    @info "Processing complete."

    # Save results
    save_segmentation_results(bm_data, gene_names, args; mol_clusts=mol_clusts, comp_segs=comp_segs, prior_polygons=prior_polygons)
    @info "All done!"

    close(log_file)

    return 0
end


text:
test_file = joinpath(@__DIR__, "CH4_631g.h5")

atoms = Molecules.parse_string("""
C   -2.131551124300    2.286168823700    0.000000000000
H   -1.061551124300    2.286168823700    0.000000000000
H   -2.488213906200    1.408104616400    0.496683911300
H   -2.488218762100    2.295059432700   -1.008766153900
H   -2.488220057000    3.155340844300    0.512081313000""")

bs = BasisSet("6-31g", atoms)
bs2 = BasisSet("sto-3g", atoms)

@testset "One-Electron Integrals" begin

    @test overlap(bs) ≈ h5read(test_file, "overlap")
    @test kinetic(bs) ≈ h5read(test_file, "kinetic")
    @test nuclear(bs) ≈ h5read(test_file, "nuclear")

    @test overlap(bs, bs2) ≈ h5read(test_file, "overlap_sto3g")
    @test kinetic(bs, bs2) ≈ h5read(test_file, "kinetic_sto3g")
    @test nuclear(bs, bs2) ./2 ≈ h5read(test_file, "nuclear_sto3g")
end

@testset "Two-Electron Four-Center" begin

    @test ERI_2e4c(bs) ≈ h5read(test_file, "denseERI")

    idx, data = sparseERI_2e4c(bs)

    uniform = zeros(Int16, length(idx)*4)
    for i in eachindex(idx)
        uniform[(1+4*(i-1)):(4+4*(i-1))] .= idx[i]
    end

    @test uniform ≈ h5read(test_file, "sparseERIidx")
    @test data ≈ h5read(test_file, "sparseERIdata")
end

@testset "Two-Electron Three-Center" begin
    @test ERI_2e3c(bs, bs2) ≈ h5read(test_file, "Pqp_aux_sto3g")
end

@testset "Two-Electron Two-Center" begin
    @test ERI_2e2c(bs) ≈ h5read(test_file, "J")
end


text:
"""
    mc_opts{T}

Storage type used for parameter options in implicit bounding routine.
* `lambda::T`: Affine weighting parameter
* `kmax::Int64`: Number of iterations run
* `style::String`: Type of contractor used
* `z_rnd::Bool`: Flag for rounding z on interval bounds
* `z_rnd_eps::T`: Amount for rounding z on interval bounds
* `z_rnd_all::Bool`: Flag for rounding z on all bounds
* `z_rnd_all_eps::T`: Amount for rounding z on all bounds
* `aff_rnd::Bool`: Flag for rounding affine relaxation on interval bounds
* `aff_rnd_eps::T`: Amount for rounding affine relaxation on interval bounds
* `aff_rnd_all::Bool`: Flag for rounding affine relaxation on all bounds
* `aff_rnd_all_eps::T`: Amount for rounding affine relaxation on all bounds
* `hhj_rnd::Bool`: Flag for rounding Jacobian on interval bounds
* `hhj_rnd_eps::T`: Amount for rounding Jacobian on interval bounds
* `hhj_rnd_all::Bool`: Flag for rounding Jacobian on all bounds
* `hhj_rnd_all_eps::T`: Amount for rounding Jacobian on all bounds
* `aff_correct_eps::T`: Affine correction tolerance
"""
type mc_opts{T<:AbstractFloat}
  lambda::T
  kmax::Int64
  style::String
  z_rnd::Bool
  z_rnd_eps::T
  z_rnd_all::Bool
  z_rnd_all_eps::T
  aff_rnd::Bool
  aff_rnd_eps::T
  aff_rnd_all::Bool
  aff_rnd_all_eps::T
  hhj_rnd::Bool
  hhj_rnd_eps::T
  hhj_rnd_all::Bool
  hhj_rnd_all_eps::T
  aff_correct_eps::T
end

"""
    mc_opts(T)

Initialization function for currently sets the weight `.lambda = 0.5`, the
contractor style `.style = KrawczykCW`, the number of iterations to `.kmax = 2`,
and the affine correction tolerance as `.aff_correct_eps = 1E-12`. Other rounding
options disabled.
"""
mc_opts(T) = mc_opts{T}(0.5,2,"KrawczykCW",
                        false,0.0,false,0.0,
                        false,0.0,false,0.0,
                        false,0.0,false,0.0,1E-12)

mc_opts() = mc_opts(Float64)
"""
    set_default!(x::mc_opts)

Restores `x` to default values. Sets the weight `.lambda = 0.5`, the
contractor style `.style = KrawczykCW`, the number of iterations to `.kmax = 2`,
and the affine correction tolerance as `.aff_correct_eps = 1E-12`. Other rounding
options disabled.
"""
function set_default!(x::mc_opts{T}) where {T<:AbstractFloat}
  x.lambda = 0.5
  x.kmax = 2
  x.style = "KrawczykCW"
  x.z_rnd = false
  x.z_rnd_eps = 0.0
  x.z_rnd_all= false
  x.z_rnd_all_eps = 0.0
  x.aff_rnd= false
  x.aff_rnd_eps = 0.0
  x.aff_rnd_all= false
  x.aff_rnd_all_eps = 0.0
  x.hhj_rnd= false
  x.hhj_rnd_eps = 0.0
  x.hhj_rnd_all= false
  x.hhj_rnd_all_eps = 0.0
  x.aff_correct_eps = 1E-12
end


text:
@testset "build" begin
    mktempdir() do dir
        files = map(1:2) do i
            without_extension = "file$i"
            file = "$(without_extension).jl"
            content = pluto_notebook_content("""
                x = begin
                    sleep(3)
                    x = 3000 + $i
                end
                """)
            path = joinpath(dir, file)
            write(path, content)
            return file
        end
        build_notebooks(BuildOptions(dir, use_distributed=false))

        html_file = joinpath(dir, "file1.html")
        @test contains(read(html_file, String), "3001")

        html_file = joinpath(dir, "file2.html")
        @test contains(read(html_file, String), "3002")
    end
end

@testset "is_pluto_file" begin
    cd(mktempdir()) do
        nb_text = """
            ### A Pluto.jl notebook ###
            # v0.14.0
            """
        write("true.jl", nb_text)
        @test PlutoStaticHTML._is_pluto_file("true.jl")

        jl_text = """
            module Foo
            end # module
            """
        write("false.jl", jl_text)
        @test !PlutoStaticHTML._is_pluto_file("false.jl")
    end
end

@testset "invalid_notebook" begin
    try
        mktempdir() do dir
            content = pluto_notebook_content("@assert false")
            file = "file.jl"
            path = joinpath(dir, file)
            write(path, content)
            build_notebooks(BuildOptions(dir), [file], use_distributed=false)
        end
        error("Test should have failed")
    catch AssertionError
        @test true # Success.
    end
end

@testset "extract_previous_output" begin
    # function extract_previous_output(html::AbstractString)::String
    html = """
        lorem
        $(PlutoStaticHTML.BEGIN_IDENTIFIER)
        ipsum
        $(PlutoStaticHTML.END_IDENTIFIER)
        dolar
        """
    actual = PlutoStaticHTML.extract_previous_output(html)
    expected = """
        $(PlutoStaticHTML.BEGIN_IDENTIFIER)
        ipsum
        $(PlutoStaticHTML.END_IDENTIFIER)
        """
    @test strip(actual) == strip(expected)
end

@testset "add_documenter_css" begin
    for add_documenter_css in (true, false)
        dir = mktempdir()
        cd(dir) do
            path = joinpath(dir, "notebook.jl")
            code = pluto_notebook_content("""
                x = 1
                """)
            write(path, code)

            use_distributed = false
            output_format = documenter_output
            bo = BuildOptions(dir; use_distributed, output_format, add_documenter_css)
            build_notebooks(bo)

            output_path = joinpath(dir, "notebook.md")
            lines = readlines(output_path)
            if add_documenter_css
                @test lines[2] == "<style>"
            else
                @test lines[2] != "<style>"
            end
        end
    end
end

@testset "EditURL" begin
    dir = joinpath(pkgdir(PlutoStaticHTML), "docs", "src", "notebooks")
    bopts = BuildOptions(dir)
    in_path = "example.jl"
    kv = [
        "GITHUB_REPOSITORY" => "",
        "GITHUB_REF" => ""
    ]
    withenv(kv...) do
        @test PlutoStaticHTML._editurl_text(bopts, in_path) == ""
    end

    kv = [
        "GITHUB_REPOSITORY" => "rikhuijzer/PlutoStaticHTML.jl",
        "GITHUB_REF" => "refs/heads/main"
    ]
    withenv(kv...) do
        url = "https://github.com/rikhuijzer/PlutoStaticHTML.jl/blob/main/docs/src/notebooks/example.jl"
        @test strip(PlutoStaticHTML._editurl_text(bopts, in_path)) == strip("""
            ```@meta
            EditURL = "$url"
            ```
            """)
    end
end

@testset "Elapsed time" begin
    n = now()
    sleep(1)
    n2 = now()
    @test PlutoStaticHTML._pretty_elapsed(n2 - n) == "1 second"
    sleep(1)
    n2 = now()
    @test PlutoStaticHTML._pretty_elapsed(n2 - n) == "2 seconds"
end


text:
export iauTaiutc
"""
Time scale transformation:  International Atomic Time, TAI, to
Coordinated Universal Time, UTC.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tai1,tai2  double   TAI as a 2-part Julian Date (Note 1)

Returned:
   utc1,utc2  double   UTC as a 2-part quasi Julian Date (Notes 1-3)

Returned (function value):
               int      status: +1 = dubious year (Note 4)
                                 0 = OK
                              -1 = unacceptable date

Notes:

   1. tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned utc1
   and utc2 form an analogous pair, except that a special convention
   is used, to deal with the problem of leap seconds - see the next
   note.

   2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.  In the 1960-1972 era
   there were smaller jumps (in either direction) each time the
   linear UTC(TAI) expression was changed, and these "mini-leaps"
   are also included in the SOFA convention.

   3. The function iauD2dtf can be used to transform the UTC quasi-JD
   into calendar date and clock time, including UTC leap second
   handling.

   4. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See iauDat for further details.

Called:
   iauUtctai    UTC to TAI

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 September 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTaiutc(double tai1, double tai2, double *utc1, double *utc2)

function iauTaiutc(tai1::Real, tai2::Real)
   ref_utc1 = Ref{Float64}(0.0)
   ref_utc2 = Ref{Float64}(0.0)

   status = ccall((:iauTaiutc, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tai1), convert(Float64, tai2), ref_utc1, ref_utc2)

   return status, ref_utc1[], ref_utc2[]
end


text:
using NetCDF, LinearAlgebra
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
#pygui(:qt)
using PyPlot
using MultivariateStats
# eofs = pyimport("eofs")
using BFloat16s
using Sonums

path = "/Users/milan/Downloads/cams/ieee"
filelist = filter(x->endswith(x,"no2.grib"),readdir(path))
gribfile = xr.open_dataset(joinpath(path,filelist[1]),engine="cfgrib")
no2 = gribfile.no2.values
level = 1:size(no2)[1]

# train sonums
trainSonum16(no2)

no2lin24 = Array{Float32}(LinQuant24Array(no2))
Q = LogQuant16Array(no2);
no2log16 = Array{Float32}(LogQuant16Array(no2))
no2bf16 = Float32.(BFloat16.(no2'))
no2s16 = Float32.(Sonum16.(no2))

## error 24bit
de = abs.(log10.(no2 ./ no2lin24))

de24m = mean(de,dims=2)[:,1]
de24mi = median(de,dims=2)[:,1]
de24min = minimum(de,dims=2)[:,1]
de24max = maximum(de,dims=2)[:,1]
de24p10 = [percentile(vec(de[i,:]),10) for i in level]
de24p90 = [percentile(vec(de[i,:]),90) for i in level]

# and with 16
de = abs.(log10.(no2 ./ no2log16))

de16m = mean(de,dims=2)[:,1]
de16mi = median(de,dims=2)[:,1]
de16min = minimum(de,dims=2)[:,1]
de16max = maximum(de,dims=2)[:,1]
de16p10 = [percentile(vec(de[i,:]),10) for i in level]
de16p90 = [percentile(vec(de[i,:]),90) for i in level]

# and for bfoat16 = bitgrooming
de = abs.(log10.(no2 ./ no2bf16))

debf16m = mean(de,dims=2)[:,1]
debf16mi = median(de,dims=2)[:,1]
debf16min = minimum(de,dims=2)[:,1]
debf16max = maximum(de,dims=2)[:,1]
debf16p10 = [percentile(vec(de[i,:]),10) for i in level]
debf16p90 = [percentile(vec(de[i,:]),90) for i in level]

# and for bfoat16 = bitgrooming
de = abs.(log10.(no2 ./ no2s16))

des16m = mean(de,dims=2)[:,1]
des16mi = median(de,dims=2)[:,1]
des16min = minimum(de,dims=2)[:,1]
des16max = maximum(de,dims=2)[:,1]
des16p10 = [percentile(vec(de[i,:]),10) for i in level]
des16p90 = [percentile(vec(de[i,:]),90) for i in level]

## plot
fig,ax = subplots(1,1,figsize=(6,10))
ax.invert_yaxis()

ax.semilogx(de24m,level,"C0",lw=3,label="24-bit lin")
ax.semilogx(de16m,level,"C1",lw=3,label="16-bit log")
ax.semilogx(debf16m,level,"C2",lw=3,label="bfloat16")
ax.semilogx(des16m,level,"C3",lw=3,label="15-bit max-entropy")

ax.semilogx(de24p90,level,"C0",lw=2)
ax.semilogx(de16p90,level,"C1",lw=2)
ax.semilogx(debf16p90,level,"C2",lw=2)
ax.semilogx(des16p90,level,"C3",lw=2)

ax.semilogx(de24max,level,"C0",lw=1)
ax.semilogx(de16max,level,"C1",lw=1)
ax.semilogx(debf16max,level,"C2",lw=1)
ax.semilogx(des16max,level,"C3",lw=1)

ax.plot(0,0,"grey",lw=3,label="mean")
ax.plot(0,0,"grey",lw=2,label="90%")
ax.plot(0,0,"grey",lw=1,label="max")

ax.set_title(L"NO$_2$ compression error")
ax.set_ylabel("model level")
ax.set_xlabel("decimal error")
ax.set_ylim(137,1)
ax.legend(loc=3,ncol=2)

# a‚x.set_xlim(1e-7,1)
tight_layout()


text:
module Endpoints

using HTTP

include("json2.jl")

function parsequerystring(query::String)
    q = Dict{String,String}()
    length(query) == 0 && return q
    for field in split(query, "&")
        keyval = split(field, "=")
        length(keyval) != 2 && throw(ArgumentError("Field '$field' did not contain an '='."))
        q[unescape(keyval[1])] = unescape(keyval[2])
    end
    return q
end

# Resource Spec
  # provided as a String
  # path segments separated by "/"
  # normal URI-encoding rules apply (resource will not be encoded)
  # function positional arguments provided by "{name::type}", "{name}" or "{::Type{T}}"
    # "type" can be one of: Int, Float64, or String (default); (could potentially allow enums or Date/DateTime, bool?)
    # typed arguments will be converted to their type via `T(string_arg)`
    # if "{::Type{T}}", the path segment corresponds to an actualy Julia type `T`, which itself will be passed as an argument to the Julia function handler

# Body spec

# URI Query Parameters
 # query parameters contained in a URI such as "google.com/?search=julia"
 # will be passed to the julia function as keyword arguments
 # the julia function can explicitly define the keyword arguments it accepts (in which case, a URI that contains additional unsupported query parameters will not dispatch to this julia function)
 # the julia function can also include `f(x; kwargs...)`, which indicates a "varargs" catchall for extra keyword arguments
 # this ensures the julia function will be dispatched to, regardless of the query parameter or number of them
 # the HTTP standard doesn't restrict the use of duplicate query parameters in a URI, but julia functions can only contain unique keyword arguments,
 # duplicate keyword arguments can, however, be passed in a `kwargs...` varargs catch-all, providing a work-around
 # query parameters are always passed as String values, leaving the julia function itself to perform any conversions

macro GET(resource, func)
    generate_dispatch("GET", __module__, resource, nothing, func)
end

macro HEAD(resource, func)
    generate_dispatch("HEAD", __module__, resource, nothing, func)
end

macro POST(resource, body, func)
    generate_dispatch("POST", __module__, resource, body, func)
end

macro PUT(resource, body, func)
    generate_dispatch("PUT", __module__, resource, body, func)
end

macro DELETE(resource, func)
    generate_dispatch("DELETE", __module__, resource, nothing, func)
end

macro CONNECT(resource, func)
    generate_dispatch("CONNECT", __module__, resource, nothing, func)
end

macro OPTIONS(resource, func)
    generate_dispatch("OPTIONS", __module__, resource, nothing, func)
end

macro TRACE(resource, func)
    generate_dispatch("TRACE", __module__, resource, nothing, func)
end

macro PATCH(resource, body, func)
    generate_dispatch("PATCH", __module__, resource, body, func)
end

function generate_dispatch(method, mod, resource_input, body, func)
    resource = typeof(resource_input) <: AbstractString ? resource_input : eval(mod, resource_input)
    # split resource into Val & args
    args = []
    vals_and_args = []
    convert_args_block = quote end
    paths = split(resource, "/"; keep=false)
    for path in paths
        path == "" && continue
        if path[1] == '{' && path[end] == '}'
            # path argument
            arg = path[2:end-1]
            # if type not specified, it's a String by default
            !contains(arg, "::") && (arg *= "::String")
            if arg[1:2] == "::"
                !contains(arg, "::Type{") && throw(ArgumentError("malformed resource argument: $arg"))
                nm = Symbol(replace(string(gensym()), "#", "_"))
                arg = string(nm) * arg
                push!(convert_args_block.args, :($nm = eval(parse($nm))))
                push!(args, nm)
            else
                nm, typ = map(Symbol, split(arg, "::"))
                if typ != :String
                    push!(convert_args_block.args, :($nm = parse($typ, $nm)))
                end
                push!(args, nm)
            end
            string_arg = replace(arg, r"::.+", "::String")
            push!(vals_and_args, parse(string_arg))
        else
            # hard-coded path value
            pathsym = Symbol(path)
            push!(vals_and_args, :(::Type{Val{$(QuoteNode(pathsym))}}))
            Endpoints.PATH_LOOKUPS[path] = Val{pathsym}
        end
    end
    method_val = Type{Val{Symbol(method)}}
    if body != nothing
        spl = split(string(body), "::")
        nm = string(spl[1])
        push!(vals_and_args, parse(nm * "::String"))
        push!(args, Symbol(nm))
        typ = length(spl) > 1 ? Symbol(spl[2]) : :String
        push!(convert_args_block.args, :($(Symbol(nm)) = JSON2.read($(eval(current_module(), typ)), IOBuffer($(Symbol(nm))))))
    end
    return quote
        function $(esc(Symbol(Endpoints))).$(:__uri_dispatch__)(::$(method_val), $(vals_and_args...); query_params...)
            # convert args to expected types
            $(convert_args_block)
            $(esc(func))($(args...); query_params...)
        end
    end
end

function __uri_dispatch__ end

const PATH_LOOKUPS = Dict{String,DataType}()

# get, head, post, put, delete, trace, connect, patch, options
const METHOD_VALS = Dict{HTTP.Method, DataType}(m=>Val{Symbol(m)} for m in instances(HTTP.Method))

function handler404(req, resp)
    resp.status = 404
    resp.body = HTTP.FIFOBuffer(HTTP.resource(HTTP.uri(req)))
    return resp
end

# user can overload to
const BASE_PATH = "/"
function setbasepath!(x::String)
    global BASE_PATH = x
    return
end

const BUF = IOBuffer()

function handler(req, resp)
    vals_and_args = Any[get(PATH_LOOKUPS, seg, seg) for seg in HTTP.splitpath(HTTP.resource(HTTP.uri(req)), length(BASE_PATH) + 1)]
    method_val = METHOD_VALS[HTTP.method(req)]
    query_params = parsequerystring(HTTP.query(HTTP.uri(req)))
    if length(HTTP.body(req)) > 0
        push!(vals_and_args, String(take!(req)))
    end
    local ret
    try
        ret = Endpoints.__uri_dispatch__(method_val, vals_and_args...; query_params...)
    catch e
        if isa(e, MethodError)
            return handler404(req, resp)
        end
        println(e)
        rethrow(e)
    end
    JSON2.write(BUF, ret)
    resp.body = HTTP.FIFOBuffer(take!(BUF))
    return resp
end

#TODO
 # figure out a good way to do user authentication; add filterfunc(request)::Bool capabilities
 # handle 405 method not allowed for paths that match, but method doesn't
 # provide default OPTIONS response
 # find a small Domo service to re-write, see what's missing here
 # Database layer!!
   # see if Persist could be useful first (generic API that other "store" packages could implement)
     # can we get rich enough expression/semantics for parent-child_array operations?
   # do we focus on auto object mapping to tuples/records?
   # @transactional macro
   # function that takes julia args, then executes an SQL query w/ args properly sanitized?
   # be able to run migrations

end # module


text:
# ---
# title: 192. Word Frequency
# id: problem192
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/word-frequency/description/>
# hidden: true
# ---
#
# Write a bash script to calculate the frequency of each word in a text file
# `words.txt`.
#
# For simplicity sake, you may assume:
#
#   * `words.txt` contains only lowercase characters and space `' '` characters.
#   * Each word must consist of lowercase characters only.
#   * Words are separated by one or more whitespace characters.
#
# **Example:**
#
# Assume that `words.txt` has the following content:
#
#
#
#     the day is sunny the the
#     the sunny is is
#
#
# Your script should output the following, sorted by descending frequency:
#
#
#
#     the 4
#     is 3
#     sunny 2
#     day 1
#
#
# **Note:**
#
#   * Don't worry about handling ties, it is guaranteed that each word's frequency count is unique.
#   * Could you write it in one-line using [Unix pipes](http://tldp.org/HOWTO/Bash-Prog-Intro-HOWTO-4.html)?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using RegressionAndOtherStories

"""

Update the ros_df DataFrame contents.

$(SIGNATURES)

## Positional arguments
```julia
* `df=ros_df # DataFrame to be updated according to `./notebooks`
```

Not exported.

"""
function update_ros_df!(df; display_actions=false)

    files = readdir(pwd(); join=true)
    for f in files
        if !(contains(f, ".DS_Store") || contains(f, "intros"))
            if isdir(f) && !contains(f, "intros")
                display_actions && println("Found dir $f")
                add_to_ros_df!(df, f)
            elseif isfile(f)
                display_actions && println("Found notebook $f")
                split_dir_path = split(f, "/")
                parts = split_dir_path[end-1:end]

                append!(df,
                    DataFrame(
                        :chapter => parts[1],
                        :section => parts[2],
                        :reset => false
                    )
                )
            end
        end
    end
end

"""

Add entries to DataFrame `df` based on notebooks found in `dir_path`.

$(SIGNATURES)

## Positional arguments
```julia
* `df::DataFrame` # DataFrame to be updated according to `./notebooks`
* `dir_path::AbstractString` # Directory holding the notebooks
```

## Optional keyword arguments
```julia
* `display_actions=false # Display action steps taken
```

Not exported.

"""
function add_to_ros_df!(df, dir_path;
    display_actions=false)

    files = readdir(dir_path; join=true)
    for f in files
        if !(contains(f, ".DS_Store") || contains(f, "intros"))
            if isdir(f)
                display_actions && println("Found dir $f")
                add_to_ros_df!(df, f)
            elseif isfile(f)
                display_actions && println("Found notebook $f")
                split_dir_path = split(f, "/")
                parts = split_dir_path[end-1:end]

                append!(df,
                    DataFrame(
                        :chapter => parts[1],
                        :section => parts[2],
                        :reset => false
                    )
                )
            end
        end
    end
end

"""

Reset a single notebook.

$(SIGNATURES)

## Positional arguments
```julia
* `fname::AbstractString` # Path to the notebook to be reset
```

## Optional keyword arguments
```julia
* `display_actions=false # Display action steps taken
* `create_pkg_files=false # Not implemented yet!
```

Not exported.

"""
function reset_notebook!(fname;
    display_actions=false,
    create_pkg_files=false)

    f = open(fname)
    c = readlines(f)
    display_actions && length(c) |> display

    d = Vector{String}()
    blocked = false

    for i in 1:length(c)
        if contains(c[i], "╔═╡ 00000000-0000-0000-0000-000000000001")
            blocked = true
            display_actions && println("Blocked at line $i")
        end

        if contains(c[i], "╔═╡ Cell order")
            blocked = false
            display_actions && println("Unblocked at line $i")
        end

        if !blocked
            if i < 3
                display_actions && println(c[i])
            end
            if !contains(c[i], "╟─00000000-0000-0000-0000-000000000001") &&
                !contains(c[i], "╟─00000000-0000-0000-0000-000000000002")

                push!(d, c[i])
            end
        end
    end
    d
end

"""

Update notebook files from DataFrame `df` with `df.reset` set to `true`.

$(SIGNATURES)

## Positional arguments
```julia
* `df=ros_df` # DataFrame with info on all notebooks
```

## Optional keyword arguments
```julia
* `display_actions=false # Display action steps taken
* `create_pkg_files=false # Not implemented yet!
```

Exported.

"""
function update_notebook_files!(df=ros_df;
    display_actions=false,
    create_package_files=false)

    if !isdir("./notebooks")
        @error "You are not in a directory that holds a `notebooks` \
            subdirectory."
        end

    if nrow(df) == 0
        @info "DataFrame ros_df is empty!\n \
            It will be recreated from directory `./notebooks`."

        oldwd = pwd()
        cd("./notebooks")
        update_ros_df!(df)
        println()
        df |> display
        println()
        cd(oldwd)
    end

    if !all(df.reset)
        @info "All ros_df.reset values are false. No actions taken. \
            \nUse `reset_all_notebooks!()`to reset all notebooks \
            \nor set some entries in ros_df.reset to `true` \
            \nand run `update_notebook_files!() again."
    end

    oldwd = pwd()
    isdir(joinpath(oldwd, "notebooks")) && cd("notebooks")
    if !contains(split(pwd(), "/")[end], "notebooks")
        @info "Notebooks directory not found."
        return
    end

    for f in eachrow(df)
        display_actions && println(f)
        if f.reset
            nb = joinpath(".", f.chapter, f.section)
            d = reset_notebook!(nb)

            fname = split(nb, "/")[end]
            new_path = joinpath(".", f.chapter, f.section)
            println("Updating file:$(new_path)")

            open(new_path, "w") do io
                for i in 1:length(d)
                    println(io, d[i])
                end
            end
            f.reset = false

        end
    end
    cd(oldwd)
end

"""

Reset all notebook files from df.

$(SIGNATURES)

## Positional arguments
```julia
* `df=ros_df` # DataFrame with info on all notebooks
```

## Optional keyword arguments
```julia
* `display_actions=false # Display action steps taken
* `create_pkg_files=false # Not implemented yet!
```

Please see [this file](https://github.com/RegressionAndOtherStoriesJulia/RegressionAndOtherStories.jl/blob/main/src/Maintenance/README.md) for how to use these maintenance functions.

Exported.

"""
function reset_all_notebooks!(df=ros_df;
    display_actions=false,
    create_package_files=false)

    if !isdir("./notebooks")
        @error "You are not in a directory that holds a `notebooks` \
            subdirectory."
        return
    end

    if nrow(ros_df) > 0
        ros_df.reset .= true
        update_notebook_files!()
    else
        @info "DataFrame `ros_df` is empty, \nupdate ros_df first by running \
            `update_notebook_files!()`"
        return
    end
end

export
    reset_all_notebooks!,
    update_notebook_files!


text:
# Example 1:
macro macint(ex)
    quote
        println("start")
        $ex
        println("after")
    end
end

@macint println("Where am I?")
#start
#Where am I?
#after

# Example 1B:
macro duplicate(ex)
    quote
        $ex
        $ex
    end
end

@duplicate println(sin(42))
# -0.9165215479156338
# -0.9165215479156338

# Example 2:
macro assert(ex)
    :($ex ? nothing : error("Assertion failed: ", $(string(ex))))
end

@assert 1 == 1.0
# @assert 1 == 42 #> ERROR: Assertion failed: 1 == 42

# Example 3:
macro unless(test, branch)
    quote
        if !$test
            $branch
        end
    end
end

arr = [3.14, 42, 'b']
@unless 41 in arr println("arr does not contain 41")
#> array does not contain 41
@unless 42 in arr println("arr does not contain 42") #> nothing
@unless isempty(arr) println("array arr has elements")
#> array arr has elements

macroexpand(:(@unless 41 in arr println("arr does not contain 41")) )
quote  # none, line 3:
    if !(41 in arr) # line 4:
        println("arr does not contain 41")
    end
end

# Example 4:
arr = ["a", "b", "c"]
macro convarr(arr, T)
    :(reshape($T[$arr...], size($arr)...))
end
@convarr arr Symbol #> returns 3-element Array{Symbol,1} :a :b :c

# Example 5:
@time factorial(10)
#> elapsed time: 3.911e-6 seconds (96 bytes allocated) 3628800

macro timeit(ex)
    quote
        local t0 = time()   # local necessary to not conflict with outside
        local val = $(esc(ex))     # variables; esc makes sure ex is not
        local t1 = time()   # expanded, but is used literally
        print("elapsed time in seconds: ")
        @printf "%.3f" t1 - t0
        val
    end
end
@timeit factorial(10) #> elapsed time in seconds: 0.0003628800
@timeit factorial(15) #> elapsed time in seconds: 0.0001307674368000
a = 42
@timeit a^3 #> elapsed time in seconds: 0.00074088

# Example 6:
macro bench(f)
    quote
        median([@elapsed sleep($f) for x=1:10])
    end
end

@bench 0.1 #> 0.1010291465
@bench 0.2 #> 0.20101891449999998
@bench 0.5 #> 0.501034175

max_stars_repo_path: Module 1/Chapter07/macros.jl

max_issues_repo_path: Module 1/Chapter07/macros.jl

max_forks_repo_path: Module 1/Chapter07/macros.jl


text:
module RembrandtExample
using Scruff
using Scruff.SFuncs
using Scruff.Operators
using Scruff.Algorithms
using Scruff.RTUtils
using Scruff.Utils

# define parameters
isAuthentic_range = [:yes, :no]
isAuthentic_prior = [0.1, 0.9]

size_range = [:small, :medium, :large]
size_CPT = Dict((:yes, :people) => [0.3, 0.2, 0.5], (:no, :people) => [1/3, 1/3, 1/3], (:yes, :landscape) => [0.1, 0.2, 0.7], (:no, :landscape) => [1/3, 1/3, 1/3])

subject_range = [:people, :landscape]
subject_CPT = Dict((:yes, ) => [0.8, 0.2], (:no, ) => [0.5, 0.5])

brightness_range = [:bright, :dark]
brightness_CPT = Dict((:yes, :people) => [0.2, 0.8], (:no, :people) => [0.5, 0.5], (:yes, :landscape) => [0.3, 0.7], (:no, :landscape) => [0.8, 0.2])

function create_network()
    isAuthentic = Cat(isAuthentic_range, isAuthentic_prior)()(:isAuthentic)
    subject = DiscreteCPT(subject_range, subject_CPT)()(:subject)
    size = DiscreteCPT(size_range, size_CPT)()(:size)
    brightness = DiscreteCPT(brightness_range, brightness_CPT)()(:brightness)
    vars = Variable[isAuthentic, subject, size, brightness]
    network = InstantNetwork(vars, VariableGraph(subject=>[isAuthentic],size=>[isAuthentic, subject],brightness => [isAuthentic, subject]))
    return network
end

function run_inference(network, alg::InstantAlgorithm)
    runtime = Runtime(network) # create runtime
    # default_initializer(runtime) # create instances

    # Set soft evidence.
    evid = SoftScore([:bright, :dark], [0.2, 0.8])
    # Set hard evidence.
    #evid = HardScore(:dark)
    # post_evidence!(runtime, inst_brightness, evid)
    # println("Evidence $(evid) applied to $(:brightness)")

    # Perform inference
    infer(alg, runtime, Dict{Symbol, Score}(:brightness => evid))

    # Get current instances.
    inst_isAuthentic = current_instance(runtime, get_node(network, :isAuthentic))
    inst_subject = current_instance(runtime, get_node(network, :subject))
    inst_size = current_instance(runtime, get_node(network, :size))
    inst_brightness = current_instance(runtime, get_node(network, :brightness))

    # Get updated beliefs
    println("Probability that it is authenthic: $(probability(alg, runtime, inst_isAuthentic, :yes))")
    println("Probability that it is small: $(probability(alg, runtime, inst_size, :small))")
    println("Probability that it is medium: $(probability(alg, runtime, inst_size, :medium))")
    println("Probability that it is large: $(probability(alg, runtime, inst_size, :large))")
    println("Probability that it is dark: $(probability(alg, runtime, inst_brightness, :dark))")
    println("Probability that it is bright: $(probability(alg, runtime, inst_brightness, :bright))")
    println("Probability that it has people: $(probability(alg, runtime, inst_subject, :people))")
    println("Probability that it has landscape: $(probability(alg, runtime, inst_subject, :landscape))")
end

# create network
network = create_network()

# run inference
println("##################################")
println("Running Three Pass Belief Propagation")
println("##################################")
run_inference(network, ThreePassBP())
println("Running Loopy Belief Propagation")
println("##################################")
run_inference(network, LoopyBP())
println("##################################")
println("Running Variable Elimination")
println("##################################")
run_inference(network, VE(get_variables(network)))
println("##################################")
println("Running Importance Sampling")
println("##################################")
prop = convert(Dict{Symbol, SFunc}, Dict(:subject => DiscreteCPT(subject_range, subject_CPT)))
run_inference(network,  Importance(make_custom_proposal(prop), 1000))
end


text:
# This file contains the Pins to connect the links

"""
    AbstractPin{T}

Abstract type of `Outpin` and `Inpin`. See also: [`Outpin`](@ref), [`Inpin`](@ref)
"""
abstract type AbstractPin{T} end

"""
    Outpin{T}()

Constructs and `OutPut` pin. The data flow from `Outpin` is outwards from the pin i.e., data is written from `OutPort` to its links.
"""
mutable struct Outpin{T} <: AbstractPin{T}
    links::Union{Vector{Link{T}}, Missing}
    id::UUID
    # NOTE: When Outpin is initialized, its links are missing.
    # The existance of links of Outpin is used to determine
    # whether the Outpin is bound or not.
    Outpin{T}() where {T} = new{T}(missing, uuid4())
end
Outpin() = Outpin{Float64}()

show(io::IO, outpin::Outpin) = print(io, "Outpin(eltype:$(eltype(outpin)), isbound:$(isbound(outpin)))")

"""
    Inpin{T}()

Constructs and `InPut` pin. The data flow from `Inpin` is inwards to the pin i.e., data is read from links of `InPort`.
"""
mutable struct Inpin{T} <: AbstractPin{T}
    link::Union{Link{T}, Missing}
    id::UUID
    # NOTE: When an Inpin is initialized, its link is missing.
    # The state of link of the Inpin is used to decide whether the Inpin is bound or not.
    Inpin{T}() where {T} = new{T}(missing, uuid4())
end
Inpin() = Inpin{Float64}()

show(io::IO, inpin::Inpin) = print(io, "Inpin(eltype:$(eltype(inpin)), isbound:$(isbound(inpin)))")

"""
    bind(link::Link, pin)

Binds `link` to `pin`. When bound, data written into or read from `pin` is written into or read from `link`.
"""
bind(link::Link, inpin::Inpin) = (inpin.link = link; link.slaveid = inpin.id)
bind(link::Link, outpin::Outpin) = (outpin.links === missing ? (outpin.links = [link]) : push!(outpin.links, link); link.masterid = outpin.id)

"""
    isbound(pin::AbstractPin)

Returns `true` if `pin` is bound to other pins.
"""
function isbound(outpin::Outpin)
    outpin.links === missing && return false
    !isempty(outpin.links)
end
isbound(inpin::Inpin) = inpin.link !== missing

"""
    eltype(pin::AbstractPin)

Returns element typef of pin.
"""
eltype(pin::AbstractPin{T}) where T = T

"""
    take!(pin::Inpin)

Takes data from `pin`. The data is taken from the links of `pin`.

!!! warning
    To take data from `pin`, a running task that puts data must be bound to `link` of `pin`.

# Example
```jldoctest
julia> ip = Inpin();

julia> l = Link();

julia> bind(l, ip);

julia> t = @async for item in 1 : 5
       put!(l, item)
       end;

julia> take!(ip)
1.0

julia> take!(ip)
2.0
```
"""
take!(pin::Inpin) = take!(pin.link)

"""
    put!(pin::Outpin, val)

Puts `val` to `pin`. `val` is put into the links of `pin`.

!!! warning
    To take data from `pin`, a running task that puts data must be bound to `link` of `pin`.

# Example
```jldoctest
julia> op = Outpin();

julia> l = Link();

julia> bind(l, op);

julia> t = @async while true
       val = take!(l)
       val === NaN && break
       println("Took " * string(val))
       end;

julia> put!(op, 1.)
Took 1.0

julia> put!(op, 3.)
Took 3.0

julia> put!(op, NaN)
```
"""
put!(pin::Outpin, val) = foreach(link -> put!(link, val), pin.links)

##### Connecting and disconnecting links
# #
# # This `iterate` function is dummy. It is defined just for `[l...]` to be written.
# #
# iterate(l::AbstractPin, i=1) = i > 1 ? nothing : (l, i + 1)

"""
    connect!(outpin::Link, inpin::Link)

Connects `outpin` to `inpin`. When connected, any element that is put into `outpin` is also put into `inpin`.

    connect!(outpin::AbstractVector{<:Link}, inpin::AbstractVector{<:Link})

Connects each link in `outpin` to each link in `inpin` one by one. See also: [`disconnect!`](@ref)

# Example
```jldoctest
julia> op, ip = Outpin(), Inpin();

julia> l = connect!(op, ip)
Link(state:open, eltype:Float64, isreadable:false, iswritable:false)

julia> l in op.links
true

julia> ip.link === l
true
```
"""
function connect!(outpin::Outpin, inpin::Inpin)
    # NOTE: The connecion of an `Outpin` to multiple `Inpin`s is possible since an `Outpin` may drive multiple
    # `Inpin`s. However, the connection of multiple `Outpin`s to the same `Inpin` is NOT possible since an `Inpin`
    # can be driven by a single `Outpin`.
    isbound(inpin) && error("$inpin is already bound. No new connections.")
    isconnected(outpin, inpin) && (@warn "$outpin and $inpin are already connected."; return)

    link = Link{promote_type(eltype(outpin), eltype(inpin))}()
    bind(link, outpin)
    bind(link, inpin)
    return link
end
connect!(outpins::AbstractVector{<:Outpin}, inpins::AbstractVector{<:Inpin}) = connect!.(outpins, inpins)

"""
    disconnect!(link1::Link, link2::Link)

Disconnects `link1` and `link2`. The order of arguments is not important. See also: [`connect!`](@ref)
"""
function disconnect!(outpin::Outpin, inpin::Inpin)
    outpin.links === missing || deleteat!(outpin.links, findall(link -> link == inpin.link, outpin.links))
    inpin.link = missing
    # inpin.link = Link{eltype(inpin)}()
end
disconnect!(outpins::AbstractVector{<:Outpin}, inpins::AbstractVector{<:Inpin}) = (disconnect!.(outpins, inpins); nothing)

"""
    isconnected(link1, link2)

Returns `true` if `link1` is connected to `link2`. The order of the arguments are not important.
See also [`connect!`](@ref), [`disconnect!`](@ref)
"""
function isconnected(outpin::Outpin, inpin::Inpin)
    if !isbound(outpin) || !isbound(inpin)
        return false
    else
        inpin.link in [link for link in outpin.links]
    end
end
isconnected(outpins::AbstractVector{<:Outpin}, inpins::AbstractVector{<:Inpin}) = all(isconnected.(outpins, inpins))

# ------------------------ Deprecations  -----------------------

#= NOTE:
The methods

    connect!(outpins, inpins) = connect!([outpins...], [inpins...])
    disconnect!(outpins, inpins) = disconnect!([outpins...], [inpins...])
    isconnected(outpins, inpins) = isconnected([outpins...], [inpins...])

are ambiguis. Since these methods throws StackOverflowError when called with `outpins` are `Inpin`s and
inpins are `Outpin`s. So, there is not need for

    iterate(l::AbstractPin, i=1) = i > 1 ? nothing : (l, i + 1)

method-
=#

# """
#     UnconnectedLinkError <: Exception

# Exception thrown when the links are not connected to each other.
# """
# struct UnconnectedLinkError <: Exception
#     msg::String
# end
# Base.showerror(io::IO, err::UnconnectedLinkError) = print(io, "UnconnectedLinkError:\n $(err.msg)")


text:
# Type for pushing and popping text states
struct CrayonStack
    incremental::Bool
    crayons::Vector{Crayon}
end

Base.print(io::IO, cs::CrayonStack) = print(io, cs.crayons[end])

function CrayonStack(; incremental::Bool = false)
    CrayonStack(incremental, [Crayon(ANSIColor(0x9, COLORS_16, !incremental),
                                     ANSIColor(0x9, COLORS_16, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental),
                                     ANSIStyle(false, !incremental))])
end

# Checks if equal disregarding active or not
_equal(a::ANSIColor, b::ANSIColor) = a.r == b.r && a.g == b.g && a.b == b.b && a.style == b.style
_equal(a::ANSIStyle, b::ANSIStyle) = a.on == b.on

# Currently we have the crayon a on the stack.
# We now want to push the crayon b.
# If in incremental mode we compute the changes needed to achive the state of b.
function _incremental_add(a::ANSIColor, b::ANSIColor, incremental::Bool)
    if b.active
        return ANSIColor(b.r, b.g, b.b, b.style, !_equal(a, b) || !incremental)
    else
        return ANSIColor(a.r, a.g, a.b, a.style, !incremental)
    end
end

# Similar to above
_incremental_add(a::ANSIStyle, b::ANSIStyle, incremental::Bool) = b.active ? ANSIStyle(b.on, !_equal(a, b) || !incremental) : ANSIStyle(a.on, !incremental)

# Have a, going to b
# State should be exactly as b!
# However, we only want to activate those that are needed, that is
# those that are active in a AND are different from b
_incremental_sub(a::ANSIColor, b::ANSIColor, incremental::Bool) = ANSIColor(b.r, b.g, b.b, b.style, !_equal(a, b) || !incremental)
_incremental_sub(a::ANSIStyle, b::ANSIStyle, incremental::Bool) = ANSIStyle(b.on, !_equal(a, b) || !incremental)

function Base.push!(cs::CrayonStack, c::Crayon)
    pc = cs.crayons[end]
    push!(cs.crayons, Crayon(
        _incremental_add(pc.fg           , c.fg           , cs.incremental),
        _incremental_add(pc.bg           , c.bg           , cs.incremental),
        _incremental_add(pc.reset        , c.reset        , cs.incremental),
        _incremental_add(pc.bold         , c.bold         , cs.incremental),
        _incremental_add(pc.faint        , c.faint        , cs.incremental),
        _incremental_add(pc.italics      , c.italics      , cs.incremental),
        _incremental_add(pc.underline    , c.underline    , cs.incremental),
        _incremental_add(pc.blink        , c.blink        , cs.incremental),
        _incremental_add(pc.negative     , c.negative     , cs.incremental),
        _incremental_add(pc.conceal      , c.conceal      , cs.incremental),
        _incremental_add(pc.strikethrough, c.strikethrough, cs.incremental)))
    return cs
end

function Base.pop!(cs::CrayonStack)
    length(cs.crayons) == 1 && throw(ArgumentError("no more Crayons left in stack"))

    c = pop!(cs.crayons)
    pc = cs.crayons[end]
    if length(cs.crayons) == 1
        pc = Crayon(ANSIColor(0x9, COLORS_16, true),
                    ANSIColor(0x9, COLORS_16, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true),
                    ANSIStyle(false, true))
    end
    cs.crayons[end] = Crayon(
        _incremental_sub(c.fg           , pc.fg           , cs.incremental),
        _incremental_sub(c.bg           , pc.bg           , cs.incremental),
        _incremental_sub(c.reset        , pc.reset        , cs.incremental),
        _incremental_sub(c.bold         , pc.bold         , cs.incremental),
        _incremental_sub(c.faint        , pc.faint        , cs.incremental),
        _incremental_sub(c.italics      , pc.italics      , cs.incremental),
        _incremental_sub(c.underline    , pc.underline    , cs.incremental),
        _incremental_sub(c.blink        , pc.blink        , cs.incremental),
        _incremental_sub(c.negative     , pc.negative     , cs.incremental),
        _incremental_sub(c.conceal      , pc.conceal      , cs.incremental),
        _incremental_sub(c.strikethrough, pc.strikethrough, cs.incremental))
    # Return the currently active crayon so we can use print(pop!(crayonstack), "bla")
    return cs
end


text:
using RangeArrays
using Base.Test

# Regular RangeMatrix tests

R = RangeMatrix(UnitRange{Int}[1:10, 11:20, 21:30, 31:40])

@test vec(collect(R)) == collect(1:40)
@test R[1,1] == 1

for i=1:40
    @test R[i] == i
end
j = 0
for i in eachindex(R)
    j+=1
    @test R[i] == j
end

@test_throws BoundsError R[0]
@test_throws BoundsError R[41]

@test_throws BoundsError R[:, 0]
@test R[:, 1] == 1:10
@test R[:, 2] == 11:20
@test R[:, 3] == 21:30
@test R[:, 4] == 31:40
@test_throws BoundsError R[:, 5]

@test R[:, 1:2] == RangeMatrix(1:10, 11:20)
@test R[:, 1:4] == R[:,:] == RangeMatrix(1:10, 11:20, 21:30, 31:40)
@test R[:, 3:4] == RangeMatrix(21:30, 31:40)

# RepeatedRangeMatrix tests
R = RepeatedRangeMatrix(1:10, 0:10:30)

@test vec(collect(R)) == collect(1:40)
@test R[1,1] == 1

for i=1:40
    @test R[i] == i
end
j = 0
for i in eachindex(R)
    j+=1
    @test R[i] == j
end

@test_throws BoundsError R[0]
@test_throws BoundsError R[41]

@test_throws BoundsError R[:, 0]
@test R[:, 1] == 1:10
@test R[:, 2] == 11:20
@test R[:, 3] == 21:30
@test R[:, 4] == 31:40
@test_throws BoundsError R[:, 5]

@test R[:, 1:2] == RangeMatrix(1:10, 11:20)
@test R[:, 1:4] == R[:,:] == RangeMatrix(1:10, 11:20, 21:30, 31:40)
@test R[:, 3:4] == RangeMatrix(21:30, 31:40)