text:
using SnoopCompile
using SnoopCompile.SnoopCompileCore
using Test
using InteractiveUtils
using Random
using Profile
using MethodAnalysis
using Core: MethodInstance
# using PyPlot: PyPlot, plt    # uncomment to test visualizations

using SnoopCompile.FlameGraphs.AbstractTrees  # For FlameGraphs tests

@testset "@snoopi_deep" begin
    # WARMUP (to compile all the small, reachable methods)
    @eval module M  # Example with some functions that include type instability
        i(x) = x+5
        h(a::Array) = i(a[1]::Integer) + 2
        g(y::Integer) = h(Any[y])
    end
    M.g(2)  # Warmup all deeply reachable functions
    M.g(true)

    # Redefine the module, so the snoop will only show these functions:
    @eval module M  # Example with some functions that include type instability
        i(x) = x+5
        h(a::Array) = i(a[1]::Integer) + 2
        g(y::Integer) = h(Any[y])
    end

    tinf = @snoopi_deep begin
        M.g(2)
        M.g(true)
    end
    @test SnoopCompile.isROOT(Core.MethodInstance(tinf))
    @test SnoopCompile.isROOT(Method(tinf))
    @test isempty(staleinstances(tinf))
    frames = flatten(tinf)
    @test length(frames) == 7  # ROOT, g(::Int), g(::Bool), h(...), i(::Integer), i(::Int), i(::Bool)
    @test issorted(frames; by=exclusive)
    names = [Method(frame).name for frame in frames]
    @test sort(names) == [:ROOT, :g, :g, :h, :i, :i, :i]
    mg = which(M.g, (Int,))
    tinfsg = collect_for(mg, tinf)
    @test length(tinfsg) == 2
    @test all(node -> Method(node) == mg, tinfsg)
    mig = MethodInstance(first(tinfsg))
    tinfsg1 = collect_for(mig, tinf)
    @test length(tinfsg1) == 1
    @test MethodInstance(tinfsg1[1]) == mig
    @test all(node -> Method(node) == mg, tinfsg)

    longest_frame_time = exclusive(frames[end])
    @test length(flatten(tinf, tmin=longest_frame_time)) == 1

    frames_unsorted = flatten(tinf; sortby=nothing)
    ifi = frames_unsorted[1].mi_info
    @test SnoopCompile.isROOT(Core.MethodInstance(ifi))
    @test SnoopCompile.isROOT(Method(ifi))
    names = [Method(frame).name for frame in frames_unsorted]
    argtypes = [MethodInstance(frame).specTypes.parameters[2] for frame in frames_unsorted[2:end]]
    @test names == [:ROOT, :g, :h,           :i,      :i,  :g,   :i]
    @test argtypes == [    Int, Vector{Any}, Integer, Int, Bool, Bool]

    timesm = accumulate_by_source(frames)
    @test length(timesm) == 4
    names = [m.name for (time, m) in timesm]
    @test sort(names) == [:ROOT, :g, :h, :i]
    longest_method_time = timesm[end][1]
    @test length(accumulate_by_source(frames; tmin=longest_method_time)) == 1

    @test SnoopCompile.isROOT(Core.MethodInstance(tinf))
    @test SnoopCompile.isROOT(Method(tinf))
    iframes = flatten(tinf; sortby=inclusive)
    @test issorted(iframes; by=inclusive)

    t = map(inclusive, frames_unsorted)
    @test t[2] >= t[3] >= t[4]
    ifi = frames_unsorted[2].mi_info
    @test Core.MethodInstance(ifi).def == Method(ifi) == which(M.g, (Int,))
    names = [Method(frame).name for frame in frames_unsorted]
    argtypes = [MethodInstance(frame).specTypes.parameters[2] for frame in frames_unsorted[2:end]]
    @test names == [:ROOT, :g, :h,           :i,      :i,  :g,   :i]
    @test argtypes == [    Int, Vector{Any}, Integer, Int, Bool, Bool]

    # Also check module-level thunks
    @eval module M  # Example with some functions that include type instability
        i(x) = x+5
        h(a::Array) = i(a[1]::Integer) + 2
        g(y::Integer) = h(Any[y])
    end
    tinfmod = @snoopi_deep begin
        @eval @testset "Outer" begin
            @testset "Inner" begin
                for i = 1:2 M.g(2) end
            end
        end
    end
    frames = flatten(tinfmod)
    timesm = accumulate_by_source(frames)
    timesmod = filter(pr -> isa(pr[2], Core.MethodInstance), timesm)
    @test length(timesmod) == 1
end

# For the higher-order function attribution test, we need to prevent `f2`
# from being passed via closure, so we define these globally.
fdouble(x) = 2x

@testset "inference_triggers" begin
    myplus(x, y) = x + y    # freshly redefined even if tests are re-run
    function f(x)
        x < 0.25 ? 1 :
        x < 0.5  ? 1.0 :
        x < 0.75 ? 0x01 : Float16(1)
    end
    g(c) = myplus(f(c[1]), f(c[2]))
    tinf = @snoopi_deep g([0.7, 0.8])
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf; exclude_toplevel=false)
    @test length(itrigs) == 2
    @test suggest(itrigs[1]).categories == [SnoopCompile.FromTestDirect]
    s = suggest(itrigs[2])
    @test SnoopCompile.FromTestCallee ∈ s.categories
    @test SnoopCompile.UnspecCall ∈ s.categories
    @test occursin("myplus", string(MethodInstance(itrigs[2].node).def.name))
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    @test filtermod(@__MODULE__, itrigs) == [itrig]
    @test isempty(filtermod(Base, itrigs))
    io = IOBuffer()
    show(io, itrig)
    str = String(take!(io))
    @test occursin(r"MethodInstance for .*myplus.*\(::UInt8, ::Float16\)", str)
    @test occursin("from g", str)
    @test occursin(r"with specialization .*::Vector\{Float64\}", str)
    mis = callerinstance.(itrigs)
    @test only(mis).def == which(g, (Any,))
    @test callingframe(itrig).callerframes[1].func === :eval
    @test_throws ArgumentError("it seems you've supplied a child node, but backtraces are collected only at the entrance to inference") inference_triggers(tinf.children[1])
    @test stacktrace(itrig) isa Vector{StackTraces.StackFrame}
    itrig0 = itrig
    counter = 0
    while !isempty(itrig.callerframes) && counter < 1000  # defensively prevent infinite loop
        itrig = callingframe(itrig)
        counter += 1
    end
    @test counter < 1000
    show(io, itrig)
    str = String(take!(io))
    @test occursin("called from toplevel", str)
    @test itrig != itrig0
    @test InferenceTrigger(itrig) == itrig0

    # Tree generation
    itree = trigger_tree(itrigs)
    @test length(itree.children) == 1
    @test isempty(itree.children[1].children)
    print_tree(io, itree)
    @test occursin(r"myplus.*UInt8.*Float16", String(take!(io)))

    # Where the caller is inlined into something else
    callee(x) = 2x
    @inline caller(c) = callee(c[1])
    callercaller(cc) = caller(cc[1])
    callercaller([Any[1]])
    cc = [Any[0x01]]
    tinf = @snoopi_deep callercaller(cc)
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    show(io, itrig)
    str = String(take!(io))
    @test occursin(r"to call MethodInstance for .*callee.*\(::UInt8\)", str)
    @test occursin("from caller", str)
    @test occursin(r"inlined into MethodInstance for .*callercaller.*\(::Vector{Vector{Any}}\)", str)
    s = suggest(itrig)
    @test !isignorable(s)
    print(io, s)
    @test occursin(r"snoopi_deep\.jl:\d+: non-inferrable or unspecialized call.*::UInt8", String(take!(io)))

    mysqrt(x) = sqrt(x)
    c = Any[1, 1.0, 0x01, Float16(1)]
    tinf = @snoopi_deep map(mysqrt, c)
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf)
    itree = trigger_tree(itrigs)
    io = IOBuffer()
    print_tree(io, itree)
    @test occursin(r"mysqrt.*Float64", String(take!(io)))
    print(io, itree)
    @test String(take!(io)) == "TriggerNode for root with 2 direct children"
    @test length(flatten(itree)) > length(c)
    length(suggest(itree).children) == 2
    loctrigs = accumulate_by_source(itrigs)
    show(io, loctrigs)
    @test any(str->occursin("4 callees from 2 callers", str), split(String(take!(io)), '\n'))
    @test filtermod(Base, loctrigs) == loctrigs
    @test isempty(filtermod(@__MODULE__, loctrigs))
    # This actually tests the suggest framework a bit, but...
    for loctrig in loctrigs
        summary(io, loctrig)
    end
    str = String(take!(io))
    @test occursin("1 callees from 1 callers, consider improving inference", str)
    @test occursin("4 callees from 2 callers, consider despecializing the callee", str)
    @test occursin("non-inferrable or unspecialized call", str)
    @test occursin("partial type call", str)
    mtrigs = accumulate_by_source(Method, itrigs)
    for mtrig in mtrigs
        show(io, mtrig)
    end
    str = String(take!(io))
    @test occursin("map(f, A::AbstractArray) in Base", str)
    @test occursin("2 callees from 1 caller", str)
    for mtrig in mtrigs
        summary(io, mtrig)
    end
    str = String(take!(io))
    @test occursin(r"map.*had 1 specialization", str)
    @test occursin(r"calling Base\.Generator", str)
    @test occursin("calling mysqrt (3 instances)", str)
    modtrigs = SnoopCompile.parcel(mtrigs)
    @test only(modtrigs).first === Base
    @test filtermod(Base, mtrigs) == mtrigs
    @test isempty(filtermod(M, mtrigs))

    # Multiple callees on the same line
    fline(x) = 2*x[]
    gline(x) = x[]
    fg(x) = fline(gline(x[]))
    cc = Ref{Any}(Ref{Base.RefValue}(Ref(3)))
    tinf = @snoopi_deep fg(cc)
    itrigs = inference_triggers(tinf)
    loctrigs = accumulate_by_source(itrigs)
    @test length(loctrigs) == 2
    loctrigs[1].tag == loctrigs[2].tag

    # Higher order function attribution
    @noinline function mymap!(f, dst, src)
        for i in eachindex(dst, src)
            dst[i] = f(src[i])
        end
        return dst
    end
    @noinline mymap(f::F, src) where F = mymap!(f, Vector{Any}(undef, length(src)), src)
    callmymap(x) = mymap(fdouble, x)
    callmymap(Any[1, 2])  # compile all for one set of types
    x = Any[1.0, 2.0]   # fdouble not yet inferred for Float64
    tinf = @snoopi_deep callmymap(x)
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    @test occursin(r"with specialization MethodInstance for .*mymap!.*\(::.*fdouble.*, ::Vector{Any}, ::Vector{Any}\)", string(itrig))
    @test occursin(r"with specialization MethodInstance for .*mymap.*\(::.*fdouble.*, ::Vector{Any}\)", string(callingframe(itrig)))
    @test occursin(r"with specialization MethodInstance for .*callmymap.*\(::Vector{Any}\)", string(skiphigherorder(itrig)))
    # Ensure we don't skip non-higher order calls
    callfdouble(c) = fdouble(c[1])
    callfdouble(Any[1])
    c = Any[Float16(1)]
    tinf = @snoopi_deep callfdouble(c)
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    @test skiphigherorder(itrig) == itrig

    # With a closure
    @eval module M
        function f(c, name)
            stringx(x) = string(x) * name
            stringx(x::Int) = string(x) * name
            stringx(x::Float64) = string(x) * name
            stringx(x::Bool) = string(x) * name

            n = 0
            for x in c
                n += length(stringx(x))
            end
            return n
        end
    end
    c = Any["hey", 7]
    tinf = @snoopi_deep M.f(c, " there")
    itrigs = inference_triggers(tinf)
    @test length(itrigs) > 1
    mtrigs = accumulate_by_source(Method, itrigs)
    summary(io, only(mtrigs))
    str = String(take!(io))
    @test occursin(r"closure.*stringx.*\{String\} at", str)
end

@testset "suggest" begin
    categories(tinf) = suggest(only(inference_triggers(tinf))).categories

    io = IOBuffer()

    # UnspecCall and relation to Test
    @eval module M
        callee(x) = 2x
        caller(c) = callee(c[1])
    end
    tinf = @snoopi_deep M.caller(Any[1])
    itrigs = inference_triggers(tinf; exclude_toplevel=false)
    @test length(itrigs) == 2
    s = suggest(itrigs[1])
    @test s.categories == [SnoopCompile.FromTestDirect]
    show(io, s)
    @test occursin(r"called by Test.*ignore", String(take!(io)))
    s = suggest(itrigs[2])
    @test s.categories == [SnoopCompile.FromTestCallee, SnoopCompile.UnspecCall]
    show(io, s)
    @test occursin(r"non-inferrable or unspecialized call.*annotate caller\(c::Vector\{Any\}\) at snoopi_deep.*callee\(::Int", String(take!(io)))

    # Same test, but check the test harness & inlineable detection
    @eval module M
        callee(x) = 2x
        @inline caller(c) = callee(c[1])
    end
    cats = categories(@snoopi_deep M.caller(Any[1]))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.CallerInlineable, SnoopCompile.UnspecCall]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin("non-inferrable or unspecialized call", String(take!(io)))

    # UnspecType
    @eval module M
        struct Container{L,T} x::T end
        Container(x::T) where {T} = Container{length(x),T}(x)
    end
    cats = categories(@snoopi_deep M.Container([1,2,3]))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.UnspecType]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin("partial type call", String(take!(io)))
    @eval module M
        struct Typ end
        struct Container{N,T} x::T end
        Container{N}(x::T) where {N,T} = Container{N,T}(x)
        typeconstruct(c) = Container{3}(c[])
    end
    c = Ref{Any}(M.Typ())
    cats = categories(@snoopi_deep M.typeconstruct(c))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.UnspecType]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    # println(String(take!(io)))
    @test occursin("partial type call", String(take!(io)))

    # Invoke
    @eval module M
        @noinline callf(@nospecialize(f::Function), x) = f(x)
        g(x) = callf(sqrt, x)
    end
    cats = categories(@snoopi_deep M.g(3))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.CallerInlineable, SnoopCompile.Invoke]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"invoked callee.*may fail to precompile", String(take!(io)))

    # CalleeVariable
    mysin(x) = 1
    mycos(x) = 2
    docall(ref, x) = ref[](x)
    function callvar(x)
        fref = Ref{Any}(rand() < 0.5 ? mysin : mycos)
        return docall(fref, x)
    end
    cats = categories(@snoopi_deep callvar(0.2))
    @test cats == [SnoopCompile.CalleeVariable]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"variable callee.*avoid assigning function to variable", String(take!(io)))
    # CalleeVariable as an Argument
    @eval module M
        mysin(x) = 1
        mycos(x) = 2
        mytan(x) = 3
        mycsc(x) = 4
        getfunc(::Int) = mysin
        getfunc(::Float64) = mycos
        getfunc(::Char) = mytan
        getfunc(::String) = mycsc
        docall(@nospecialize(f), x) = f(x)
        function callvar(ref, f=nothing)
            x = ref[]
            if f === nothing
                f = getfunc(x)
            end
            return docall(f, x)
        end
    end
    tinf = @snoopi_deep M.callvar(Ref{Any}(0.2))
    cats = suggest(inference_triggers(tinf)[end]).categories
    @test cats == [SnoopCompile.CalleeVariable]
    # CalleeVariable & varargs
    @eval module M
        f1va(a...) = 1
        f2va(a...) = 2
        docallva(ref, x) = ref[](x...)
        function callsomething(args)
            fref = Ref{Any}(rand() < 0.5 ? f1va : f2va)
            docallva(fref, args)
        end
    end
    cats = categories(@snoopi_deep M.callsomething(Any['a', 2]))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.CallerInlineable, SnoopCompile.CalleeVariable]
    @eval module M
        f1va(a...) = 1
        f2va(a...) = 2
        @noinline docallva(ref, x) = ref[](x...)
        function callsomething(args)
            fref = Ref{Any}(rand() < 0.5 ? f1va : f2va)
            docallva(fref, args)
        end
    end
    cats = categories(@snoopi_deep M.callsomething(Any['a', 2]))
    @test cats == [SnoopCompile.CalleeVariable]

    # CallerVararg
    @eval module M
        f1(x) = 2x
        c1(x...) = f1(x[2])
    end
    c = Any['c', 1]
    cats = categories(@snoopi_deep M.c1(c...))
    @test SnoopCompile.CalleeVararg ∉ cats
    @test SnoopCompile.CallerVararg ∈ cats
    @test SnoopCompile.UnspecCall ∈ cats
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"non-inferrable or unspecialized call with vararg caller.*annotate", String(take!(io)))

    # CalleeVararg
    @eval module M
        f2(x...) = 2*x[2]
        c2(x) = f2(x...)
    end
    cats = categories(@snoopi_deep M.c2(c))
    @test SnoopCompile.CallerVararg ∉ cats
    @test SnoopCompile.CalleeVararg ∈ cats
    @test SnoopCompile.UnspecCall ∈ cats
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"non-inferrable or unspecialized call with vararg callee.*annotate", String(take!(io)))

    # InvokeCalleeVararg
    @eval module M
        struct AType end
        struct BType end
        Base.show(io::IO, ::AType) = print(io, "A")
        Base.show(io::IO, ::BType) = print(io, "B")
        @noinline doprint(ref) = print(IOBuffer(), "a", ref[], 3.2)
    end
    cats = categories(@snoopi_deep M.doprint(Ref{Union{M.AType,M.BType}}(M.AType())))
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.InvokedCalleeVararg]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"invoked callee is varargs.*homogenize", String(take!(io)))

    # Vararg that resolves to a UnionAll
    @eval module M
        struct ArrayWrapper{T,N,A,Args} <: AbstractArray{T,N}
            data::A
            args::Args
        end
        ArrayWrapper{T}(data, args...) where T = ArrayWrapper{T,ndims(data),typeof(data),typeof(args)}(data, args)
        @noinline makewrapper(data::AbstractArray{T}, args) where T = ArrayWrapper{T}(data, args...)
    end
    # run and redefine for reproducible results
    M.makewrapper(rand(2,2), ["a", 'b', 5])
    @eval module M
        struct ArrayWrapper{T,N,A,Args} <: AbstractArray{T,N}
            data::A
            args::Args
        end
        ArrayWrapper{T}(data, args...) where T = ArrayWrapper{T,ndims(data),typeof(data),typeof(args)}(data, args)
        @noinline makewrapper(data::AbstractArray{T}, args) where T = ArrayWrapper{T}(data, args...)
    end
    tinf = @snoopi_deep M.makewrapper(rand(2,2), ["a", 'b', 5])
    itrigs = inference_triggers(tinf)
    @test length(itrigs) == 2
    s = suggest(itrigs[1])
    @test s.categories == [SnoopCompile.FromTestCallee, SnoopCompile.UnspecType]
    print(io, s)
    @test occursin("partial type call", String(take!(io)))
    s = suggest(itrigs[2])
    @test s.categories == [SnoopCompile.CallerVararg, SnoopCompile.UnspecType]
    print(io, s)
    @test occursin(r"partial type call with vararg caller.*ignore.*annotate", String(take!(io)))
    mtrigs = accumulate_by_source(Method, itrigs)
    for mtrig in mtrigs
        summary(io, mtrig)
    end
    str = String(take!(io))
    @test occursin("makewrapper(data", str)
    @test occursin("ArrayWrapper{Float64", str)
    @test occursin("Tuple{String", str)

    # ErrorPath
    @eval module M
        struct MyType end
        struct MyException <: Exception
            info::Vector{MyType}
        end
        MyException(obj::MyType) = MyException([obj])
        @noinline function checkstatus(b::Bool, info)
            if !b
                throw(MyException(info))
            end
            return nothing
        end
    end
    cats = categories(@snoopi_deep try M.checkstatus(false, M.MyType()) catch end)
    @test cats == [SnoopCompile.FromTestCallee, SnoopCompile.ErrorPath]
    SnoopCompile.show_suggest(io, cats, nothing, nothing)
    @test occursin(r"error path.*ignore", String(take!(io)))

    # Core.Box
    @eval module M
        struct MyInt <: Integer end
        Base.:(*)(::MyInt, r::Int) = 7*r
        function abmult(r::Int, z)  # from https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
            if r < 0
                r = -r
            end
            f = x -> x * r
            return f(z)
        end
    end
    z = M.MyInt()
    tinf = @snoopi_deep M.abmult(3, z)
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    s = suggest(itrig)
    @test SnoopCompile.HasCoreBox ∈ s.categories
    print(io, s)
    @test occursin(r"Core\.Box.*fix this.*http", String(take!(io)))
    mtrigs = accumulate_by_source(Method, itrigs)
    summary(io, only(mtrigs))
    @test occursin(r"Core\.Box.*fix this.*http", String(take!(io)))

    # Test one called from toplevel
    fromtask() = (while false end; 1)
    tinf = @snoopi_deep wait(@async fromtask())
    @test isempty(staleinstances(tinf))
    itrigs = inference_triggers(tinf)
    itrig = only(itrigs)
    s = suggest(itrig)
    @test s.categories == [SnoopCompile.NoCaller]
    itree = trigger_tree(itrigs)
    print_tree(io, itree)
    @test occursin(r"{var\"#fromtask", String(take!(io)))
    print(io, s)
    occursin(r"unknown caller.*Task", String(take!(io)))
    mtrigs = accumulate_by_source(Method, itrigs)
    @test isempty(mtrigs)

    # Empty
    SnoopCompile.show_suggest(io, SnoopCompile.Suggestion[], nothing, nothing)
    @test occursin("Unspecialized or unknown for", String(take!(io)))

    # Printing says *something* for any set of categories
    annots = instances(SnoopCompile.Suggestion)
    iter = [1:2 for _ in 1:length(annots)]
    cats = SnoopCompile.Suggestion[]
    for state in Iterators.product(iter...)
        empty!(cats)
        for (i, s) in enumerate(state)
            if s == 2
                push!(cats, annots[i])
            end
        end
        SnoopCompile.show_suggest(io, cats, nothing, nothing)
        @test !isempty(String(take!(io)))
    end
end

@testset "flamegraph_export" begin
    @eval module M  # Take another tinf
        i(x) = x+5
        h(a::Array) = i(a[1]::Integer) + 2
        g(y::Integer) = h(Any[y])
    end

    tinf = @snoopi_deep begin
        M.g(2)
    end
    @test isempty(staleinstances(tinf))
    frames = flatten(tinf; sortby=inclusive)

    fg = SnoopCompile.flamegraph(tinf)
    @test length(collect(AbstractTrees.PreOrderDFS(fg))) == 5
    # Test that the span covers the whole tree.
    for leaf in AbstractTrees.PreOrderDFS(fg)
        @test leaf.data.span.start in fg.data.span
        @test leaf.data.span.stop in fg.data.span
    end

    frame1, frame2 = frames[1], frames[2]
    t1, t2 = inclusive(frame1), inclusive(frame2)
    # Ensure there's a tinf gap, and that cutting off the fastest-to-infer won't leave the tree headless
    if t1 != t2 && Method(frame1).name !== :g
        cutoff_bottom_frame = (t1 + t2) / 2
        fg2 = SnoopCompile.flamegraph(tinf, tmin = cutoff_bottom_frame)
        @test length(collect(AbstractTrees.PreOrderDFS(fg2))) == (length(collect(AbstractTrees.PreOrderDFS(fg))) - 1)
    end
    fg1 = flamegraph(tinf.children[1])
    @test endswith(string(fg.child.data.sf.func), "M.g") && endswith(string(fg1.child.data.sf.func), "M.h")
    fg2 = flamegraph(tinf.children[2])
    @test endswith(string(fg2.child.data.sf.func), "M.i")

    # Printing
    @eval module M
        i(x) = x+5
        h(a::Array) = i(a[1]::Integer) + 2
        g(y::Integer) = h(Any[y])
    end
    tinf = @snoopi_deep begin
        M.g(2)
        M.g(true)
    end
    @test isempty(staleinstances(tinf))
    fg = SnoopCompile.flamegraph(tinf)
    @test endswith(string(fg.child.data.sf.func), "M.g")
    counter = Dict{Method,Int}()
    visit(getfield(@__MODULE__, :M)) do item
        if isa(item, Core.MethodInstance)
            m = item.def
            if isa(m, Method)
                counter[m] = get(counter, m, 0) + 1
            end
            return false
        end
        return true
    end
    fg = SnoopCompile.flamegraph(tinf; mode=counter)
    @test endswith(string(fg.child.data.sf.func), "M.g (2)")
end

@testset "demos" begin
    # Just ensure they run
    @test SnoopCompile.itrigs_demo() isa SnoopCompile.InferenceTimingNode
    @test SnoopCompile.itrigs_higherorder_demo() isa SnoopCompile.InferenceTimingNode
end

include("testmodules/SnoopBench.jl")
@testset "parcel" begin
    a = SnoopBench.A()
    tinf = @snoopi_deep SnoopBench.f1(a)
    @test isempty(staleinstances(tinf))
    ttot, prs = SnoopCompile.parcel(tinf)
    mod, (tmod, tmis) = only(prs)
    @test mod === SnoopBench
    t, mi = only(tmis)
    @test ttot == tmod == t  # since there is only one
    @test mi.def.name === :f1
    ttot2, prs = SnoopCompile.parcel(tinf; tmin=10.0)
    @test isempty(prs)
    @test ttot2 == ttot

    A = [a]
    tinf = @snoopi_deep SnoopBench.mappushes(identity, A)
    @test isempty(staleinstances(tinf))
    ttot, prs = SnoopCompile.parcel(tinf)
    mod, (tmod, tmis) = only(prs)
    @test mod === SnoopBench
    @test ttot == tmod  # since there is only one
    @test length(tmis) == 2
    io = IOBuffer()
    SnoopCompile.write(io, tmis; tmin=0.0)
    str = String(take!(io))
    @test occursin(r"typeof\(mappushes\),Any,Vector\{A\}", str)
    @test occursin(r"typeof\(mappushes!\),typeof\(identity\),Vector\{Any\},Vector\{A\}", str)

    list = Any[1, 1.0, Float16(1.0), a]
    tinf = @snoopi_deep SnoopBench.mappushes(isequal(Int8(1)), list)
    @test isempty(staleinstances(tinf))
    ttot, prs = SnoopCompile.parcel(tinf)
    @test length(prs) == 2
    _, (tmodBase, tmis) = prs[findfirst(pr->pr.first === Base, prs)]
    tw, nw = SnoopCompile.write(io, tmis; tmin=0.0)
    @test 0.0 <= tw <= tmodBase * (1+10*eps())
    @test 0 <= nw <= length(tmis)
    str = String(take!(io))
    @test !occursin(r"Base.Fix2\{typeof\(isequal\).*SnoopBench.A\}", str)
    @test length(split(chomp(str), '\n')) == nw
    _, (tmodBench, tmis) = prs[findfirst(pr->pr.first === SnoopBench, prs)]
    @test sum(inclusive, tinf.children[1:end-1]) <= tmodBench + tmodBase # last child is not precompilable
    tw, nw = SnoopCompile.write(io, tmis; tmin=0.0)
    @test nw == 2
    str = String(take!(io))
    @test occursin(r"typeof\(mappushes\),Any,Vector\{Any\}", str)
    @test occursin(r"typeof\(mappushes!\),Base.Fix2\{typeof\(isequal\).*\},Vector\{Any\},Vector\{Any\}", str)

    td = joinpath(tempdir(), randstring(8))
    SnoopCompile.write(td, prs; tmin=0.0, ioreport=io)
    str = String(take!(io))
    @test occursin(r"Base: precompiled [\d\.]+ out of [\d\.]+", str)
    @test occursin(r"SnoopBench: precompiled [\d\.]+ out of [\d\.]+", str)
    file_base = joinpath(td, "precompile_Base.jl")
    @test isfile(file_base)
    @test occursin("ccall(:jl_generating_output", read(file_base, String))
    rm(td, recursive=true, force=true)
    SnoopCompile.write(td, prs; ioreport=io, header=false)
    str = String(take!(io))  # just to clear it in case we use it again
    @test !occursin("ccall(:jl_generating_output", read(file_base, String))
    rm(td, recursive=true, force=true)

    # issue #197
    f197(::Vector{T}) where T<:Integer = zero(T)
    g197(@nospecialize(x::Vector{<:Number})) = f197(x)
    g197([1,2,3])
    @test SnoopCompile.get_reprs([(rand(), mi) for mi in methodinstances(f197)])[1] isa AbstractSet
end

@testset "Specialization" begin
    Ts = subtypes(Any)
    tinf_unspec = @snoopi_deep SnoopBench.mappushes(SnoopBench.spell_unspec, Ts)
    tf_unspec = flatten(tinf_unspec)
    # To ensure independent data, invalidate all compiled CodeInstances
    mis = map(last, accumulate_by_source(MethodInstance, tf_unspec))
    for mi in mis
        SnoopCompile.isROOT(mi) && continue
        visit(mi) do item
            isa(item, Core.CodeInstance) || return true
            item.max_world = 0
            return true
        end
    end
    tinf_spec = @snoopi_deep SnoopBench.mappushes(SnoopBench.spell_spec, Ts)
    tf_spec = flatten(tinf_spec)
    @test length(tf_unspec) < length(Ts) ÷ 5
    @test any(tmi -> occursin("spell_unspec(::Any)", repr(MethodInstance(tmi))), tf_unspec)
    @test length(tf_spec) >= length(Ts)
    @test !any(tmi -> occursin("spell_spec(::Any)", repr(MethodInstance(tmi))), tf_unspec)
    @test !any(tmi -> occursin("spell_spec(::Type)", repr(MethodInstance(tmi))), tf_unspec)

    # fig, axs = plt.subplots(1, 2)

    nruns = 10^3
    SnoopBench.mappushes(SnoopBench.spell_spec, Ts)
    @profile for i = 1:nruns
        SnoopBench.mappushes(SnoopBench.spell_spec, Ts)
    end
    rit = runtime_inferencetime(tinf_spec)
    @test !any(rit) do (ml, _)
        endswith(string(ml.file), ".c")   # attribute all costs to Julia code, not C code
    end
    m = @which SnoopBench.spell_spec(first(Ts))
    dspec = rit[findfirst(pr -> pr.first == m, rit)].second
    @test dspec.tinf > dspec.trun      # more time is spent on inference than on runtime
    @test dspec.nspec >= length(Ts)
    # Check that much of the time in `mappushes!` is spend on runtime dispatch
    mp = @which SnoopBench.mappushes!(SnoopBench.spell_spec, [], first(Ts))
    dmp = rit[findfirst(pr -> pr.first == mp, rit)].second
    @test dmp.trtd >= 0.5*dmp.trun
    # pgdsgui(axs[1], rit; bystr="Inclusive", consts=true, interactive=false)

    Profile.clear()
    SnoopBench.mappushes(SnoopBench.spell_unspec, Ts)
    @profile for i = 1:nruns
        SnoopBench.mappushes(SnoopBench.spell_unspec, Ts)
    end
    rit = runtime_inferencetime(tinf_unspec)
    m = @which SnoopBench.spell_unspec(first(Ts))
    dunspec = rit[findfirst(pr -> pr.first == m, rit)].second # trunspec, trdtunspec, tiunspec, nunspec
    @test dunspec.tinf < dspec.tinf/10
    @test dunspec.trun < 10*dspec.trun
    @test dunspec.nspec == 1
    # Test that no runtime dispatch occurs in mappushes!
    dmp = rit[findfirst(pr -> pr.first == mp, rit)].second
    @test dmp.trtd == 0
    # pgdsgui(axs[2], rit; bystr="Inclusive", consts=true, interactive=false)
end


text:
module Tests

using JuMP
using ParameterJuMP
using Test

import SparseArrays

macro test_expression(expr)
    return esc(quote
        @test JuMP.isequal_canonical(@expression(m, $expr), $expr)
    end)
end

macro test_expression_with_string(expr, str)
    return esc(quote
        @test string(@inferred $expr) == $str
        @test_expression $expr
    end)
end

function test_basic(args...)
    m_slave = Model(args...)

    @variable(m_slave, x[1:2] == 4.0, Param())
    @test x == all_parameters(m_slave)
    @variable(m_slave, y[1:6])

    @constraint(m_slave, ctr1, 3*y[1] >= 2 - 7*x[1])

    @objective(m_slave, Min, 5*y[1])

    JuMP.optimize!(m_slave)

    @test 5/3 ≈ JuMP.dual(ctr1) atol=1e-3
    @test [-35/3, 0.0] ≈ JuMP.dual.(x) atol=1e-3
    @test [-26/3, 0.0, 0.0, 0.0, 0.0, 0.0] ≈ JuMP.value.(y) atol=1e-3
    @test -130/3 ≈ JuMP.objective_value(m_slave) atol=1e-3
    @test parametrized_dual_objective_value(m_slave) ≈ -35/3 * x[1] + 10/3 atol=1e-3
end

function test_multiple_parameters(args...)
    m_slave = Model(args...)

    @variable(m_slave, x[1:6] == 4.0, Param())
    @variable(m_slave, y[1:6])

    @constraint(m_slave, ctr1, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr2, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr3, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr4, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr5, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr6, 3*y[1] >= 2 - 7*x[3])
    @constraint(m_slave, ctr7, sum(3*y[i]+x[i] for i in 2:4) >= 2 - 7*x[3])
    @constraint(m_slave, ctr8, sum(3*y[i]+7.0*x[i]-x[i] for i in 2:4) >= 2 - 7*x[3])

    @objective(m_slave, Min, 5*y[1])

    JuMP.optimize!(m_slave)

    @test 5/3 ≈ JuMP.dual(ctr1) + JuMP.dual(ctr2) + JuMP.dual(ctr3) + JuMP.dual(ctr4) + JuMP.dual(ctr5) + JuMP.dual(ctr6) atol=1e-3
    @test 0.0 ≈ JuMP.dual(ctr7) atol=1e-3
    @test 0.0 ≈ JuMP.dual(ctr8) atol=1e-3
    @test [0.0, 0.0, -35/3, 0.0, 0.0, 0.0] ≈ JuMP.JuMP.dual.(x) atol=1e-3
    @test [-26/3, 0.0, 0.0, 0.0, 0.0, 0.0] ≈ JuMP.JuMP.value.(y) atol=1e-3
    @test -130/3 ≈ JuMP.objective_value(m_slave) atol=1e-3
    @test parametrized_dual_objective_value(m_slave) ≈ -35/3 * x[3] + 10/3 atol=1e-3
end

function test_lessthan(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    @test value(α) == 1.0
    set_value(α, -1.0)
    @test value(α) == -1.0
    @variable(model, x)
    cref = @constraint(model, x ≤ α)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α - 2 atol=1e-3

    set_value(α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α + 4 atol=1e-3
end

function test_equalto(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x == α)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α - 2 atol=1e-3

    set_value(α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α + 4 atol=1e-3
end

function test_greaterthan(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x >= α)
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ 1α

    set_value(α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ 1α
end

function test_multiple_parameters2(args...)
    model = Model(args...)
    @variable(model, α[1:10] == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, x == sum(2 * α[i] for i in 1:10))
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 20.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α[3]) == 2.0
    @test parametrized_dual_objective_value(model) ≈ 2sum(α)
end

function test_mixing_params_and_vars_1(args...)
    model = Model(args...)
    @variable(model, α[1:5] == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, sum(x for i in 1:5) == sum(2 * α[i] for i in 1:5))
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == 1/5
    @test JuMP.dual(α[3]) == 2/5
    @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α)
end

function test_mixing_params_and_vars_2(args...)
    model = Model(args...)
    @variable(model, α[1:5] == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, 0.0 == sum(-x + 2 * α[i] for i in 1:5))
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == 1/5
    @test JuMP.dual(α[3]) == 2/5
    @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α)
end

function test_mixing_params_and_vars_3(args...)
    model = Model(args...)
    @variable(model, α[1:5] == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, 0.0 == sum(-x + 2.0 + 2 * α[i] for i in 1:5))
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 4.0
    @test JuMP.dual(cref) == 1/5
    @test JuMP.dual(α[3]) == 2/5
    @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α) + 2
end

function test_lazy_duals_error(args...)
    model_2 = Model(args...)
    ParameterJuMP.set_lazy_duals(model_2)
    ParameterJuMP.set_lazy_duals(model_2) # warn
    ParameterJuMP.set_not_lazy_duals(model_2)
    ParameterJuMP.set_not_lazy_duals(model_2) # warn
    @variable(model_2, y == 1.0, Param())
    @test_throws ErrorException ParameterJuMP.set_lazy_duals(model_2)

    model_3 = Model(args...)
    ParameterJuMP.set_lazy_duals(model_3)
    @variable(model_3, y == 1.0, Param())
    @test_throws ErrorException ParameterJuMP.set_not_lazy_duals(model_3)
    @test !ParameterJuMP.is_sync(model_3)
end

function test_add_after_solve(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x <= α)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α - 2

    @variable(model, b == -2.0, Param())
    @test all_parameters(model) == [α, b]
    cref = @constraint(model, x <= b)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -2.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 0.0
    @test JuMP.dual(b) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -b - 4
end

function test_set_coefficient(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, x <= 0.0)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 0.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 0.0
    @test parametrized_dual_objective_value(model) ≈ 0α

    set_coefficient(cref, α, 1.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ α - 2
end

function test_change_coefficient(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x >= α)
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ 1α

    ParameterJuMP.set_coefficient(cref, α, -2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -2.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 2.0
    @test parametrized_dual_objective_value(model) ≈ 2α
end

function test_set_coefficient_lazy(args...)
    model = Model(args...)
    ParameterJuMP.set_lazy_duals(model)
    @variable(model, α == 1.0, Param())
    @variable(model, x)
    cref = @constraint(model, x <= 0.0)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 0.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 0.0
    @test parametrized_dual_objective_value(model) ≈ 0α

    set_coefficient(cref, α, 1.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ α - 2
end

function test_set_coefficient_lazy2(args...)
    model = Model(args...)
    ParameterJuMP.set_lazy_duals(model)
    @variable(model, α == 1.0, Param())
    @variable(model, b == 0.0, Param())
    @variable(model, x)
    cref = @constraint(model, x <= b)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 0.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 0.0
    @test JuMP.dual(b) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -b

    set_coefficient(cref, α, 1.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == 1.0
    @test JuMP.dual(b) == -1.0
    @test parametrized_dual_objective_value(model) ≈ α - b - 2
end

function test_remove_parameter_constraint(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x >= α)
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ 1α

    ParameterJuMP.delete_from_constraint(cref, α)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 0.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 0.0
    @test parametrized_dual_objective_value(model) ≈ 0α
end

function test_remove_parameter_all_constraints(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x >= α)
    @objective(model, Min, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 1.0
    @test parametrized_dual_objective_value(model) ≈ 1α

    ParameterJuMP.delete_from_constraints(α)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 0.0
    @test JuMP.dual(cref) == 1.0
    @test JuMP.dual(α) == 0.0
    @test parametrized_dual_objective_value(model) ≈ 0α
end

function test_no_duals(args...)
    model = Model(args...)
    ParameterJuMP.set_no_duals(model)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    cref = @constraint(model, x == α)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test isnan(JuMP.dual(α))

    model_2 = Model(args...)
    @variable(model_2, y == 1.0, Param())
    @test_throws ErrorException ParameterJuMP.set_no_duals(model_2)
end

function test_ParameterizedAffExpr(args...)
    m = Model()
    @variable(m, x)
    @variable(m, y)
    a = @variable(m, variable_type = Param())
    @test name(a) == ""
    set_name(a, "a")
    @test name(a) == "a"
    b = @variable(m, variable_type = Param())
    set_name(b, "b")
    @test all_parameters(m) == [a, b]

    exp1 = x + y + a
    @test typeof(exp1) == ParameterJuMP.PAE{Float64}
    @test length(exp1.v.terms) == 2
    exp1 = exp1 + y
    @test length(exp1.v.terms) == 2

    @test iszero(ParameterJuMP.PAE{Float64}())
    @test iszero(zero(exp1))
    @test iszero(one(exp1) - one(ParameterJuMP.PAE{Float64}))
    @test iszero(SparseArrays.dropzeros((exp1 - copy(exp1))))

    empty_func(empty_arg) = 0
    exp2 = map_coefficients(Base.zero, exp1)
    @test iszero(SparseArrays.dropzeros(exp2))
    @test iszero(constant(exp2))

    @test isequal_canonical(exp1, copy(exp1))

    exp4 = exp1 - copy(exp1)
    @test iszero(SparseArrays.dropzeros(exp4))

    # var + num
    @test_expression_with_string 7.1 * x + 2.5 "7.1 x + 2.5"
    @test_expression_with_string 1.2 * y + 1.2 "1.2 y + 1.2"
    # par + num
    @test_expression_with_string 7.1 * a + 2.5 "7.1 a + 2.5"
    @test_expression_with_string 1.2 * b + 1.2 "1.2 b + 1.2"

    # par + par + num
    @test_expression_with_string b + a + 1.2 "b + a + 1.2"
    # par - par + num
    @test_expression_with_string b - a + 1.2 "b - a + 1.2"
    # par - (par - num)
    @test_expression_with_string b - (a - 1.2) "b - a + 1.2"
    # par - par - num
    @test_expression_with_string b - a - 1.2 "b - a - 1.2"
    # var + par + num
    @test_expression_with_string x + a + 1.2 "x + a + 1.2"
    # var + par + num * num
    @test_expression_with_string x + a + 1.2 * 2.0 "x + a + 2.4"
    # var + var + par + num
    @test_expression_with_string x + y + a + 1.2 "x + y + a + 1.2"
    # var + var - par + num
    @test_expression_with_string x + y - a + 1.2 "x + y - a + 1.2"
    # var + var - par + par + num
    @test_expression_with_string x + y - a + b + 1.2 "x + y + b - a + 1.2"
    # par + par - var + var + num
    @test_expression_with_string a + b - x + y + 1.2 "y - x + a + b + 1.2"

    exp5 = x + y
    @test_expression_with_string x + y "x + y"
    @test_expression_with_string convert(ParameterJuMP.PAE{Float64}, exp5) "x + y"
end

function test_add_ctr_alaternative(args...)
    model = Model(args...)
    @variable(model, α == 1.0, Param())
    set_value(α, -1.0)
    @variable(model, x)
    exp = x - α
    cref = @constraint(model, exp ≤ 0)
    @objective(model, Max, x)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α - 2
end

function test_macro(args...)
    m = Model()
    @variable(m, x)
    @variable(m, y)
    @variable(m, a, ParameterJuMP.Param())
    @variable(m, b, ParameterJuMP.Param())
    # TODO - escaping
    # @variable(m, b, Param())
    # TODO - anonimous
    # b = @variable(m, Param())
    # set_name(b, "b")
    @test_expression_with_string a + b - x + y + 1.2 "y - x + a + b + 1.2"
end

function test_deletion_constraint(args...)
    model = Model(args...)
    @variable(model, α == -1.0, Param())
    @variable(model, x)
    cref1 = @constraint(model, x ≤ α/2)
    cref2 = @constraint(model, x ≤ α)
    cref3 = @constraint(model, x ≤ 2α)
    @objective(model, Max, x)
    JuMP.delete(model, cref3)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -1.0
    @test JuMP.dual(cref1) == 0.0
    @test JuMP.dual(cref2) == -1.0
    @test JuMP.dual(α) == -1.0
    @test parametrized_dual_objective_value(model) ≈ -α - 2

    JuMP.delete(model, cref2)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -0.5
    @test JuMP.dual(cref1) == -1.0
    @test JuMP.dual(α) == -0.5
    @test parametrized_dual_objective_value(model) ≈ -0.5α - 1
end

function test_add_to_expression(args...)
    model = Model(args...)
    @variable(model, a == 1.0, Param())
    @variable(model, x)
    ex = @expression(model, x + a)
    @test isequal(ex, x + a)
    # Number
    add_to_expression!(ex, 1.0)
    @test isequal(ex, x + a + 1.0)
    # ParameterRef
    add_to_expression!(ex, a)
    @test isequal(ex, x + 2.0 * a + 1.0)
    # PAE
    add_to_expression!(ex, copy(ex))
    @test isequal(ex, 2.0 * x + 4.0 * a + 2.0)
    # Number * Number
    add_to_expression!(ex, 2.0, 3.0)
    @test isequal(ex, 2.0 * x + 4.0 * a + 8.0)
    # Number * Variable
    add_to_expression!(ex, 3.0, x)
    @test isequal(ex, 5.0 * x + 4.0 * a + 8.0)
    # Variable * Number
    add_to_expression!(ex, x, 2.0)
    @test isequal(ex, 7.0 * x + 4.0 * a + 8.0)
    # Number * ParameterRef
    add_to_expression!(ex, 3.0, a)
    @test isequal(ex, 7.0 * x + 7.0 * a + 8.0)
    # ParameterRef * Number
    add_to_expression!(ex, a, 2.0)
    @test isequal(ex, 7.0 * x + 9.0 * a + 8.0)
    # PAE * Number
    add_to_expression!(ex, copy(ex), 2.0)
    @test isequal(ex, 3 * (7.0 * x + 9.0 * a + 8.0))
    # Number * PAE
    add_to_expression!(ex, 2.0, copy(ex))
    @test isequal(ex, 9 * (7.0 * x + 9.0 * a + 8.0))
end

function test_mutable_operate(args...)
    model = Model()
    x = @variable(model)
    @variable(model, p == 1.0, Param())
    exv = @expression(model, x + 1.0)
    exp = @expression(model, 2.0 * p)
    # -(PAE, Number, Number, Parameter)
    ex = @expression(model, x - 1 * (1 * p))
    @test isequal(ex, x - p)
    # -(PAE, Number, Number, Variable)
    ex = @expression(model, p - 1 * (1 * x))
    @test isequal(ex, p - x)
    # -(PAE, Number, Number, GAEv)
    ex = @expression(model, p - 1 * (1 * exv))
    @test isequal(ex, p - exv)
    # -(PAE, Number, Number, GAEp)
    ex = @expression(model, x - 1 * (1 * exp))
    a = x - exp; a.p.constant = abs(a.p.constant); a.v.constant = abs(a.v.constant)
    @test isequal(ex, a)

    # -(PAE, Number, Parameter, Number)
    ex = @expression(model, x - 1 * (p * 1))
    @test isequal(ex, x - p)
    # -(PAE, Number, Variable, Number)
    ex = @expression(model, p - 1 * (x * 1))
    @test isequal(ex, p - x)
    # -(PAE, Number, GAEv, Number)
    ex = @expression(model, p - 1 * (exv * 1))
    @test isequal(ex, p - exv)
    # -(PAE, Number, GAEp, Number)
    ex = @expression(model, x - 1 * (exp * 1))
    a = x - exp; a.p.constant = abs(a.p.constant); a.v.constant = abs(a.v.constant)
    @test isequal(ex, a)
end

function test_narg_mutable_operate(args...)
    model = Model()
    @variable(model, x[1:3])
    @variable(model, p[1:3] == 1.0, Param())
    @expression(model, ex[i=1:3], i * x[i] + 2.0 * p[i] + sqrt(i))
    @test isequal(
        @expression(model, sum(ex)),
        @expression(model, sum(i * x[i] + 2.0 * p[i] + sqrt(i) for i = 1:3))
    )
end

function test_float(args...)
    model = Model()
    @variable(model, p == 1, Param())
    @variable(model, x)
    @expression(model, ex, 2 * p + 1)
    @constraint(model, ex <= 0)
    @constraint(model, ex >= x)
    @constraint(model, ex <= 2x + 1)
end

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

end

ext: jl


text:
# Run this file from docs/
# by julia make.jl
using Documenter, Classic

makedocs(sitename="Classic.jl")


text:
function main()
    return 1,
            2
end

a,
b = main()
println(a)
println(b)

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using FunctionOperators, LinearAlgebra, Test, BenchmarkTools

data = [sin(i+j+k)^2 for i=1:300, j=1:300, k=1:50]
bOp₁ = FunctionOperator{Float64}(name="Op₁",
    forw = (buffer, x) -> buffer .= x.^2,
    backw = (buffer, x) -> broadcast!(sqrt, buffer, x),
    inDims = (300, 300, 50), outDims = (300, 300, 50))
weights = [sin((i-j)*l) + 1 for i=1:300, j=1:300, k=1:50, l=1:10]
bOp₂ = FunctionOperator{Float64}(name="Op₂",
    forw = (buffer,x) -> buffer .= reshape(x, 300, 300, 50, 1) .* weights,
    backw = (buffer,x) -> dropdims(sum!(reshape(buffer, 300, 300, 50, 1), x ./ weights), dims=4),
    inDims=(300, 300, 50), outDims=(300, 300, 50, 10))
combined = bOp₂ * (bOp₁ - 2.5*I) * bOp₁'
output = combined * data
function getAggregatedFunction()
    weights = [sin((i-j)*l) + 1 for i=1:300, j=1:300, k=1:50, l=1:10]
    buffer2 = Array{Float64}(undef, 300, 300, 50)
    buffer3 = Array{Float64}(undef, 300, 300, 50)
    buffer4 = Array{Float64}(undef, 300, 300, 50)
    (buffer, x) -> begin
        broadcast!(sqrt, buffer2, x)  # Of course, this two lines can be optimized to
        buffer3 .= buffer2 .^ 2       # (√x)^2 = |x|, but let's now avoid this fact
        broadcast!(-, buffer3, buffer3, broadcast!(*, buffer4, 2.5, buffer2))
        buffer .= reshape(buffer3, 300, 300, 50, 1) .* weights
    end
end
aggrFun = getAggregatedFunction()

@testset "Performance" begin
    t1 = @belapsed mul!(output, combined, data)
    t2 = @belapsed aggrFun(output, data)
    @test t1 / t2 ≈ 1 atol=0.05
end


text:
module WebIO

using Observables
using Requires
using AssetRegistry
using Base64: stringmime
import Widgets
import Widgets: node, AbstractWidget
using Logging
using UUIDs
using Base64
using RelocatableFolders

include("../deps/bundlepaths.jl")
include("../deps/mimetypes.jl")
"""
The filesystem path where the WebIO frontend packages lives.
"""
const packagepath = normpath(joinpath(@__DIR__, "..", "packages"))

"""
The path of the main WebIO JavaScript bundle file.
"""
const bundlepath = @path CORE_BUNDLE_PATH

"""
The MIME type for WebIO nodes.

This is used when serializing a WebIO node tree to a frontend.
For example, when using Jupyter, the frontend will render the serialized JSON
via this MIME bundle rather than using the data sent as `text/html`.
"""
const WEBIO_NODE_MIME = MIME"application/vnd.webio.node+json"
Base.Multimedia.istextmime(::WEBIO_NODE_MIME) = true

const WEBIO_APPLICATION_MIME = MIME"application/vnd.webio.application+html"
Base.Multimedia.istextmime(::WEBIO_APPLICATION_MIME) = true

include("util.jl")
include("connection.jl")
include("syntax.jl")
include("asset.jl")
include("node.jl")
include("iframe.jl")
include("observable.jl")
include("scope.jl")
include("render.jl")
include("messaging.jl")
include("rpc.jl")

# Extra "non-core" functionality
include("devsetup.jl")
include("../deps/jupyter.jl")

"""
    setup_provider(s::Union{Symbol, AbstractString})

Perform any side-effects necessary to set up the given provider. For example,
in IJulia, this causes the frontend to load the webio javascript bundle.
"""
setup_provider(s::Union{Symbol, AbstractString}) = setup_provider(Val(Symbol(s)))
export setup_provider

const baseurl = Ref{String}("")

"""
A union of all WebIO renderable types.
"""
const Application = Union{Scope, Node, AbstractWidget, Observable}
export Application

function setbaseurl!(str)
    baseurl[] = str
    if :ijulia in providers_initialised
        # re-run IJulia setup with the new base URL
        WebIO.setup_provider(Val{:ijulia}())
    end
end

const providers_initialised = Set{Symbol}()

function setup(provider::Symbol)
    haskey(ENV, "WEBIO_DEBUG") && println("WebIO: setting up $provider")
    haskey(ENV, "JULIA_WEBIO_BASEURL") && (baseurl[] = ENV["JULIA_WEBIO_BASEURL"])
    setup_provider(provider)
    push!(providers_initialised, provider)
    re_register_renderables()
end
setup(provider::AbstractString) = setup(Symbol(provider))

function prefetch_provider_file(basename)
  filepath = @path joinpath(@__DIR__, "providers", basename)
  code = read(filepath, String)
  (file = filepath, code = code)
end

provider_mux = prefetch_provider_file("mux.jl")
provider_blink = prefetch_provider_file("blink.jl")
provider_ijulia = prefetch_provider_file("ijulia.jl")
provider_generic_http = prefetch_provider_file("generic_http.jl")

function __init__()
    push!(Observables.addhandler_callbacks, WebIO.setup_comm)
    @require Mux="a975b10e-0019-58db-a62f-e48ff68538c9" begin
        include_string(@__MODULE__, provider_mux.code, provider_mux.file)
    end
    @require Blink="ad839575-38b3-5650-b840-f874b8c74a25" begin
        # The latest version of Blink defines their own WebIO integration
        # (after https://github.com/JunoLab/Blink.jl/pull/201).
        if isdefined(Blink.AtomShell, :initwebio!)
            return
        end
        Base.depwarn(
            "Please upgrade Blink for a smoother integration with WebIO.",
            :webio_blink_upgrade,
        )
        include_string(@__MODULE__, provider_blink.code, provider_blink.file)
    end
    @require IJulia="7073ff75-c697-5162-941a-fcdaad2a7d2a" begin
        include_string(@__MODULE__, provider_ijulia.code, provider_ijulia.file)
    end
    @require WebSockets="104b5d7c-a370-577a-8038-80a2059c5097" begin
        include_string(@__MODULE__, provider_generic_http.code, provider_generic_http.file)
    end

end

end # module


text:
"""
    submit_limit_order!(
        ob::OrderBook{Sz,Px,Oid,Aid},
        orderid::Oid,
        side::OrderSide,
        limit_price::Real,
        limit_size::Real,
        [, acct_id::Aid, fill_mode::OrderTraits ]
    )

Enter limit order with size `limit_size`, price `limit_price` with `side::OrderSide` into `ob::OrderBook`.

If an account if `acct_id` is provided, account holdings are tracked in `ob.acct_map`.

Order execution logic can be modified according to the argument `fill_mode::`[`OrderTraits`](@ref) which
defaults to `fill_mode=VANILLA_FILLTYPE`, representing the default order matching mode.

`submit_limit_order!` returns tuple of
 - `new_open_order::Order` representing the order left in the book after matching. Is `nothing` if no order was inserted
 - `order_match_lst::Vector{Order}` representing the matched orders if the order crosses the book.
 - `left_to_trade::Sz` representing the size of the portion of the order which could neither inserted nor matched.

"""
function submit_limit_order!(
    ob::OrderBook{Sz,Px,Oid,Aid},
    orderid::Oid,
    side::OrderSide,
    limit_price::Real,
    limit_size::Real,
    acct_id::Union{Nothing,Aid}=nothing,
    fill_mode::OrderTraits=VANILLA_FILLTYPE,
) where {Sz,Px,Oid,Aid}
    # Part 0 - Check Arguments
    if !((limit_price > zero(limit_price)) && (limit_size > zero(limit_size)))
        error("Both limit_price and limit_size must be positive")
    end
    ## Part 1 - Cross the Order if appropriate
    best_bid, best_ask = best_bid_ask(ob) # check that price in right range
    if allows_cross(fill_mode) &&
       isbuy(side) &&
       !isnothing(best_ask) &&
       (limit_price >= best_ask) # order is bid (buy) and can cross
        cross_match_lst, remaining_size = _walk_order_book_bysize!(
            ob.ask_orders, Sz(limit_size), Px(limit_price), fill_mode
        )
    elseif allows_cross(fill_mode) &&
           !isbuy(side) &&
           !isnothing(best_ask) &&
           (limit_price >= best_ask) # order is ask (sell) and can cross
        cross_match_lst, remaining_size = _walk_order_book_bysize!(
            ob.bid_orders, Sz(limit_size), Px(limit_price), fill_mode
        )
    else # order can or does not cross, return empty matches
        cross_match_lst, remaining_size = Vector{Order{Sz,Px,Oid,Aid}}(), Sz(limit_size)
    end

    ## Part 2 - Rest the remaining order in the book if possible
    if allows_book_insert(fill_mode) && !iszero(remaining_size) # if there are remaining shares, try to add remaining to the LOB
        best_bid, best_ask = best_bid_ask(ob) # new best bid and ask
        if isbuy(side) && (isnothing(best_ask) || (limit_price < best_ask)) # if order is a buy and limit price is valid for resting order
            # create and insert order object into BID book
            new_open_order = Order{Sz,Px,Oid,Aid}(
                side, remaining_size, limit_price, orderid, acct_id
            )
            insert_order!(ob.bid_orders, new_open_order)
            # if account_id present, add order account map
            isnothing(acct_id) || _add_order_acct_map!(ob.acct_map, acct_id, new_open_order)
            # set remaining size to zero
            remaining_size = zero(Sz)
        elseif !isbuy(side) && (isnothing(best_bid) || (limit_price > best_bid)) # if order is a sell and limit price is valid for resting order
            # create and insert order object into ASK book
            new_open_order = Order{Sz,Px,Oid,Aid}(
                side, remaining_size, limit_price, orderid, acct_id
            )
            insert_order!(ob.ask_orders, new_open_order)
            # if account_id present, add order account map
            !isnothing(acct_id) &&
                _add_order_acct_map!(ob.acct_map, acct_id, new_open_order)
            # set remaining size to zero
            remaining_size = zero(Sz)
        else # if not appropriate to insert new order, register new order as nothing
            new_open_order = nothing
        end
    else
        new_open_order = nothing
    end

    ## Part 3 - Return information
    return (
        new_open_order, cross_match_lst, remaining_size
    )::Tuple{Union{Order{Sz,Px,Oid,Aid},Nothing},Vector{Order{Sz,Px,Oid,Aid}},Sz}
end

@inline _is_best_price_inside_limit(::OneSidedBook, ::Nothing) = true

@inline function _is_best_price_inside_limit(
    sb::OneSidedBook{Sz,Px,Oid,Aid}, limit_price::Px
) where {Sz,Px,Oid,Aid}
    if isbidbook(sb)
        return sb.best_price >= limit_price
    else
        return sb.best_price <= limit_price
    end
end

"""
    _walk_order_book_bysize!(
        sb::OneSidedBook{Sz,Px,Oid,Aid},
        order_size::Sz,
        limit_price::Union{Px,Nothing},
        order_mode::OrderTraits,
    )

Cross (limit or market) order with opposite single side of book.
Order size is specified in number of shares.
Return matches and outstanding size.

__Notes__
 - If `limit_price::Nothing`, order is treated as Market Order
 - Function expects `order_size>0` and `limit_price>0`

"""
function _walk_order_book_bysize!(
    sb::OneSidedBook{Sz,Px,Oid,Aid},
    order_size::Sz,
    limit_price::Union{Px,Nothing},
    order_mode::OrderTraits,
)::Tuple{Vector{Order{Sz,Px,Oid,Aid}},Sz} where {Sz,Px,Oid,Aid}
    # Allocate memory for order output
    order_match_lst = Vector{Order{Sz,Px,Oid,Aid}}()
    shares_left = order_size # remaining quantity to trade
    # Perform initial available liquidity check
    if isallornone(order_mode) && (_size_available(sb, limit_price) < order_size)
        return order_match_lst, shares_left
    end
    # Perform matching logic
    limit_price_check = Base.Fix2(_is_best_price_inside_limit, limit_price)
    while !isempty(sb.book) && !iszero(shares_left) && limit_price_check(sb) # while book not empty, order not done and best price within limit price
        price_queue::OrderQueue = _popfirst_queue!(sb)
        if price_queue.total_volume[] <= shares_left # If entire queue is to be wiped out
            append!(order_match_lst, price_queue.queue) # Add all of the orders to the match list
            shares_left -= price_queue.total_volume[] # decrement what's left to trade
        else
            while !isempty(price_queue) && !iszero(shares_left) # while not done and queue not empty
                best_ord::Order = popfirst!(price_queue) # pop out best order
                if shares_left >= best_ord.size # Case 1: Limit order gets wiped out
                    # Add best_order to match list & decrement outstanding MO
                    push!(order_match_lst, best_ord)
                    shares_left -= best_ord.size
                else
                    # shares_left < best_ord.size # Case 2: Market Order gets wiped out
                    # Return the difference: LO.size-MO.size back to the order book
                    return_ord = copy_modify_size(best_ord, best_ord.size - shares_left)
                    pushfirst!(price_queue, return_ord)
                    # Add remainder to match list & decrement outstanding MO
                    best_ord = copy_modify_size(best_ord, shares_left)
                    push!(order_match_lst, best_ord)
                    shares_left -= best_ord.size
                end
            end
            if !isempty(price_queue) # If price queue wasn't killed, put it back into the OneSidedBook
                _insert_queue!(sb, price_queue)
            end
        end
    end
    # Return results
    return order_match_lst, shares_left
end

"""
Fill market order by walking the book using Price/Arrival priority.
__This version walks the book with funds rather than trade size.__
"""
function _walk_order_book_byfunds!(
    sb::OneSidedBook{Sz,Px,Oid,Aid},
    order_funds::Real,
    limit_price::Union{Px,Nothing},
    order_mode::OrderTraits,
)::Tuple{Vector{Order{Sz,Px,Oid,Aid}},Real} where {Sz,Px,Oid,Aid}
    # Allocate memory for order output
    order_match_lst = Vector{Order{Sz,Px,Oid,Aid}}()
    funds_left = order_funds # remaining quantity to trade
    # Perform initial available liquidity check
    if isallornone(order_mode) && (_size_available(sb, limit_price) < order_size)
        return order_match_lst, funds_left
    end
    # Perform matching logic
    limit_price_check = Base.Fix2(_is_best_price_inside_limit, limit_price)
    while !isempty(sb.book) && !iszero(funds_left) && limit_price_check(sb) # while book not empty, order not done and best price within limit price
        price_queue::OrderQueue = popfirst_queue!(sb)
        if (price_queue.total_volume[] * price_queue.price) <= funds_left # If entire queue is to be wiped out
            append!(order_match_lst, price_queue.queue) # Add all of the orders to the match list
            funds_left -= price_queue.total_volume[] * price_queue.price # decrement what's left to trade
        else
            while !isempty(price_queue) && !iszero(funds_left) # while not done and queue not empty
                best_ord::Order = popfirst!(price_queue) # pop out best order
                if funds_left >= (best_ord.size * best_ord.price) # Case 1: Limit order gets wiped out
                    # Add best_order to match list & decrement outstanding MO
                    push!(order_match_lst, best_ord)
                    funds_left -= (best_ord.size * best_ord.price)
                else # funds_left < best_ord.size # Case 2: Market Order gets wiped out
                    # Return the difference: LO.size-MO.size back to the order book
                    rem_match_size = floor(Sz, funds_left / best_ord.price)
                    return_ord = copy_modify_size(best_ord, best_ord.size - rem_match_size)
                    pushfirst!(price_queue, return_ord)
                    # Add remainder to match list & decrement outstanding MO
                    best_ord = copy_modify_size(best_ord, rem_match_size)
                    push!(order_match_lst, best_ord)
                    funds_left -= (best_ord.size * best_ord.price)
                end
            end
            if !isempty(price_queue) # If price queue wasn't killed, put it back into the OneSidedBook
                _insert_queue!(sb, price_queue)
            end
        end
    end
    # Return results
    return order_match_lst, funds_left
end

"""
    submit_market_order!(ob::OrderBook,side::OrderSide,mo_size[,fill_mode::OrderTraits])

Submit market order to `ob::OrderBook` with `side::OrderSide` and size `mo_size`.
Optionally `fill_mode::OrderTraits` may be provided to modify fill logic.
Market orders are filled by price-time priority.

Returns tuple `( ord_lst::Vector{Order}, left_to_trade::Sz )`
where
 - `ord_lst` is a list of _limit orders_ that _market order_ matched with
 - `left_to_trade` is the remaining size of un-filled order ( `==0` if order is complete, `>0` if incomplete)

__Note:__ Only `fill_mode.allornone` will be used from `fill_mode::OrderTraits`.
All other entries will be ignored.

"""
function submit_market_order!(
    ob::OrderBook{Sz,Px,Oid,Aid},
    side::OrderSide,
    mo_size::Real;
    fill_mode::OrderTraits=VANILLA_FILLTYPE,
) where {Sz,Px,Oid,Aid}
    mo_size > zero(mo_size) || error("market order argument mo_size must be positive")
    if isbuy(side)
        return _walk_order_book_bysize!(ob.ask_orders, Sz(mo_size), nothing, fill_mode)
    else
        return _walk_order_book_bysize!(ob.bid_orders, Sz(mo_size), nothing, fill_mode)
    end
end

"""
    submit_market_order_byfunds!(ob::OrderBook,side::Symbol,funds[,mode::OrderTraits])

Submit market order to `ob::OrderBook` `side::OrderSide` and available funds `funds::Real`.
Optionally `fill_mode::OrderTraits` may be provided to modify fill logic.
Market orders are filled by price-time priority.

Functionality is exactly the same as `submit_market_order!` except _available funds_ (max total price paid on order)
is provided, rather than _number of shares_ (order size).

Returns tuple `( ord_lst::Vector{Order}, funds_leftover )`
where
 - `ord_lst` is a list of _limit orders_ that _market order_ matched with
 - `funds_leftover` is the amount of remaining funds if not enough liquidity was available ( `==0` if order is complete, `>0` if incomplete)

__Note:__ Only `fill_mode.allornone` will be considered from `fill_mode::OrderTraits`.
All other entries will be ignored.

"""
function submit_market_order_byfunds!(
    ob::OrderBook, side::OrderSide, funds::Real, fill_mode::OrderTraits=VANILLA_FILLTYPE
)
    funds > zero(funds) || error("market order argument funds must be positive")
    if isbuy(side)
        return _walk_order_book_byfunds!(ob.ask_orders, funds, nothing, fill_mode)
    else
        return _walk_order_book_byfunds!(ob.bid_orders, funds, nothing, fill_mode)
    end
end

# Order Cancellation functions

"""
    cancel_order!(ob::OrderBook, o::Order)
    cancel_order!(ob::OrderBook, orderid, side, price [, acct_id=nothing])

Cancels Order `o`, or order with matching information from OrderBook.

Provide `acct_id` if known to guarantee correct account tracking.
"""
function cancel_order!(
    ob::OrderBook{Sz,Px,Oid,Aid}, orderid::Oid, side::OrderSide, price
) where {Sz,Px,Oid,Aid}
    # Delete order from bid (buy) or ask (sell) book
    if isbuy(side)
        popped_ord = pop_order!(ob.bid_orders, Px(price), orderid)
    else
        popped_ord = pop_order!(ob.ask_orders, Px(price), orderid)
    end
    # Delete order from account maps
    if !isnothing(popped_ord) && !isnothing(popped_ord.acctid)
        _delete_order_acct_map!(ob.acct_map, popped_ord.acctid, popped_ord.orderid)
    end
    return popped_ord
end

cancel_order!(ob::OrderBook, o::Order) = cancel_order!(ob, o.orderid, o.side, o.price)


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #
function _build_call_contract_object(data::Dict{String,Any})

    # grab stuff from data -
    # "symbol":"AMD",
    # "type":"call",
    # "sense":"sell",
    # "strike_price":60.0,
    # "premium_value":1.50,
    # "number_of_contracts":1,
    # "expiration":"2020-08-07"

    symbol = data["ticker_symbol"]
    sense = data["sense"]
    strike_price = data["strike_price"]
    premium_value = data["premium_value"]
    number_of_contracts = data["number_of_contracts"]
    expiration_date_string = data["expiration"]
    contract_multuplier = data["contract_multiplier"]

    # TODO - check sense, is this legit?
    # TODO - check is date string formatting correct?

    # build -
    callOptionContract = PSCallOptionContract(symbol, Date(expiration_date_string), strike_price, premium_value, number_of_contracts;
        sense=Symbol(sense), contractMultiplier=contract_multuplier)

    # return -
    return VLResult(callOptionContract)
end

function _build_put_contract_object(data::Dict{String,Any})

    symbol = data["ticker_symbol"]
    sense = data["sense"]
    strike_price = data["strike_price"]
    premium_value = data["premium_value"]
    number_of_contracts = data["number_of_contracts"]
    expiration_date_string = data["expiration"]
    contract_multiplier = data["contract_multiplier"]

    # TODO - check sense, is this legit?
    # TODO - check is date string formatting correct?

    # build -
    putOptionContract = PSPutOptionContract(symbol, Date(expiration_date_string), strike_price, premium_value, number_of_contracts;
        sense=Symbol(sense),contractMultiplier=contract_multiplier)

    # return -
    return VLResult(putOptionContract)
end

function _build_equity_object(data::Dict{String,Any})

    try

        # get stuff from dictionary -
        ticker_symbol = data["ticker_symbol"]
        purchase_price = data["purchase_price_per_share"]
        number_of_shares = data["number_of_shares"]
        purchase_date_string = data["purchase_date"]

        # build -
        equityObject = PSEquityAsset(ticker_symbol, purchase_price, number_of_shares, Date(purchase_date_string))

        # return -
        return VLResult(equityObject)
    catch error
        return VLResult(error)
    end
end
# ----------------------------------------------------------------------------------------------------------- #

# --- PUBLIC METHODS ---------------------------------------------------------------------------------------- #
function build_simulation_contract_set(simulation_dictionary::Dict{String,Any})::VLResult

    try

        # initialize -
        asset_set = Set{PSAbstractAsset}()

        # TODO: check - do we have the correct keys?

        # grab the list of asset dictionaries -
        asset_dictionary_array = simulation_dictionary["contract_set_parameters"]
        for (_, asset_dictionary) in enumerate(asset_dictionary_array)

            # initialize -
            local result = nothing

            # ok, so lets grab data from the asset_dictionary, and build each asset type -
            type_string = asset_dictionary["type"]
            type_symbol = Symbol(type_string)
            if (type_symbol == :call)
                result = _build_call_contract_object(asset_dictionary) |> check
            elseif (type_symbol == :put)
                result = _build_put_contract_object(asset_dictionary) |> check
            elseif (type_symbol == :equity)
                result = _build_equity_object(asset_dictionary) |> check
            end

            # grab -
            push!(asset_set, result)
        end

        # return -
        return VLResult(asset_set)

    catch error
        return VLResult(error)
    end
end

function build_simulation_contract_set(pathToSimulationFile::String)::VLResult

    try

        # load the experimet file -
        simulation_dictionary = JSON.parsefile(pathToSimulationFile)

        # build the contract set -
        return build_simulation_contract_set(simulation_dictionary)

    catch error
        return VLResult(error)
    end
end
# ----------------------------------------------------------------------------------------------------------- #


text:
using RDatasets, DataFrames, InteractiveFixedEffectModels, Distances, Base.Test

precision = 2e-1

for method in [:levenberg_marquardt, :dogleg, :gauss_seidel]
	println(method)

	df = dataset("plm", "Cigar")
	df[:pState] = pool(df[:State])
	df[:pYear] = pool(df[:Year])

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 1), Sales ~ Price, df, method = method, save = true)
	@test norm(result.coef ./ [328.1653237715761, -1.0415042260420706] .- 1)  < precision
	@test norm(abs(result.augmentdf[1, :factors1]) / 0.228 - 1) < precision
	@test norm(abs(result.augmentdf[1, :loadings1]) / 851.514 - 1) < precision
	@test norm(result.augmentdf[1, :residuals] / -9.7870 - 1) < precision
	@test result.r2_within > 0.0

	show(result)
	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price, df, method = method, save = true)
	@test norm(result.coef ./ [163.01350, -0.40610] - 1) < precision
	@test norm(abs(result.augmentdf[1, :factors1]) / 0.227 - 1) < precision
	@test norm(abs(result.augmentdf[1, :loadings1]) / 184.1897 - 1) < precision
	@test norm(result.augmentdf[1, :residuals]  / -1.774 - 1) < precision
	@test result.r2_within > 0.0

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price |> pState, df, method = method, save = true)
	@test norm(result.coef / -0.425389 - 1) < precision
	@test norm(abs(result.augmentdf[1, :factors1])  / 0.2474 - 1) < precision
	@test norm(abs(result.augmentdf[1, :loadings1]) / 123.460 - 1) < precision
	@test norm(result.augmentdf[1, :residuals] / -5.222 - 1) < precision
	@test result.r2_within > 0.0

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price |> pYear, df, method = method, save = true)
	@test norm(result.coef / -0.3744296120563005 -1 ) < precision
	@test norm(abs(result.augmentdf[1, :factors1])  / 0.1918 - 1) < precision
	@test norm(abs(result.augmentdf[1, :loadings1]) / 102.336 - 1) < precision
	@test norm(result.augmentdf[1, :residuals] / -2.2153 - 1) < precision
	@test result.r2_within > 0.0

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price |> pState + pYear, df, method = method, save = true)
	@test norm(result.coef / -0.524157 - 1) < precision
	@test norm(abs(result.augmentdf[1, :factors1]) / 0.256 - 1) < precision
	@test norm(abs(result.augmentdf[1, :loadings1]) / 60.0481 - 1) < precision
	@test norm(result.augmentdf[1, :residuals] / -5.614 - 1) < precision
	@test result.r2_within > 0.0

	# subset
	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price |> pState, df, method = method, subset = (df[:State] .<= 30), save = true)
	@test size(result.augmentdf, 1) == size(df, 1)
	@test norm(abs(result.augmentdf[:factors1][1]) /0.25965 - 1) < precision
	@test norm(abs(result.augmentdf[:loadings1][1]) /107.832 - 1) < precision
	@test norm(abs(result.augmentdf[:factors2][1])  /0.2655  - 1) < precision
	@test norm(abs(result.augmentdf[:loadings2][1])    /15.551 - 1) < precision
	@test norm(result.augmentdf[:residuals][1] /-3.8624  - 1) < precision
	@test norm(result.augmentdf[:pState][1] /131.6162 - 1) < precision
end

# local minima
using RDatasets, DataFrames, InteractiveFixedEffectModels, Distances, Base.Test
precision = 2e-1

for method in [:levenberg_marquardt, :dogleg]
	println(method)

	df = dataset("plm", "Cigar")
	df[:pState] = pool(df[:State])
	df[:pYear] = pool(df[:Year])
	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 1), Sales ~ Price |> pState, df, method = method,
	weight = :Pop, save = true)
	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ Price |> pState, df, method = method, weight = :Pop, save = true)

	df = dataset("plm", "EmplUK")
	df[:id1] = df[:Firm]
	df[:id2] = df[:Year]
	df[:pid1] = pool(df[:id1])
	df[:pid2] = pool(df[:id2])
	df[:y] = df[:Wage]
	df[:x1] = df[:Emp]
	df[:w] = df[:Output]
	result = fit(InteractiveFixedEffectModel(:pid1, :pid2, 2), y ~ x1, df, method = method, save = true)
	@test norm(result.coef ./ [4.53965, -0.0160858] - 1) < precision
	result = fit(InteractiveFixedEffectModel(:pid1, :pid2, 2), y ~ x1, df, method = method, save = true, weight = :w)
	@test norm(result.coef ./ [3.47551,-0.017366] - 1) < precision
	result = fit(InteractiveFixedEffectModel(:pid1, :pid2, 1), y ~ x1, df, method = method, save = true, weight = :w)
	@test norm(result.coef ./ [ -2.62105, -0.0470005] - 1) < precision
end


text:
# This basically just defines an API and utilities for
# JS RPCs – the implementation is left to specific objects.

# Note that this API is very limited. Although it would
# be perfectly possible to have full interop with references
# to JS objects etc., it would be an enormously leaky
# abstraction over a TCP connection with any latency.

# Convert Julia `Expr`s to strings of JS

include("jsexprs.jl")
include("callbacks.jl")

# RPC API

export js, js_, @js, @js_, @var, @new

msg(o, m) = error("$(typeof(o)) object doesn't support JS messages")

function js(o, js::JSString; callback = true)
  cmd = @d(:type => :eval,
           :code => js.s)
  if callback
    id, cond = callback!()
    cmd[:callback] = id
  end
  msg(o, cmd)
  return callback ? wait(cond) : o
end

js(o, s; callback = true) = js(o, jsexpr(s); callback = callback)

js_(args...) = js(args..., callback = false)

macro js(o, ex)
  :(js($(esc(o)), $(Expr(:quote, ex))))
end

macro js_(o, ex)
  :(js_($(esc(o)), $(Expr(:quote, ex))))
end


text:
using Oceananigans.Grids: solid_interface

"""
    abstract type AbstractScalarDiffusivity <: AbstractTurbulenceClosure end

Abstract type for closures with *isotropic* diffusivities.
"""
abstract type AbstractScalarBiharmonicDiffusivity{F} <: AbstractTurbulenceClosure{ExplicitTimeDiscretization} end

@inline formulation(::AbstractScalarBiharmonicDiffusivity{F}) where {F} = F()

#####
##### Stress divergences
#####

const AIBD = AbstractScalarBiharmonicDiffusivity{<:ThreeDimensionalFormulation}
const AHBD = AbstractScalarBiharmonicDiffusivity{<:HorizontalFormulation}
const AVBD = AbstractScalarBiharmonicDiffusivity{<:VerticalFormulation}

@inline viscous_flux_ux(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), ∂xᶜᶜᶜ, biharmonic_mask_x, ∇²ᶠᶜᶜ, U.u)
@inline viscous_flux_vx(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_x, ∂xᶠᶠᶜ, ∇²ᶜᶠᶜ, U.v)
@inline viscous_flux_wx(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_x, ∂xᶠᶜᶠ, ∇²ᶜᶜᶠ, U.w)
@inline viscous_flux_uy(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_y, ∂yᶠᶠᶜ, ∇²ᶠᶜᶜ, U.u)
@inline viscous_flux_vy(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), ∂yᶜᶜᶜ, biharmonic_mask_y, ∇²ᶜᶠᶜ, U.v)
@inline viscous_flux_wy(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_y, ∂yᶜᶠᶠ, ∇²ᶜᶜᶠ, U.w)
@inline viscous_flux_uz(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_z, ∂zᶠᶜᶠ, ∇²ᶠᶜᶜ, U.u)
@inline viscous_flux_vz(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_z, ∂zᶜᶠᶠ, ∇²ᶜᶠᶜ, U.v)
@inline viscous_flux_wz(i, j, k, grid, closure::AIBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), ∂zᶜᶜᶜ, biharmonic_mask_z, ∇²ᶜᶜᶠ, U.w)

@inline viscous_flux_ux(i, j, k, grid, closure::AHBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), δ★ᶜᶜᶜ, U.u, U.v)
@inline viscous_flux_vx(i, j, k, grid, closure::AHBD, clock, U, args...) = + ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), ζ★ᶠᶠᶜ, U.u, U.v)
@inline viscous_flux_wx(i, j, k, grid, closure::AHBD, clock, U, args...) = + ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_x, ∂xᶠᶜᶠ, ∇²ᶜᶜᶠ, U.w)
@inline viscous_flux_uy(i, j, k, grid, closure::AHBD, clock, U, args...) = - ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), ζ★ᶠᶠᶜ, U.u, U.v)
@inline viscous_flux_vy(i, j, k, grid, closure::AHBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), δ★ᶜᶜᶜ, U.u, U.v)
@inline viscous_flux_wy(i, j, k, grid, closure::AHBD, clock, U, args...) = + ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_y, ∂yᶜᶠᶠ, ∇²ᶜᶜᶠ, U.w)

@inline viscous_flux_uz(i, j, k, grid, closure::AVBD, clock, U, args...) = + ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_z, ∂zᶠᶜᶠ, ∂²zᶠᶜᶜ, U.u)
@inline viscous_flux_vz(i, j, k, grid, closure::AVBD, clock, U, args...) = + ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), biharmonic_mask_z, ∂zᶜᶠᶠ, ∂²zᶜᶠᶜ, U.v)
@inline viscous_flux_wz(i, j, k, grid, closure::AVBD, clock, U, args...) = + ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), ∂zᶜᶜᶜ, biharmonic_mask_z, ∂²zᶜᶜᶠ, U.w)

#####
##### Diffusive fluxes
#####

@inline diffusive_flux_x(i, j, k, grid, closure::AIBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶠᶜᶜ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_x, ∂xᶠᶜᶜ, ∇²ᶜᶜᶜ, c)
@inline diffusive_flux_y(i, j, k, grid, closure::AIBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶜᶠᶜ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_y, ∂yᶜᶠᶜ, ∇²ᶜᶜᶜ, c)
@inline diffusive_flux_z(i, j, k, grid, closure::AIBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶜᶜᶠ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_z, ∂zᶜᶜᶠ, ∇²ᶜᶜᶜ, c)

@inline diffusive_flux_x(i, j, k, grid, closure::AHBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶠᶜᶜ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_x, ∂x_∇²h_cᶠᶜᶜ, c)
@inline diffusive_flux_y(i, j, k, grid, closure::AHBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶜᶠᶜ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_y, ∂y_∇²h_cᶜᶠᶜ, c)
@inline diffusive_flux_z(i, j, k, grid, closure::AVBD, c, ::Val{tracer_index}, clock, args...) where tracer_index = κ_σᶜᶜᶠ(i, j, k, grid, clock, closure.κ[tracer_index], biharmonic_mask_z, ∂³zᶜᶜᶠ, c)

#####
##### Zero out not used fluxes
#####

for (dir, closure) in zip((:x, :y, :z), (:AVBD, :AVBD, :AHBD))
    diffusive_flux = Symbol(:diffusive_flux_, dir)
    viscous_flux_u = Symbol(:viscous_flux_u, dir)
    viscous_flux_v = Symbol(:viscous_flux_v, dir)
    viscous_flux_w = Symbol(:viscous_flux_w, dir)
    @eval begin
        @inline $diffusive_flux(i, j, k, grid, closure::$closure, c, c_idx, clock, args...) = zero(eltype(grid))
        @inline $viscous_flux_u(i, j, k, grid, closure::$closure, c, clock, U, args...)     = zero(eltype(grid))
        @inline $viscous_flux_v(i, j, k, grid, closure::$closure, c, clock, U, args...)     = zero(eltype(grid))
        @inline $viscous_flux_w(i, j, k, grid, closure::$closure, c, clock, U, args...)     = zero(eltype(grid))
    end
end

#####
##### Biharmonic-specific viscous operators
#####

# See https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#horizontal-dissipation
@inline function δ★ᶜᶜᶜ(i, j, k, grid, u, v)

    # These closures seem to be needed to help the compiler infer types
    # (either of u and v or of the function arguments)
    @inline Δy_∇²u(i, j, k, grid, u) = Δy_qᶠᶜᶜ(i, j, k, grid, biharmonic_mask_x, ∇²hᶠᶜᶜ, u)
    @inline Δx_∇²v(i, j, k, grid, v) = Δx_qᶜᶠᶜ(i, j, k, grid, biharmonic_mask_y, ∇²hᶜᶠᶜ, v)

    return 1 / Azᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Δy_∇²u, u) +
                                       δyᵃᶜᵃ(i, j, k, grid, Δx_∇²v, v))
end

@inline function ζ★ᶠᶠᶜ(i, j, k, grid, u, v)

    # These closures seem to be needed to help the compiler infer types
    # (either of u and v or of the function arguments)
    @inline Δy_∇²v(i, j, k, grid, v) = Δy_qᶜᶠᶜ(i, j, k, grid, biharmonic_mask_y, ∇²hᶜᶠᶜ, v)
    @inline Δx_∇²u(i, j, k, grid, u) = Δx_qᶠᶜᶜ(i, j, k, grid, biharmonic_mask_x, ∇²hᶠᶜᶜ, u)

    return 1 / Azᶠᶠᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Δy_∇²v, v) -
                                       δyᵃᶠᵃ(i, j, k, grid, Δx_∇²u, u))
end

#####
##### Biharmonic-specific diffusion operators
#####

@inline ∂x_∇²h_cᶠᶜᶜ(i, j, k, grid, c) = 1 / Azᶠᶜᶜ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶜᶜ, ∇²hᶜᶜᶜ, c)
@inline ∂y_∇²h_cᶜᶠᶜ(i, j, k, grid, c) = 1 / Azᶜᶠᶜ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, Δx_qᶜᶜᶜ, ∇²hᶜᶜᶜ, c)

#####
##### Biharmonic-specific operators that enforce "no-flux" boundary conditions and "0-value" boundary conditions for the Laplacian operator
##### biharmonic_mask(∂(∇²(var))) ensures that fluxes are 0 on the boundaries
##### ∂(biharmonic_mask(∇²(var))) ensures that laplacians are 0 on the boundaries
#####

biharmonic_mask_x(i, j, k, grid, f, args...) = ifelse(solid_x(i, j, k, grid), zero(eltype(grid)), f(i, j, k, grid, args...))
biharmonic_mask_y(i, j, k, grid, f, args...) = ifelse(solid_y(i, j, k, grid), zero(eltype(grid)), f(i, j, k, grid, args...))
biharmonic_mask_z(i, j, k, grid, f, args...) = ifelse(solid_z(i, j, k, grid), zero(eltype(grid)), f(i, j, k, grid, args...))

solid_x(i, j, k, grid) = solid_interface(Face(), Center(), Center(), i, j, k, grid)
solid_y(i, j, k, grid) = solid_interface(Center(), Face(), Center(), i, j, k, grid)
solid_z(i, j, k, grid) = solid_interface(Center(), Center(), Face(), i, j, k, grid)


text:
include("../src/components/sealevelrise_component.jl")

@testset "sealevelrise" begin

    Precision = 1.0e-11
    T = length(MimiDICE2010.model_years)
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, sealevelrise, :sealevelrise)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :sealevelrise, :TATM, read_params(f, "B121:BI121", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :sealevelrise, :therm0, p[:therm0])
    set_param!(m, :sealevelrise, :gsic0, p[:gsic0])
    set_param!(m, :sealevelrise, :gis0, p[:gis0])
    set_param!(m, :sealevelrise, :ais0, p[:ais0])
    set_param!(m, :sealevelrise, :therm_asym, p[:therm_asym])
    set_param!(m, :sealevelrise, :gsic_asym, p[:gsic_asym])
    set_param!(m, :sealevelrise, :gis_asym, p[:gis_asym])
    set_param!(m, :sealevelrise, :ais_asym, p[:ais_asym])
    set_param!(m, :sealevelrise, :thermrate, p[:thermrate])
    set_param!(m, :sealevelrise, :gsicrate, p[:gsicrate])
    set_param!(m, :sealevelrise, :gisrate, p[:gisrate])
    set_param!(m, :sealevelrise, :aisrate, p[:aisrate])
    set_param!(m, :sealevelrise, :slrthreshold, p[:slrthreshold])

# Run the one-component model
    run(m)

# Extract the generated variables
    ThermSLR    = m[:sealevelrise, :ThermSLR]
    GSICSLR     = m[:sealevelrise, :GSICSLR]
    GISSLR      = m[:sealevelrise, :GISSLR]
    AISSLR      = m[:sealevelrise, :AISSLR]
    TotSLR      = m[:sealevelrise, :TotSLR]

# Extract the true values
    True_ThermSLR    = read_params(f, "B178:BI178", T)
    True_GSICSLR    = read_params(f, "B179:BI179", T)
    True_GISSLR    = read_params(f, "B180:BI180", T)
    True_AISSLR    = read_params(f, "B181:BI181", T)
    True_TotSLR    = read_params(f, "B182:BI182", T)

# Test that the values are the same
    @test maximum(abs, ThermSLR .- True_ThermSLR) ≈ 0. atol = Precision
    @test maximum(abs, GSICSLR .- True_GSICSLR) ≈ 0. atol = Precision
    @test maximum(abs, GISSLR .- True_GISSLR) ≈ 0. atol = Precision
    @test maximum(abs, AISSLR .- True_AISSLR) ≈ 0. atol = Precision
    @test maximum(abs, TotSLR .- True_TotSLR) ≈ 0. atol = Precision

end


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, van

testTol = 0.2

for prob in [prob_ode_linear, prob_ode_2Dlinear]
  sim21 = test_convergence(1 .//2 .^(6:-1:3),prob,RadauIIA5())
  @test sim21.𝒪est[:final] ≈ 5 atol=testTol
end

# test adaptivity
for iip in (true, false)
  if iip
    vanstiff = ODEProblem{iip}(van, [0;sqrt(3)], (0.0,1.0), 1e6)
  else
    vanstiff = ODEProblem{false}((u,p,t)->van(u,p,t), [0;sqrt(3)], (0.0,1.0), 1e6)
  end
  sol = solve(vanstiff, RadauIIA5())
  if iip
    @test sol.destats.naccept + sol.destats.nreject > sol.destats.njacs # J reuse
    @test sol.destats.njacs < sol.destats.nw # W reuse
  end
  @test length(sol) < 150
  @test length(solve(remake(vanstiff, p=1e7), RadauIIA5())) < 150
  @test length(solve(remake(vanstiff, p=1e7), reltol=[1e-4, 1e-6], RadauIIA5())) < 170
  @test length(solve(remake(vanstiff, p=1e7), RadauIIA5(), reltol=1e-9, abstol=1e-9)) < 870
  @test length(solve(remake(vanstiff, p=1e9), RadauIIA5())) < 170
  @test length(solve(remake(vanstiff, p=1e10), RadauIIA5())) < 190
end

##Tests fot RadauIIA3
for prob in [prob_ode_linear, prob_ode_2Dlinear]
  dts = 1 ./2 .^(8:-1:1)
  sim = test_convergence(dts, prob, RadauIIA3())
  @test sim.𝒪est[:final] ≈ 3 atol = 0.25
end

# test adaptivity
for iip in (true, false)
  if iip
    vanstiff = ODEProblem{iip}(van, [0;sqrt(3)], (0.0,1.0), 1e6)
  else
    vanstiff = ODEProblem{false}((u,p,t)->van(u,p,t), [0;sqrt(3)], (0.0,1.0), 1e6)
  end
  sol = solve(vanstiff, RadauIIA3())
  if iip
    @test sol.destats.naccept + sol.destats.nreject > sol.destats.njacs # J reuse
    @test sol.destats.njacs < sol.destats.nw # W reuse
  end
  @test length(sol) < 5000 # the error estimate is not very good
end


text:
@testset "Unconstrained LTI system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    C = Matrix([0.5 1.5])
    D = Matrix([1.0]')
    s = LinearTimeInvariantSystem(A, B, C, D)
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test outputmap(s) isa LinearControlMap

    # check alias
    s = LTISystem(A, B, C, D)

    # check the LTI system's outputmap fields
    @test outputmap(s).A == C
    @test outputmap(s).B == D

    # check the LTI system's system (s) field
    @test s.s.A == A
    @test s.s.B == B

    # check that the system is linear
    @test islinear(s) && isaffine(s)
end

@testset "Constrained LTI system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    C = Matrix([0.5 1.5])
    D = Matrix([1.0]')
    X = Line([1., -1], 0.)
    U = Interval(0.9, 1.1)
    s = LinearTimeInvariantSystem(A, B, C, D, X, U)
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test outputmap(s) isa ConstrainedLinearControlMap

    # check alias
    s = LTISystem(A, B, C, D, X, U)

    # check that the system is linear
    @test islinear(s) && isaffine(s)
end


text:
z = dfs_tree(g5,1)

@test ne(z) == 3 && nv(z) == 4
@test !has_edge(z, 1, 3)

@test topological_sort_by_dfs(g5) == [1, 2, 3, 4]
@test !is_cyclic(g5)
g = DiGraph(3)

add_edge!(g,1,2); add_edge!(g,2,3); add_edge!(g,3,1)

@test is_cyclic(g)


text:
using ConvOptDL

include("data_loader_tests.jl")
include("util_tests.jl")
include("qp_tests.jl")
include("adjoin_tests.jl")
include("drop_block_tests.jl")


text:
function update_r!(s::StateFS, c::ConstantsFS, d::DataFS;
                   marg_lam::Bool=true, joint_update::Bool=true)
  if joint_update  # update r jointly. Proposed r flipped by one bit.
    idx = CartesianIndices(s.r)
    idx_to_flip = Distributions.sample(idx)
    i, k = Tuple(idx_to_flip)
    if marg_lam
      update_r_marg_lam!(i, k, s, c, d)
    else
      update_r!(i, k, s, c, d)
    end
    # NOTE: Make sure to always update W immediately after updating r or
    #       W_star!
    update_W!(s, c, d)
  else
    # update each r_{i, k}
    for i in 1:d.data.I
      for k in 1:c.constants.K
        if marg_lam
          update_r_marg_lam!(i, k, s, c, d)
        else
          update_r!(i, k, s, c, d)
        end
        # NOTE: Make sure to always update W immediately after updating r or
        #       W_star!
        update_W!(s, c, d)
      end
    end
  end
end

function compute_w(ws::Vector{Float64}, r::Vector{Bool})
  wr = ws .* r
  return wr ./ sum(wr)
end

function logdmixture_r(i::Integer, n::Integer, k::Integer,
                       s::State, c::Constants, d::Data)::Float64
  return sum(logdmixture(s.Z[j, k], i, n, j, s, c, d) for j in 1:d.J)
end

function update_r_marg_lam!(i::Integer, k::Integer,
                            s::StateFS, c::ConstantsFS, d::DataFS)

  # Log likelihood as a function of `r_{i, k}`
  function loglike(r_ik::Bool)::Float64
    # Make a copy of the current (k-dim) vector r_i
    ri = s.r[i, :]

    # Replace r_{i, k} with the provided `r_ik`
    ri[k] = r_ik

    # Compute W given W* and the updated r_i
    wi = compute_w(s.W_star[i, :], ri)

    # Get the indices of the W_i such that W_{i, k} > 0
    non_zero_wri_idx = findall(wi .> 0)

    # Initialize log likelihood
    ll = 0.0

    if length(non_zero_wri_idx) > 0  # i.e. if any W_i > 0
      for n in 1:d.data.N[i]
        ldmix = [logdmixture_r(i, n, k_prime, s.theta, c.constants, d.data)
                 for k_prime in non_zero_wri_idx]

        # Add to log likelihood
        ll += MCMC.logsumexp(log.(wi[non_zero_wri_idx]) + ldmix)
      end
    else
      ll = -Inf
    end

    return ll
  end

  # Log prior as a function of `r_{i, k}`
  function logprior(r_ik::Bool)::Float64
    p_xi = compute_p(d.X[i, :], s.omega)
    return logpdf(Bernoulli(p_xi), r_ik)
  end

  # Log full conditional as a function of `r_{i, k}`
  log_prob(r_ik::Bool)::Float64 = logprior(r_ik) + loglike(r_ik)

  # update r_{i, k} with a metropolis step
  rik_metropolis_update!(i, k, log_prob, s)
end

function rik_metropolis_update!(i::Integer, k::Integer,
                                log_prob::Function, s::StateFS)
  cand = !s.r[i, k]  # Flip bit
  log_acceptance_ratio = log_prob(cand) - log_prob(s.r[i, k])

  # accept with probability `accecptance_ratio`
  if log_acceptance_ratio > log(rand())
    s.r[i, k] = cand
  end
end

log1m(x::T) where {T <: Real} = log1p(-x)
sf_lgamma(x::T) where T = Distributions.SpecialFunctions.lgamma(x)
# FIXME: spitting out NaN's for W!
#        pretty sure this is wrong!
function update_r!(i::Integer, k::Integer,
                   s::StateFS, c::ConstantsFS, d::DataFS; verbose::Int=0)
  p_xi = compute_p(d.X[i, :], s.omega)
  aw = shape(c.W_star_prior)
  aw_ri = aw * s.r[i, :]
  aw_ri_sum = sum(aw_ri)
  u = aw_ri_sum - aw_ri[k]
  log_numer = log1m(p_xi)
  if u > 0
    log_numer += sf_lgamma(u) + sf_lgamma(aw)
  else
    log_numer = -Inf
  end
  log_denom = log(p_xi) + sf_lgamma(aw + u)
  if s.theta.W[i, k] > 0
    log_denom += aw * log(s.theta.W[i, k])
  end
  if verbose > 0
    println("log_numer: $log_numer")
    println("log_denom: $log_denom")
  end
  prob = 1 / (1 + exp(log_numer - log_denom))
  s.r[i, k] = prob > rand()
end


text:
using Symplectic
using LinearAlgebra
using Optim
using Plots
using LaTeXStrings
using Plots.PlotMeasures

cr = CoupledResonators()

cr = addPassive(cr, 1.0, 1, 1)
cr = addPassive(cr, 1.0, 2, 2)
cr = addPassive(cr, 1.0, 3, 3)
cr = addPassive(cr, 1.0, 4, 4)
cr = addPassive(cr, 1.0, 5, 5)

cr = addPassive(cr, 1.0, 1, 2 )
cr = addActive(cr, 1.0, 1, 2 )

cr = addPassive(cr, 1.0, 1, 3 )
cr = addActive(cr, 1.0, 1, 3 )

cr = addPassive(cr, 1.0, 1, 4 )
cr = addActive(cr, 1.0, 1, 4 )

cr = addPassive(cr, 1.0, 5, 2 )
cr = addActive(cr, 1.0, 5, 2 )

cr = addPassive(cr, 1.0, 5, 3 )
cr = addActive(cr, 1.0, 5, 3 )

cr = addPassive(cr, 1.0, 5, 4 )
cr = addActive(cr, 1.0, 5, 4 )

S = scatteringMatrix( 1.0, cr  )

const bs = beamSplitter

function evaluation(θs::Vector)
    S = solution(θs)
    T = [0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] # [0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] # bs(π/2) # I(4) # bs(π/3) * Diagonal([2, 1/2, 1/3, 3]) * bs(π/3)
    return sqrt(tr((S - T) * transpose(S - T)))
end

function solution(θs::Vector)
    return teleportation(phaseShifting([1, 1, θs[1:8]...]...)*Matrix(S)*phaseShifting([1, 1, θs[9:16]...]...), [1, 2], [1, 2])
end

res = optimize(evaluation,
                2*π*rand(16),
                method = BFGS(),
                iterations = 1000000,
                autodiff = :forward)
round.(solution(Optim.minimizer(res)), digits = 4)

# function transmission(θ)
#     res = teleportation(
#         phaseShifting(
#             [1, 1, θ, θ]...)*Matrix(S)*phaseShifting([1, 1, θ + π/3, θ + π/6]...),
#             [1, 2],
#             [1, 2])
#     return det(res[3:4, 1:2])
# end

# plot(transmission, -π, π,
#         ylims = (-2, 2),
#         xguide = L"$\theta$",
#         xguide_position = :center,
#         xticks = (π*[-1, -2/3,  -1/3,  0,  1/3,  2/3, 1],
#                     [L"\pi", L"-\frac{2\pi}{3}",  L"-\frac{\pi}{3}",  L"0",  L"\frac{\pi}{3}", L"\frac{2\pi}{3}", L"\pi"]),
#         xtickfontsize = 11,
#         yticks = ([-2, -1, 0, 1, 2],
#                  [L"-2", L"-1", L"0", L"1", L"2"]),
#         ytickfontsize = 11,
#         label = "Transmission",
#         bottom_margin = 10px )

# savefig("~/test.pdf")


text:
@testset "FOS Constructor" begin
    fo1 = PFO(FO(:c,2,Up,false),symbols(:t1))
    fo2 = PFO(FO(:c,-1,Dn,true),symbols(:t1))
    fos = FOS([fo1;fo2])
    @test isa(fos,FOS)
end