text:
using InteractiveUtils: subtypes
using Test

using DataStructures

# generic
function test_interface(::Type{T}) where T
    for S in subtypes(T)
        _test_interface(S)
    end
end
function test_interface(::Type{T}, ::Type{U}) where {T,U}
    for S in subtypes(T)
        _test_interface(S{U})
    end
end

# run tests

include("./list.jl")
test_interface(AbstractList, Int)
test_interface(AbstractQueue, Int)
test_interface(AbstractStack, Int)

include("./graph.jl")
test_interface(AbstractGraph)

include("./graph_algorithms.jl")

include("./tree.jl")

include("./sort.jl")


text:
module GroupsTests

using BenchmarkExt
using BenchmarkExt: TrialEstimate, Parameters
using ReTest

seteq(a, b) = length(a) == length(b) == length(intersect(a, b))

function setup_vals()
    g1 = BenchmarkGroup(["1", "2"])

    t1a = TrialEstimate(Parameters(time_tolerance = .05, memory_tolerance = .05), 32, 1, 2, 3)
    t1b = TrialEstimate(Parameters(time_tolerance = .40, memory_tolerance = .40), 4123, 123, 43, 9)
    tc = TrialEstimate(Parameters(time_tolerance = 1.0, memory_tolerance = 1.0), 1, 1, 1, 1)

    g1["a"] = t1a
    g1["b"] = t1b
    g1["c"] = tc

    g1copy = copy(g1)
    g1similar = similar(g1)

    g2 = BenchmarkGroup(["2", "3"])

    t2a = TrialEstimate(Parameters(time_tolerance = .05, memory_tolerance = .05), 323, 1, 2, 3)
    t2b = TrialEstimate(Parameters(time_tolerance = .40, memory_tolerance = .40), 1002, 123, 43, 9)

    g2["a"] = t2a
    g2["b"] = t2b
    g2["c"] = tc

    trial = BenchmarkExt.Trial(Parameters(), [1, 2, 5], [0, 1, 1], [3, 4, 5], [56, 58, 62])

    gtrial = BenchmarkGroup([], Dict("t" => trial))

    (g1, g1copy, g1similar, g2, gtrial, t1a, t1b, tc, t2a, t2b)
end

function setup_extra(g1, g2, gtrial)
    groupsa = BenchmarkGroup()
    groupsa["g1"] = g1
    groupsa["g2"] = g2
    g3a = addgroup!(groupsa, "g3", ["3", "4"])
    g3a["c"] = TrialEstimate(Parameters(time_tolerance = .05, memory_tolerance = .05), 6341, 23, 41, 536)
    g3a["d"] = TrialEstimate(Parameters(time_tolerance = .13, memory_tolerance = .13), 12341, 3013, 2, 150)

    groups_copy = copy(groupsa)
    groups_similar = similar(groupsa)

    groupsb = BenchmarkGroup()
    groupsb["g1"] = g1
    groupsb["g2"] = g2
    g3b = addgroup!(groupsb, "g3", ["3", "4"])
    g3b["c"] = TrialEstimate(Parameters(time_tolerance = .05, memory_tolerance = .05), 1003, 23, 41, 536)
    g3b["d"] = TrialEstimate(Parameters(time_tolerance = .23, memory_tolerance = .23), 25341, 3013, 2, 150)

    groupstrial = BenchmarkGroup()
    groupstrial["g"] = gtrial

    (groupsa, g3a, groupsb, g3b, groupstrial, groups_copy, groups_similar)
end

struct Bar end
@testset "BenchmarkGroup" begin
    g1, g1copy, g1similar, g2, gtrial, t1a, t1b, tc, t2a, t2b = setup_vals()

    @test BenchmarkGroup() == BenchmarkGroup([], Dict())
    @test length(g1) == 3
    @test g1["a"] == t1a
    @test g1["b"] == t1b
    @test g1["c"] == tc
    @test haskey(g1, "a")
    @test !(haskey(g1, "x"))
    @test seteq(keys(g1), ["a", "b", "c"])
    @test seteq(values(g1), [t1a, t1b, tc])
    @test iterate(g1) == iterate(g1.data)
    @test iterate(g1, 1) == iterate(g1.data, 1)
    @test seteq([x for x in g1], Pair["a"=>t1a, "b"=>t1b, "c"=>tc])

    @test g1 == g1copy
    @test seteq(keys(delete!(g1copy, "a")), ["b", "c"])
    @test isempty(delete!(delete!(g1copy, "b"), "c"))
    @test isempty(g1similar)
    @test g1similar.tags == g1.tags

    @test time(g1).data == Dict("a" => time(t1a), "b" => time(t1b), "c" => time(tc))
    @test gctime(g1).data == Dict("a" => gctime(t1a), "b" => gctime(t1b), "c" => gctime(tc))
    @test memory(g1).data == Dict("a" => memory(t1a), "b" => memory(t1b), "c" => memory(tc))
    @test allocs(g1).data == Dict("a" => allocs(t1a), "b" => allocs(t1b), "c" => allocs(tc))
    @test params(g1).data == Dict("a" => params(t1a), "b" => params(t1b), "c" => params(tc))

    @test max(g1, g2).data == Dict("a" => t2a, "b" => t1b, "c" => tc)
    @test min(g1, g2).data == Dict("a" => t1a, "b" => t2b, "c" => tc)
    @test ratio(g1, g2).data == Dict("a" => ratio(t1a, t2a), "b" => ratio(t1b, t2b), "c" => ratio(tc, tc))
    @test (judge(g1, g2; time_tolerance = 0.1, memory_tolerance = 0.1).data ==
           Dict("a" => judge(t1a, t2a; time_tolerance = 0.1, memory_tolerance = 0.1),
                "b" => judge(t1b, t2b; time_tolerance = 0.1, memory_tolerance = 0.1),
                "c" => judge(tc, tc; time_tolerance = 0.1, memory_tolerance = 0.1)))
    @test (judge(ratio(g1, g2); time_tolerance = 0.1, memory_tolerance = 0.1) ==
           judge(g1, g2; time_tolerance = 0.1, memory_tolerance = 0.1))
    @test ratio(g1, g2) == ratio(judge(g1, g2))

    @test isinvariant(judge(g1, g1))
    @test isinvariant(time, judge(g1, g1))
    @test isinvariant(memory, judge(g1, g1))
    @test !(isregression(judge(g1, g1)))
    @test !(isregression(time, judge(g1, g1)))
    @test !(isregression(memory, judge(g1, g1)))
    @test !(isimprovement(judge(g1, g1)))
    @test !(isimprovement(time, judge(g1, g1)))
    @test !(isimprovement(memory, judge(g1, g1)))

    @test BenchmarkExt.invariants(judge(g1, g2)).data == Dict("c" => judge(tc, tc))
    @test BenchmarkExt.invariants(time, (judge(g1, g2))).data == Dict("c" => judge(tc, tc))
    @test BenchmarkExt.invariants(memory, (judge(g1, g2))).data == Dict("a" => judge(t1a, t2a), "b" => judge(t1b, t2b), "c" => judge(tc, tc))
    @test BenchmarkExt.regressions(judge(g1, g2)).data == Dict("b" => judge(t1b, t2b))
    @test BenchmarkExt.regressions(time, (judge(g1, g2))).data == Dict("b" => judge(t1b, t2b))
    @test BenchmarkExt.regressions(memory, (judge(g1, g2))).data == Dict()
    @test BenchmarkExt.improvements(judge(g1, g2)).data == Dict("a" => judge(t1a, t2a))
    @test BenchmarkExt.improvements(time, (judge(g1, g2))).data == Dict("a" => judge(t1a, t2a))
    @test BenchmarkExt.improvements(memory, (judge(g1, g2))).data == Dict()

    @test isinvariant(judge(g1, g1))
    @test !(isinvariant(judge(g1, g2)))
    @test isregression(judge(g1, g2))
    @test !(isregression(judge(g1, g1)))
    @test isimprovement(judge(g1, g2))
    @test !(isimprovement(judge(g1, g1)))
    @test invariants(judge(g1, g2)).data == Dict("c" => judge(tc, tc))
    @test regressions(judge(g1, g2)).data == Dict("b" => judge(t1b, t2b))
    @test improvements(judge(g1, g2)).data == Dict("a" => judge(t1a, t2a))

    @testset "struct Bar" begin
        @test BenchmarkExt.invariants(Bar()) == Bar()
        @test BenchmarkExt.invariants(time, (Bar())) == Bar()
        @test BenchmarkExt.invariants(memory, (Bar())) == Bar()
        @test BenchmarkExt.regressions(Bar()) == Bar()
        @test BenchmarkExt.regressions(time, (Bar())) == Bar()
        @test BenchmarkExt.regressions(memory, (Bar())) == Bar()
        @test BenchmarkExt.improvements(Bar()) == Bar()
        @test BenchmarkExt.improvements(time, (Bar())) == Bar()
        @test BenchmarkExt.improvements(memory, (Bar())) == Bar()
    end

    @test minimum(gtrial)["t"] == minimum(gtrial["t"])
    @test median(gtrial)["t"] == median(gtrial["t"])
    @test mean(gtrial)["t"] == mean(gtrial["t"])
    @test maximum(gtrial)["t"] == maximum(gtrial["t"])
    @test params(gtrial)["t"] == params(gtrial["t"])
end

@testset "BenchmarkGroups of BenchmarkGroups" begin
    g1, g1copy, g1similar, g2, gtrial, t1a, t1b, tc, t2a, t2b = setup_vals()
    groupsa, g3a, groupsb, g3b, groupstrial, groups_copy, groups_similar = setup_extra(g1, g2, gtrial)

    @test time(groupsa).data == Dict("g1" => time(g1), "g2" => time(g2), "g3" => time(g3a))
    @test gctime(groupsa).data == Dict("g1" => gctime(g1), "g2" => gctime(g2), "g3" => gctime(g3a))
    @test memory(groupsa).data == Dict("g1" => memory(g1), "g2" => memory(g2), "g3" => memory(g3a))
    @test allocs(groupsa).data == Dict("g1" => allocs(g1), "g2" => allocs(g2), "g3" => allocs(g3a))
    @test params(groupsa).data == Dict("g1" => params(g1), "g2" => params(g2), "g3" => params(g3a))

    for (k, v) in BenchmarkExt.leaves(groupsa)
        @test groupsa[k] == v
    end

    @test max(groupsa, groupsb).data == Dict("g1" => max(g1, g1), "g2" => max(g2, g2), "g3" => max(g3a, g3b))
    @test min(groupsa, groupsb).data == Dict("g1" => min(g1, g1), "g2" => min(g2, g2), "g3" => min(g3a, g3b))
    @test ratio(groupsa, groupsb).data == Dict("g1" => ratio(g1, g1), "g2" => ratio(g2, g2), "g3" => ratio(g3a, g3b))
    @test (judge(groupsa, groupsb; time_tolerance = 0.1, memory_tolerance = 0.1).data ==
           Dict("g1" => judge(g1, g1; time_tolerance = 0.1, memory_tolerance = 0.1),
                "g2" => judge(g2, g2; time_tolerance = 0.1, memory_tolerance = 0.1),
                "g3" => judge(g3a, g3b; time_tolerance = 0.1, memory_tolerance = 0.1)))
    @test (judge(ratio(groupsa, groupsb); time_tolerance = 0.1, memory_tolerance = 0.1) ==
           judge(groupsa, groupsb; time_tolerance = 0.1, memory_tolerance = 0.1))
    @test ratio(groupsa, groupsb) == ratio(judge(groupsa, groupsb))

    @test isinvariant(judge(groupsa, groupsa))
    @test !(isinvariant(judge(groupsa, groupsb)))
    @test isregression(judge(groupsa, groupsb))
    @test !(isregression(judge(groupsa, groupsa)))
    @test isimprovement(judge(groupsa, groupsb))
    @test !(isimprovement(judge(groupsa, groupsa)))
    @test invariants(judge(groupsa, groupsb)).data == Dict("g1" => judge(g1, g1), "g2" => judge(g2, g2))
    @test regressions(judge(groupsa, groupsb)).data == Dict("g3" => regressions(judge(g3a, g3b)))
    @test improvements(judge(groupsa, groupsb)).data == Dict("g3" => improvements(judge(g3a, g3b)))

    @test minimum(groupstrial)["g"]["t"] == minimum(groupstrial["g"]["t"])
    @test maximum(groupstrial)["g"]["t"] == maximum(groupstrial["g"]["t"])
    @test median(groupstrial)["g"]["t"] == median(groupstrial["g"]["t"])
    @test mean(groupstrial)["g"]["t"] == mean(groupstrial["g"]["t"])
    @test params(groupstrial)["g"]["t"] == params(groupstrial["g"]["t"])

end

@testset "Tagging" begin
    g1, g1copy, g1similar, g2, gtrial, t1a, t1b, tc, t2a, t2b = setup_vals()
    groupsa, g3a, groupsb, g3b, groupstrial, groups_copy, groups_similar = setup_extra(g1, g2, gtrial)

    @test groupsa[@tagged "1"] == BenchmarkGroup([], "g1" => g1)
    @test groupsa[@tagged "2"] == BenchmarkGroup([], "g1" => g1, "g2" => g2)
    @test groupsa[@tagged "3"] == BenchmarkGroup([], "g2" => g2, "g3" => g3a)
    @test groupsa[@tagged "4"] == BenchmarkGroup([], "g3" => g3a)
    @test groupsa[@tagged "3" && "4"] == groupsa[@tagged "4"]
    @test groupsa[@tagged ALL && !("2")] == groupsa[@tagged !("2")]
    @test groupsa[@tagged "1" || "4"] == BenchmarkGroup([], "g1" => g1, "g3" => g3a)
    @test groupsa[@tagged ("1" || "4") && !("2")] == groupsa[@tagged "4"]
    @test groupsa[@tagged !("1" || "4") && "2"] == BenchmarkGroup([], "g2" => g2)
    @test groupsa[@tagged ALL] == groupsa
    @test groupsa[@tagged !("1" || "3") && !("4")] == similar(groupsa)

    gnest = BenchmarkGroup(["1"],
                           "2" => BenchmarkGroup(["3"], 1 => 1),
                           4 => BenchmarkGroup(["3"], 5 => 6),
                           7 => 8,
                           "a" => BenchmarkGroup(["3"], "a" => :a, (11, "b") => :b),
                           9 => BenchmarkGroup(["2"],
                                               10 => BenchmarkGroup(["3"]),
                                               11 => BenchmarkGroup()))

    @test sort(leaves(gnest), by=string) ==
    Any[(Any["2",1],1), (Any["a","a"],:a), (Any["a",(11,"b")],:b), (Any[4,5],6), (Any[7],8)]

    @test gnest[@tagged 11 || 10] == BenchmarkGroup(["1"],
                                                    "a" => BenchmarkGroup(["3"],
                                                                          (11, "b") => :b),
                                                    9 => gnest[9])

    @test gnest[@tagged "3"] == BenchmarkGroup(["1"], "2" => gnest["2"], 4 => gnest[4], "a" => gnest["a"],
                                               9 => BenchmarkGroup(["2"], 10 => BenchmarkGroup(["3"])))

    @test gnest[@tagged "1" && "2" && "3"] == BenchmarkGroup(["1"], "2" => gnest["2"],
                                                             9 => BenchmarkGroup(["2"], 10 => BenchmarkGroup(["3"])))

    k = 3 + im
    gnest = BenchmarkGroup(["1"], :hi => BenchmarkGroup([], 1 => 1, k => BenchmarkGroup(["3"], 1 => 1)), 2 => 1)

    @test gnest[@tagged "1"] == gnest
    @test gnest[@tagged "1" && !(:hi)] == BenchmarkGroup(["1"], 2 => 1)
    @test gnest[@tagged :hi && !("3")] == BenchmarkGroup(["1"], :hi => BenchmarkGroup([], 1 => 1))
    @test gnest[@tagged k] == BenchmarkGroup(["1"], :hi => BenchmarkGroup([], k => BenchmarkGroup(["3"], 1 => 1)))
end

@testset "Indexing by BenchmarkGroup" begin
    g = BenchmarkGroup()
    d = Dict("1" => 1, "2" => 2, "3" => 3)
    g["a"] = BenchmarkGroup([], copy(d))
    g["b"] = BenchmarkGroup([], copy(d))
    g["c"] = BenchmarkGroup([], copy(d))
    g["d"] = BenchmarkGroup([], copy(d))
    g["e"] = BenchmarkGroup([], "1" => BenchmarkGroup([], copy(d)),
                            "2" => BenchmarkGroup([], copy(d)),
                            "3" => BenchmarkGroup([], copy(d)))

    x = BenchmarkGroup()
    x["a"] = BenchmarkGroup([], "1" => '1', "3" => '3')
    x["c"] = BenchmarkGroup([], "2" => '2')
    x["d"] = BenchmarkGroup([], "1" => '1', "2" => '2', "3" => '3')
    x["e"] = BenchmarkGroup([], "1" => x["a"], "3" => x["c"])

    gx = BenchmarkGroup()
    gx["a"] = BenchmarkGroup([], "1" => 1, "3" => 3)
    gx["c"] = BenchmarkGroup([], "2" => 2)
    gx["d"] = BenchmarkGroup([], "1" => 1, "2" => 2, "3" => 3)
    gx["e"] = BenchmarkGroup([], "1" => g["e"]["1"][x["a"]], "3" => g["e"]["3"][x["c"]])

    @test g[x] == gx
end

@testset "Indexing by vector" begin
    g1 = BenchmarkGroup(1 => BenchmarkGroup("a" => BenchmarkGroup()))
    g1[[1, "a", :b]] = "hello"
    @test g1[[1, "a", :b]] == "hello"

    g2 = BenchmarkGroup()
    g2[[1, "a", :b]] = "hello"  # should create higher levels on the fly
    @test g2[[1, "a", :b]] == "hello"

    @test g1 == g2

    @testset "benchmarkset" begin
        g1 = @benchmarkset "test set" begin
            @case "test case 1" 1 + 1
            @case "test case 2" 2 + 2
        end

        @test haskey(g1, "test set")
        @test haskey(g1["test set"], "test case 1")
        @test haskey(g1["test set"], "test case 2")
    end
end

@testset "Pretty printing" begin
    g1, g1copy, g1similar, g2, gtrial, t1a, t1b, tc, t2a, t2b = setup_vals()

    g1 = BenchmarkGroup(["1", "2"])
    g1["a"] = t1a
    g1["b"] = t1b
    g1["c"] = tc

    data = read(joinpath(@__DIR__, "data", "test03_pretty.txt"), String)
    pp = strip.(split(data, "\n\n"))

    @test sprint(show, g1) == pp[1]
    @test sprint(show, g1; context = :boundto => 1) == pp[2]
    @test sprint(show, g1; context = :limit => false) == pp[3]
    @test @test_deprecated(sprint(show, g1; context = :limit => 1)) == pp[4]
end

end # module

ext: jl


text:
function getsymbols(C::ClimGrid)
    latsymbol = Symbol(C.dimension_dict["lat"])
    lonsymbol = Symbol(C.dimension_dict["lon"])

    return latsymbol, lonsymbol
end

"""
    merge(A::ClimGrid, B::ClimGrid)

Combines two ClimGrid. Based on the AxisArrays method.
"""
function Base.merge(A::ClimGrid, B::ClimGrid)
    axisArray = merge(A.data, B.data)
    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=A.experiment, run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:+(A::ClimGrid, B::ClimGrid)
    axisArraytmp = A.data .+ B.data

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " + ", B.experiment), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:+(A::ClimGrid, k)
    axisArraytmp = A.data .+ k

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " + ", k), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:-(A::ClimGrid, B::ClimGrid)

    axisArraytmp = A.data .- B.data

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " - ", B.experiment), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:-(A::ClimGrid, k)
    axisArraytmp = A.data .- k

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " - ", k), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:*(A::ClimGrid, B::ClimGrid)

    axisArraytmp = A.data .* B.data

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " * ", B.experiment), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:*(A::ClimGrid, k)
    axisArraytmp = A.data .* k

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " * ", k), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:/(A::ClimGrid, B::ClimGrid)
    axisArraytmp = A.data ./ B.data

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " / ", B.experiment), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

function Base.:/(A::ClimGrid, k)
    axisArraytmp = A.data ./ k

    axisArray = buildarrayinterface(axisArraytmp, A)

    ClimGrid(axisArray, longrid=A.longrid, latgrid=A.latgrid, msk=A.msk, grid_mapping=A.grid_mapping, dimension_dict=A.dimension_dict, timeattrib=A.timeattrib, model=A.model, frequency=A.frequency, experiment=string(A.experiment, " / ", k), run=A.run, project=A.project, institute=A.institute, filename=A.filename, dataunits=A.dataunits, latunits=A.latunits, lonunits=A.lonunits, variable=A.variable, typeofvar=A.typeofvar, typeofcal=A.typeofcal, varattribs=A.varattribs, globalattribs=A.globalattribs)
end

"""
    mean(A::ClimGrid)

Compute the mean of `ClimGrid` A
"""
mean(A::ClimGrid) = Statistics.mean(A[1])

"""
    minimum(A::ClimGrid)

Compute the minimum value of `ClimGrid` A
"""
Base.minimum(A::ClimGrid) = minimum(A[1])

"""
    maximum(A::ClimGrid)

Compute the maximum value of `ClimGrid` A
"""
Base.maximum(A::ClimGrid) = maximum(A[1])

"""
    std(A::ClimGrid)

Compute the standard deviation of `ClimGrid` A
"""
std(A::ClimGrid) = Statistics.std(A[1])

"""
    var(A::ClimGrid)

Compute the variance of `ClimGrid` A
"""
var(A::ClimGrid) = Statistics.var(A[1])

function getindex(C::ClimGrid,i::Int)
    if i == 1
        return C.data
    elseif i == 2
        return C.dataunits
    elseif i == 3
        return C.model
    elseif i == 4
        return C.experiment
    elseif i == 5
        return C.run
    elseif i == 6
        return C.lonunits
    elseif i == 7
        return C.latunits
    elseif i == 8
        return C.filename
    elseif i == 9
        return C.variable
    elseif i == 10
        return C.typeofvar
    elseif i == 11
        return C.typeofcal
    elseif i == 12
        return C.globalattribs
    else
        throw(error("You can't index like that!"))
    end
end

# Base.IndexStyle{T<:ClimGrid}(::Type{T}) = Base.IndexLinear()
Base.length(C::ClimGrid) = length(fieldnames(typeof(C)))
Base.size(C::ClimGrid) = (length(C),)
Base.size(C::ClimGrid, n::Int) = n==1 ? length(C) : error("Only dimension 1 has a well-defined size.")
#Base.endof(C::ClimGrid) = length(C)
Base.ndims(::ClimGrid) = 1

Base.show(io::IO, ::MIME"text/plain", C::ClimGrid) = print(io, "ClimGrid struct with data:\n   ", summary(C[1]), "\n",
"Project: ", C.project, "\n",
"Institute: ", C.institute, "\n",
"Model: ", C[3], "\n",
"Experiment: ", C[4], "\n",
"Run: ", C[5], "\n",
"Variable: ", C[9], "\n",
# "Variable CF standard name: ", C.varattribs["standard_name"], "\n",
"Data units: ", C[2], "\n",
"Frequency: ", C.frequency, "\n",
"Global attributes: ", summary(C[12]), "\n",
"Filename: ", C[8])

# Base.show(io::IO, ::MIME"text/plain", ITP::TransferFunction) = print(io, "TransferFunction type with fields *itp*, *method* and *detrend*", "\n",
#     "Interpolation array: ", size(ITP.itp), " transfer functions", "\n",
#     "Method: ", ITP.method, "\n",
#     "Detrended: ", ITP.detrend)


text:
# special Requires.jl syntax
import .Polyester: splitloop, combine, NoLoop

splitloop(t::LazyTree) = NoLoop(), eachindex(t), t
combine(t::LazyTree, ::NoLoop, j) = LazyEvent(innertable(t), j)

splitloop(e::Base.Iterators.Enumerate{LazyTree{T}}) where T = NoLoop(), eachindex(e.itr), e
combine(e::Iterators.Enumerate{LazyTree{T}}, ::NoLoop, j) where T =  @inbounds e[j]


text:
# Semantics Analysis

function semanticanalysis!(ctx)
    # 1st pass : Register global scope functions
    for decl in ctx.decls
        @alphaassert decl.type.kind == k_proc_t || decl.type.kind == k_fn_t decl.location "Unsupported global declaration type $(decl.type)"

        ctx_newsymglobal!(ctx, decl.type, decl.id, decl.location)
    end

    # 2nd pass : Verify types and push symbols to scopes
    for decl in ctx.decls
        ctx.current_func = decl
        decl_resolve!(ctx, decl, false)
    end
end

# Semantic analysis for a declaration
function decl_resolve!(ctx, decl, pushsym=true)
    if decl.type.kind == k_proc_t || decl.type.kind == k_fn_t
        @alphaassert length(ctx.scopes) == 1 decl.location "Can't declare functions in a local scope"

        ctx_pushscope!(ctx, decl)

        args_resolve!(ctx, decl.type.args, decl.location)
        stmt_resolve!(ctx, decl.body)

        if decl.id == "main"
            @alphaassert decl.type.kind == k_fn_t decl.location "main must be a function returning an int (the exit status), not a procedure"
        end

        if decl.type.kind == k_fn_t
            @alphaassert(isreturning(ctx, decl.body),
                        decl.location,
                        "The function $(decl.id) must return a value, use proc instead of fun" *
                        " if this function must return nothing")
        end

        ctx_popscope!(ctx)
    elseif decl.type.kind == k_int_t
        exp_resolve!(ctx, decl.value)

        @alphaassert decl.type == decl.value.type decl.location "Invalid type of expression, must be int"
    else
        alphaerror("Unsupported type '$(decl.type)' for declaration", ctx, decl.location)
    end

    # Declare this symbol to the inner scope
    if pushsym
        decl.sym = ctx_newsymlocal!(ctx, decl, decl.location)
    end
end

function args_resolve!(ctx, args, location)
    for (i, arg) in enumerate(args)
        ctx_newsymarg!(ctx, t_int, arg.id, i, location)
    end
end

function stmt_resolve!(ctx, stmt)
    if stmt == nothing
        return
    end

    if stmt.kind == k_stmt_decl
        decl_resolve!(ctx, stmt.decl)
    elseif stmt.kind == k_stmt_exp
        exp_resolve!(ctx, stmt.exp)
    elseif stmt.kind == k_stmt_return
        exp_resolve!(ctx, stmt.exp)

        # Returns have only ints (empty return not allowed)
        @alphaassert stmt.exp.type.kind == k_int_t stmt.exp.location "Functions must return ints"

        # Check this is not a procedure
        @alphaassert ctx.current_func.type.kind == k_fn_t stmt.location "Cannot return within procedures"
    elseif stmt.kind == k_stmt_ifelse
        exp_resolve!(ctx, stmt.exp)
        stmt_resolve!(ctx, stmt.ifbody)
        stmt_resolve!(ctx, stmt.elsebody)

        @alphaassert stmt.exp.type.kind == k_int_t stmt.location "Invalid condition expression"
    elseif stmt.kind == k_stmt_loop
        ctx_pushscope!(ctx)

        stmt_resolve!(ctx, stmt.initbody)
        exp_resolve!(ctx, stmt.exp)
        stmt_resolve!(ctx, stmt.iterbody)
        stmt_resolve!(ctx, stmt.loopbody)

        ctx_popscope!(ctx)

        @alphaassert stmt.exp.type.kind == k_int_t stmt.location "Invalid condition expression"
    elseif stmt.kind == k_stmt_block
        ctx_pushscope!(ctx)

        for child in stmt.stmts
            stmt_resolve!(ctx, child)
        end

        ctx_popscope!(ctx)
    elseif stmt.kind == k_stmt_chain
        for child in stmt.stmts
            stmt_resolve!(ctx, child)
        end

        # Generate this statement as a block
        stmt.kind = k_stmt_block
    elseif stmt.kind == k_stmt_printint
        exp_resolve!(ctx, stmt.exp)

        @alphaassert stmt.exp.type.kind == k_int_t stmt.exp.location "Cannot display this value"
    elseif stmt.kind == k_stmt_scan
        exp_resolve!(ctx, stmt.exp)

        @alphaassert stmt.exp.type.kind == k_int_t stmt.exp.location "Invalid type for scan"
    elseif stmt.kind ∈ (k_stmt_printstr,)
        # Ignored
    else
        alphaerror("Unsupported statement kind $(stmt.kind)", ctx, stmt.location)
    end
end

function exp_resolve!(ctx, exp)
    if exp == nothing
        return
    end

    # Resolve left and right
    if exp.left != nothing
        exp_resolve!(ctx, exp.left)
        exp_resolve!(ctx, exp.right)

        if exp.right != nothing && exp.right.type != exp.left.type
            alphaerror("Invalid type for this expression", ctx, stmt.location)
        end

        exp.type = exp.left.type
    end

    if exp.kind == k_exp_id
        exp.sym = ctx_fetchscope(ctx, exp.id)

        @alphaassert exp.sym != nothing exp.location "Not found variable named '$(exp.id)'"

        exp.type = exp.sym.type
    elseif exp.kind == k_exp_int
        exp.type = t_int
    elseif exp.kind == k_exp_set
        @alphaassert exp.left.sym != nothing exp.location "Can't assign an rvalue"
    elseif exp.kind == k_exp_call
        # Check args
        for (i, arg) in enumerate(exp.args)
            exp_resolve!(ctx, arg)

            @alphaassert arg.type.kind == k_int_t exp.location "Invalid arg #$(i + 1) type, must be int"
        end

        fun = ctx_fetchglobal(ctx, exp.id)
        exp.type = fun.type.kind == k_fn_t ? t_int : t_void

        @alphaassert fun != nothing exp.location "Function $(exp.id) not declared"
        @alphaassert length(fun.type.args) == length(exp.args) exp.location "Invalid number of arguments to call $(fun.id), $(length(fun.type.args)) args required"
    elseif exp.kind == k_exp_test
        @alphaassert haskey(test_operators, exp.operator) exp.location "Invalid conditional operator $(exp.operator)"
    end
end

# Returns whether the function returns in all cases an int
function isreturning(ctx, stmt)
    if stmt == nothing
        return false
    elseif stmt.kind == k_stmt_return
        return true
    elseif stmt.kind ∈ (k_stmt_block, k_stmt_chain)
        for s in stmt.stmts
            isreturning(ctx, s) && (return true)
        end
    elseif stmt.kind == k_stmt_ifelse
        return isreturning(ctx, stmt.ifbody) && isreturning(ctx, stmt.elsebody)
    end

    return false
end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
#using Revise
using CutCell
include("useful_routines.jl")

polyorder = 1
numqp = 2

basis = TensorProductBasis(2, polyorder)
levelset = InterpolatingPolynomial(1, basis)
nf = CutCell.number_of_basis_functions(basis)
mesh = CutCell.Mesh([0.0, 0.0], [3.0, 1.0], [3, 1], nf)
nodalcoordinates = CutCell.nodal_coordinates(mesh)

normal = [1.0, 0.0]
x0 = [0.5, 0.0]
levelsetcoeffs = plane_distance_function(nodalcoordinates, normal, x0)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

lambda, mu = (1.0, 2.0)
stiffness = CutCell.HookeStiffness(lambda,mu,lambda,mu)

bilinearforms =
    CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)

@test length(bilinearforms.cellmatrices) == 4
testcelltomatrix = [
    3 1 1
    4 0 0
]
@test allequal(testcelltomatrix, bilinearforms.celltomatrix)


text:
import Base.getindex
import Base.setindex!
import Base.+
import Base.eltype
struct DiskFourTensor{T}
    """
    Data structure for Rank 4 tensors
    """
    fname::String
    dname::String
    sz1::Int
    sz2::Int
    sz3::Int
    sz4::Int
    dtype::Type
end
function squeeze(A::AbstractArray)
    #singleton_dims = tuple((d for d in 1:ndims(A) if size(A, d) == 1)...)
    return dropdims(A, dims = (findall(size(A) .== 1)...,))
end
function DiskFourTensor{T}(
    fname::String,
    sz1::Int,
    sz2::Int,
    sz3::Int,
    sz4::Int,
    mode::String = "r+",
) where T
    """
    Constructor for DiskFourTensor objects
    """
    file = h5open(fname, mode)
    if (mode == "w") | (mode == "w+")
        dataset = d_create(file, "data", datatype(dtype), dataspace(sz1))
    else
        dataset = file["data"]
    end
    close(file)
    DiskFourTensor(fname, "data", sz1, sz2, sz3, sz4, dtype)
end
# >>> overload getindex ( A[i,j,k,l] ) syntax

function getindex(
    dtens::DiskFourTensor,
    i1::Union{UnitRange{Int64},Int64,Colon},
    i2::Union{UnitRange{Int64},Int64,Colon},
    i3::Union{UnitRange{Int64},Int64,Colon},
    i4::Union{UnitRange{Int64},Int64,Colon},
)
    h5open(dtens.fname, "r") do fid
        squeeze(fid["$dtens.dname"][ranger(i1), ranger(i2), ranger(i3), ranger(i4)])
    end
end
# <<<

# >>> overload setindex! ( A[i,j,k,l] = 2.0 )
function setindex!(
    dtens::DiskFourTensor,
    val,
    i1::UnitRange{Int64},
    i2::UnitRange{Int64},
    i3::UnitRange{Int64},
    i4::UnitRange{Int64},
)
    h5open(dtens.fname, "r+") do fid
        fid["$dtens.dname"][i1, i2, i3, i4] = val
    end
end
function setindex!(
    dtens::DiskFourTensor,
    val,
    i1::Union{UnitRange{Int64},Int64,Colon},
    i2::Union{UnitRange{Int64},Int64,Colon},
    i3::Union{UnitRange{Int64},Int64,Colon},
    i4::Union{UnitRange{Int64},Int64,Colon},
)
    h5open(dtens.fname, "r+") do fid
        fid["$dtens.dname"][ranger(i1), ranger(i2), ranger(i3), ranger(i4)] = val
    end
end
# <<< overload setindex!
# >>> overload eltype
function eltype(dtens::DiskFourTensor)
    return dtens.dtype
end
# <<< overload eltype

function blockfill!(dtens::DiskFourTensor, val)
    """
    Fill a DiskFourTensor with a single value.
    """
    A = zeros(Float64, dtens.sz1, dtens.sz2, dtens.sz3, dtens.sz4)
    A .= val
    h5write(dtens.fname, "$dtens.dname", A)
end


text:
module CoinbasePro

using HTTP
using JSON
using Dates
using DataFrames
using TimesDates

const ENDPOINT = "https://api.pro.coinbase.com/"

include("trades.jl")
include("utils.jl")
include("products.jl")
include("stats.jl")
include("ticker.jl")
include("book.jl")
include("candles.jl")

export trades, lasttrades, products, product, stats, ticker, book, candles

end


text:
@info "FunctionalStateMachine.jl is adding Graphs.jl related tools (Visualization)."

export
  histGraphStateMachineTransitions,
  drawStateTransitionStep,
  drawStateMachineHistory,
  animateStateMachineHistoryByTime,
  animateStateMachineHistoryByTimeCompound,
  animateStateMachineHistoryIntervalCompound

# bad import warning
# import Graphs: incdict
# incdict(::Type{V}, ::Type{E}; is_directed::Bool = true) where {V,E} = incdict(Dict{Int,V}(), E{V}; is_directed=is_directed)

"""
    $SIGNATURES

Create a `Graphs.incdict` object and populate with nodes (states) and edges (transitions)
according to the contents of parameters passed in.

Notes:
- Current implementation repeats duplicate transitions as new edges.
"""
function histGraphStateMachineTransitions(stateVisits, allStates::Vector{Symbol};
                                          maxpenwidth::Real=5, minpenwidth::Real=0.75 )
  #
  g = Graphs.incdict(Dict{Int,Graphs.ExVertex}(), Graphs.ExEdge{Graphs.ExVertex}; is_directed=true)
  # g = Graphs.incdict(Graphs.ExVertex,Graphs.ExEdge,is_directed=true)
  lookup = Dict{Symbol, Int}()

  # add all required states as nodes to the visualization graph
  fid = 0
  for state in allStates
    fid += 1
    lbl = string(state)
    exvert = Graphs.ExVertex(fid, lbl)
    exvert.attributes["label"] = lbl
    Graphs.add_vertex!(g, exvert)
    lookup[state] = fid
  end

  alledges = []

  # add all edges to graph
  maxtransedge = 0
  count = 0
  for (from, tos) in stateVisits
    for to in tos
      count += 1
      exvf = g.vertices[lookup[from]]
      exvt = g.vertices[lookup[to]]
      # add the edge from one to the next state
      # TODO, don't add if already there.
      addedge = true
      for oun in Graphs.out_neighbors(exvf, g)
        if oun.index == exvt.index
          addedge = false
          # increase penwidth+=1 on that edge
          for ed in Graphs.out_edges(exvf, g)
            ed.attributes["penwidth"] += 1
            # ed.attributes["penwidth"] = minimum([maxpenwidth;ed.attributes["penwidth"]])
            maxtransedge = maxtransedge < ed.attributes["penwidth"] ? ed.attributes["penwidth"] : maxtransedge
          end
          break
        end
      end
      if addedge
        edge = Graphs.make_edge(g, exvf, exvt)
        edge.attributes["penwidth"] = 1.0
        Graphs.add_edge!(g, edge)
        push!(alledges, edge)
      end
    end
  end

  # normalize edge penwidth to maxpenwidth
  normwidth = maxtransedge/maxpenwidth
  for ed in alledges
    ed.attributes["penwidth"] = maximum([minpenwidth; ed.attributes["penwidth"]/normwidth])

  end

  return g, lookup
end

function renderStateMachineFrame( vg,
                                  frame::Int;
                                  title::String="",
                                  viewerapp::String="eog",
                                  fext::String="png",
                                  dpi::Int=100,
                                  engine::String="dot",
                                  show::Bool=true,
                                  folder::String="fsm_animation",
                                  folderpath = "/tmp/$folder/",
                                  timest::String="",
                                  rmfirst::Bool=false  )
  #
  if rmfirst
    @warn "removing contents of $(folderpath)"
    Base.rm(folderpath, recursive=true, force=true)
  end
  mkpath(folderpath)
  dotfile = "$folderpath/csm_$frame.dot"
  filepath= "$folderpath/csm_$frame.$(fext)"

  # dot file handle
  fid = open(dotfile,"w")
  write(fid,Graphs.to_dot(vg))
  close(fid)

  # build the dot file somewhat manually
  fid = open("$folderpath/dotscript.sh","w")
  str = "head -n `wc -l $dotfile | awk '{print \$1-1}'` $dotfile > $folderpath/tmpdot.dot"
  println(fid, str)
  println(fid, "echo \"graph [label=\\\"$title, $(timest)\\\", labelloc=t];\" >> $folderpath/tmpdot.dot")
  println(fid, "echo \"}\" >> $folderpath/tmpdot.dot")
  close(fid)
  run(`chmod u+x $folderpath/dotscript.sh`)
  run(`sh $folderpath/dotscript.sh`)
  Base.rm(dotfile)
  Base.rm("$folderpath/dotscript.sh")
  run(`mv $folderpath/tmpdot.dot $dotfile`)

  # compile output and maybe show to user
  run(`$(engine) $(dotfile) -T$(fext) -Gdpi=$dpi -o $(filepath)`)
  show ? (@async run(`$viewerapp $filepath`)) : nothing
  return filepath
end

function setVisGraphOnState!( vg, vertid;
                              xlabel::String="",
                              appendxlabel::String="",
                              vertColor::AbstractString="red" )
  #
  vg.vertices[vertid].attributes["fillcolor"] = vertColor
  vg.vertices[vertid].attributes["style"] = "filled"
  if length(xlabel) > 0
    vg.vertices[vertid].attributes["xlabel"] = xlabel
  end
  if haskey(vg.vertices[vertid].attributes, "xlabel")
    vg.vertices[vertid].attributes["xlabel"] = vg.vertices[vertid].attributes["xlabel"]*appendxlabel
  elseif length(appendxlabel) > 0
    vg.vertices[vertid].attributes["xlabel"] = appendxlabel
  end
end

function clearVisGraphAttributes!(vg)
  for (vid,vert) in vg.vertices
    haskey(vert.attributes, "fillcolor") ? delete!(vert.attributes, "fillcolor") : nothing
    haskey(vert.attributes, "style") ? delete!(vert.attributes, "style") : nothing
    haskey(vert.attributes, "xlabel") ? delete!(vert.attributes, "xlabel") : nothing
  end
  nothing
end

function drawStateTransitionStep( hist,
                                  step::Int,
                                  vg,
                                  lookup::Dict{Symbol,Int};
                                  title::String="",
                                  viewerapp::String="eog",
                                  fext::String="png",
                                  engine::String="dot",
                                  show::Bool=true,
                                  folder::String="",
                                  frame::Int=step,
                                  vertColor::AbstractString="red"  )
  #

  lbl = getStateLabel(hist[step][3])
  vertid = lookup[lbl]
  vert = vg.vertices[vertid]

  fillcolorbefore = haskey(vert.attributes, "fillcolor") ? deepcopy(vg.vertices[vertid].attributes["fillcolor"]) : nothing
  stylebefore = haskey(vert.attributes, "style") ? deepcopy(vg.vertices[vertid].attributes["style"]) : nothing
  xlabelbefore = haskey(vg.vertices[vertid].attributes,"xlabel") ? deepcopy(vg.vertices[vertid].attributes["xlabel"]) : nothing
  # delete!(vert.attributes, "fillcolor")
  # delete!(vert.attributes, "style")

  # identify and set the node
  xlabel = length(title) > 0 ? (xlabelbefore !== nothing ? xlabelbefore*"," : "")*title : ""
  setVisGraphOnState!(vg, vertid, xlabel=xlabel, vertColor=vertColor )

  # render state machine frame
  filepath = renderStateMachineFrame( vg,
                                      frame,
                                      title=title,
                                      viewerapp=viewerapp,
                                      fext=fext,
                                      engine=engine,
                                      show=show,
                                      folder=folder,
                                      timest=string(split(string(hist[step][1]),'T')[end]),
                                      rmfirst=false)
  #

  # clean up the vg structure
  fillcolorbefore === nothing ? delete!(vert.attributes, "fillcolor") : (vert.attributes["fillcolor"]=fillcolorbefore)
  stylebefore === nothing ? delete!(vert.attributes, "style") : (vert.attributes["style"]=stylebefore)
  xlabelbefore === nothing ? delete!(vert.attributes, "xlabel") : (vert.attributes["xlabel"]=xlabelbefore)

  return filepath
end

function drawStateMachineHistory(hist; show::Bool=false, folder::String="" )

  stateVisits, allStates = histStateMachineTransitions(hist)

  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  for i in 1:length(hist)
    drawStateTransitionStep(hist, i, vg, lookup, folder=folder, show=show  )
  end

  return nothing
end

"""
    $SIGNATURES

Draw simultaneously separate time synchronized frames from each of the desired
state machines.  These images can be produced into synchronous side-by-side videos
which allows for easier debugging and comparison of concurrently running state
machines.
"""
function animateStateMachineHistoryByTime(hist::Vector{Tuple{DateTime, Int, <: Function, T}};
                                          frames::Int=100,
                                          folder="animatestate",
                                          title::String="",
                                          show::Bool=false,
                                          startT=hist[1][1],
                                          stopT=hist[end][1],
                                          rmfirst::Bool=true,
                                          vertColor::AbstractString="red"  ) where T
  #
  stateVisits, allStates = histStateMachineTransitions(hist)

  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  totT = stopT - startT
  totT = Millisecond(round(Int, 1.05*totT.value))
  # totT *= 1.05

  step = 1
  len = length(hist)
  @showprogress "exporting state machine images, $title " for i in 1:frames
    aniT = Millisecond(round(Int, i/frames*totT.value)) + startT
    # aniT = i/frames*totT + startT
    if hist[step][1] < aniT && step < len
      step += 1
    end
    drawStateTransitionStep(hist, step, vg, lookup, title=title, folder=folder, show=show, frame=i, vertColor=vertColor  )
  end

  nothing
end

function animateStateMachineHistoryByTimeCompound(hists::Dict{Symbol, Vector{Tuple{DateTime, Int, <: Function, T}}},
                                                  startT,
                                                  stopT;
                                                  frames::Int=100,
                                                  folder="animatestate",
                                                  title::String="",
                                                  show::Bool=false,
                                                  clearstale::Bool=true,
                                                  rmfirst::Bool=true,
                                                  fsmColors::Dict{Symbol,String}=Dict{Symbol,String}(),
                                                  defaultColor::AbstractString="gray",
                                                  autocolor_cb::Function=(histstep,csym,aniT)->(haskey(fsmColors, csym) ? fsmColors[csym] : defaultColor)  ) where T
  #
  # Dict{Symbol, Vector{Symbol}}
  stateVisits = Dict{Symbol, Vector{Symbol}}()
  allStates = Vector{Symbol}()
  for (csym,hist) in hists
    stateVisits, allStates = histStateMachineTransitions(hist,allStates=allStates, stateVisits=stateVisits  )
  end

  #
  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  # total draw time and step initialization
  totT = stopT - startT
  totT = Millisecond(round(Int, 1.05*totT.value))
  histsteps = ones(Int, length(hists))

  # clear any stale state
  clearstale ? clearVisGraphAttributes!(vg) : nothing

  # loop across time
  @showprogress "exporting state machine images, $title " for i in 1:frames
    # calc frame time
    aniT = Millisecond(round(Int, i/frames*totT.value)) + startT

    # loop over all state machines
    histidx = 0
    for (csym, hist) in hists
      histidx += 1
      step = histsteps[histidx]
      len = length(hist)
      if hist[step][1] < aniT && step < len
        histsteps[histidx] += 1
      end
      # redefine after +1
      step = histsteps[histidx]

      # modify vg for each history
      lbl = getStateLabel(hist[step][3])
      vertid = lookup[lbl]
      vertColor=autocolor_cb(hist[step], csym, aniT)
      # vertColor=haskey(fsmColors, csym) ? fsmColors[csym] : defaultColor
      setVisGraphOnState!(vg, vertid, appendxlabel=string(csym)*",", vertColor=vertColor )
    end

    # finally render one frame
    renderStateMachineFrame(vg,
                            i,
                            title=title,
                            show=false,
                            folder=folder,
                            timest=string(split(string(aniT),' ')[1]),
                            rmfirst=false  )
    #
    clearVisGraphAttributes!(vg)
  end

end

# count the total number of transitions contained in hists
function getTotalNumberSteps( hists::Dict{Symbol, Vector{Tuple{DateTime, Int, <: Function, T}}} ) where T
  totSteps = 0
  for (whId, hist) in hists, hi in hist
    totSteps += 1
  end
  return totSteps
end

# point to the start step among all history steps
function getFirstStepHist( hists::Dict{Symbol, Vector{Tuple{DateTime, Int, <: Function, T}}} ) where T
  startTime = now()
  maxTime = DateTime(0)
  # NOTE, this whichId=:null is super important to ensure rendering loop can exit properly
  whichId, whichStep = :null, 0
  for (whId, hist) in hists, (st,hi) in enumerate(hist)
    if hi[1] < startTime
      # new starting point indicator
      whichId = whId
      whichStep = st
      startTime = hi[1]
    end
    if maxTime < hi[1]
      maxTime = hi[1]
    end
  end
  return whichId, whichStep, startTime, maxTime
end

# give the next step, closest in time and that has not previously been added to `prevList`.
# Also update prevList
function getNextStepHist!(hists,
                          intuple::Tuple{Symbol, Int, DateTime},
                          maxTime::DateTime,
                          prevList::Dict{Symbol, Vector{Int}} )
  #
  oldId, oldStep, oldT = intuple

  whichId, whichStep, newT = :null, 0, maxTime
  for (whId, hist) in hists, (st,hi) in enumerate(hist)
    # make sure all options are populated in previous list tracker
    if !haskey(prevList, whId)  prevList[whId] = Int[]; end
    if oldT < hi[1] && Millisecond(0) <= (hi[1] - oldT) < (newT-oldT) &&
        !(st in prevList[whId])        # must be a different step than before
      # new closest next step
      whichId = whId
      whichStep = st
      newT = hi[1]
    end
  end

  # register this step has previously been taken
  if !haskey(prevList, whichId)
    prevList[whichId] = Int[]
  end
  push!(prevList[whichId], whichStep)

  return whichId, whichStep, newT
end

# for slower movies, use a slower fps
# run(`ffmpeg -r 10 -i /tmp/caesar/csmCompound/csm_%d.png -c:v libtheora -vf fps=5 -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" -q 10 /tmp/caesar/csmCompound/out.ogv`)
# @async run(`totem /tmp/caesar/csmCompound/out.ogv`)
# draw_more_cb(::Tuple, ::Int, ::String)
function animateStateMachineHistoryIntervalCompound(hists::Dict{Symbol, Vector{Tuple{DateTime, Int, <: Function, T}}};
                                                    easyNames::Dict{Symbol,N}=Dict{Symbol,Nothing}(),
                                                    interval::Int=2,
                                                    folderpath="/tmp/animatestate",
                                                    dpi::Int=100,
                                                    title::String="",
                                                    show::Bool=false,
                                                    clearstale::Bool=true,
                                                    rmfirst::Bool=true,
                                                    draw_more_cb::Function=(x...)->(),
                                                    fsmColors::Dict{Symbol,String}=Dict{Symbol,String}(),
                                                    defaultColor::AbstractString="red",
                                                    autocolor_cb::Function=(histstep,csym,aniT)->(haskey(fsmColors, csym) ? fsmColors[csym] : defaultColor)  ) where {T, N}
  #
  # Dict{Symbol, Vector{Symbol}}
  stateVisits = Dict{Symbol, Vector{Symbol}}()
  allStates = Vector{Symbol}()
  for (csym,hist) in hists
    stateVisits, allStates = histStateMachineTransitions(hist,allStates=allStates, stateVisits=stateVisits  )
  end

  #
  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  # total draw time and step initialization
  # totT = stopT - startT
  # totT = Millisecond(round(Int, 1.05*totT.value))
  # histsteps = ones(Int, length(hists))

  # clear any stale state
  clearstale ? clearVisGraphAttributes!(vg) : nothing

  totSteps = getTotalNumberSteps(hists)
  whId, fsmStep, aniT, maxTime = getFirstStepHist(hists)
  prevList = Dict{Symbol, Vector{Int}}()
  latestList = Dict{Symbol, Int}(whId => fsmStep)

  prevT = aniT
  frameCount = 0
  # loop across time
  @showprogress "exporting state machine images, $title " for stepCount in 1:totSteps
    # which step among the hist fsms is next
    if 1 < stepCount
      # skip first would-be repeat
      whId, fsmStep, aniT = getNextStepHist!(hists, (whId, fsmStep, aniT), maxTime, prevList)
      latestList[whId] = fsmStep
    end

    # loop over all state "known" machines
    for (csym, lstep) in latestList
      # terminate at end of drawing sequence
      csym == :null ? break : nothing
      # modify vg for each history
      lbl = getStateLabel(hists[csym][lstep][3])
      vertid = lookup[lbl]
      vertColor=autocolor_cb(hists[csym][lstep], csym, aniT)
      easyn = haskey(easyNames, csym) ? easyNames[csym] : csym
      setVisGraphOnState!(vg, vertid, appendxlabel="($easyn.$lstep),", vertColor=vertColor )
    end

    # and draw as many frames for that setup
    for itr in 1:interval
      # increment frame counter
      frameCount += 1
      deltaT = (aniT - prevT).value
      prevT = aniT
      # finally render one frame
      renderStateMachineFrame(vg,
                              frameCount,
                              title=title*" || dt=$deltaT ms ||",
                              show=false,
                              dpi=dpi,
                              folderpath=folderpath,
                              timest=string(split(string(aniT),' ')[1]),
                              rmfirst=false  )
      #
      # terminate at end of drawing sequence
      whId == :null ? break : nothing
      draw_more_cb(hists[whId][fsmStep], frameCount, folderpath)
    end
    # clear current frame in prep for the next interval
    clearVisGraphAttributes!(vg)
  end

end


text:
using Flux.Optimise:update!

function get_energy(K::Matrix{T}, samples) where T <: Real
    energy = sum(samples .* (K*samples), dims=1)
end

function free_energy(K::Matrix{T}, model::AbstractSampler, samples) where T <: Real
    return mean(get_energy(K, samples) .+ get_logp(model, samples))
end

function loss(K::Matrix{T}, model::AbstractSampler, nbatch::Int) where T <: Real
    samples = gen_samples(model, nbatch)
    free_energy(K, model, samples)
end

function loss_reinforce(K::Matrix{T}, model::AbstractSampler, samples) where T <: Real
    e = get_energy(K, samples)
    logp = get_logp(model, samples)
    f = e .+ logp
    b = mean(f)
    return mean(logp.* (f .- b))
end

function grad_model(K::Matrix{T}, model::AbstractSampler, samples) where T <: Real
    grad = gradient(loss_reinforce, K, model, samples)[2]
    unpack_gradient(model, grad)
end

function train(K::Matrix{T}, model::AbstractSampler; optimizer=ADAM(0.1), nbatch::Int=100, niter::Int=100) where T <: Real
    θ = model_parameters(model)
    for i = 1:niter
        _, gθ, _ = gradient(loss, K, model, nbatch)
        update!.(Ref(optimizer), θ, gθ)
        model_dispatch!(model, θ)
        println("$i, Free Energy = ", loss(K, model, nbatch))
    end
end


text:
"""
Shock and Detonation Toolbox
"cv" module

Calculates constant-volume explosions.

This module defines the following functions:

    cvsolve

################################################################################
Theory, numerical methods and applications are described in the following report:

    SDToolbox Numerical Tools for Shock and Detonation Wave Modeling,
    Explosion Dynamics Laboratory, Contributors:
    S. Browne, J. Ziegler, N. Bitter, B. Schmidt, J. Lawson and J. E. Shepherd,
    GALCIT Technical Report FM2018.001 Revised January 2021.
    California Institute of Technology, Pasadena, CA USA

Please cite this report and the website if you use these routines.

Please refer to LICENCE.txt or the above report for copyright and disclaimers.

http://shepherd.caltech.edu/EDL/PublicResources/sdt/

"""
module CV

export cvsolve

using OrdinaryDiffEq
using PyCall

const ct = PyNULL()

function __init__()
    copy!(ct, pyimport_conda("cantera","cantera","cantera"))
    global R̄ = ct.gas_constant
end

"""
    cv!(dy::Vector{Float64},y::Vector{Float64},params,t::Real)::Vector{Float64}

Evaluates the system of ordinary differential equations for an adiabatic,
constant-volume, zero-dimensional reactor.
It assumes that the "gas" object represents a reacting ideal gas mixture.

INPUT:
    t = time
    y = solution array [temperature, species mass 1, 2, ...]
    gas = working gas object

OUTPUT:
    An array containing time derivatives of:
        temperature and species mass fractions,
    formatted in a way that the integrator in cvsolve can recognize.

"""
function cv!(dy::Vector{Float64},y::Vector{Float64},params,t::Real)::Vector{Float64}

    # current gas state and constant parameters
    gas::PyObject,ρ₁::Real,Mᵢ::Vector{Float64} = params

    # make sure physical bounds are respected to avoid cantera errors
    if y[1] >= 0
        T, Y = y[1], y[2:end]
        gas.TDY = T, ρ₁, Y
    else # if unphysical, just continue with latest working values and hope step size is adjusted
        T, _, Y = gas.TDY::Tuple{Float64, Float64, Vector{Float64}}
        println("out of bounds for T")
    end

    ω̇  = gas.net_production_rates::Vector{Float64}
    hₛ_RT = gas.standard_enthalpies_RT
    cᵥ = gas.cv_mass::Float64
    Ṫ = -T * R̄ /(ρ₁*cᵥ) * sum((hₛ_RT .- 1) .* ω̇ )
    Ẏ = ω̇  .* Mᵢ ./ ρ₁

    dy[1] = Ṫ
    dy[2:end] = Ẏ

    return dy
end

"""
    cvsolve(gas::PyObject;t_end::Real=1e-6,max_step::Real=1e-5,
                t_eval=nothing,relTol::Real=1e-5,absTol::Real=1e-8)

Solves the ODE system defined in cv!, taking the gas object input as the
initial state.

FUNCTION SYNTAX:
    output = cvsolve(gas,**kwargs)

INPUT:
    gas = working gas object

OPTIONAL INPUT:
    t_end = end time for integration, in sec
    max_step = maximum time step for integration, in sec
    t_eval = array of time values to evaluate the solution at.
                If left as "None", solver will select values.
                Sometimes these may be too sparse for good-looking plots.
    relTol = relative tolerance
    absTol = absolute tolerances

OUTPUT:
    output = a dictionary containing the following results:
        time = time array
        T = temperature profile array
        P = pressure profile array
        speciesY = species mass fraction array
        speciesX = species mole fraction array

        gas = working gas object

        exo_time = pulse width (in secs) of temperature gradient (using 1/2 max)
        ind_time = time to maximum temperature gradient
        ind_time_10 = time to 10% of maximum temperature gradient
        ind_time_90 = time to 90% of maximum temperature gradient

"""
function cvsolve(gas::PyObject;t_end::Real=1e-6,max_step::Real=1e-5,
                t_eval=nothing,relTol::Real=1e-5,absTol::Real=1e-8)

    global ρ₁ = gas.density::Float64
    y₀ = vcat(gas.T::Float64,gas.Y::Vector{Float64})
    global Mᵢ = gas.molecular_weights::Vector{Float64}

    tel = [0.,t_end] # Timespan

    # determine adiabatic flame temperature
    old_state = gas.state
    gas.equilibrate("UV")
    T_ad = gas.T::Float64
    gas.state = old_state

    function approaches_equilibrium(y::Vector{Float64},t::Real,integrator)
        """ Returns true if T >= 0.999 * T_ad """
        if y[1] >= 0.999*T_ad && get_du(integrator)[1] <= 1e-6
            println("Reached thermal equilibrium. Terminating!")
            return true
        else
            return false
        end
    end

    params = [gas,ρ₁,Mᵢ]

    # Discrete Callback to terminate ODE solver when approaching equilibrium
    cb = DiscreteCallback(approaches_equilibrium,terminate!)
    prob = ODEProblem(cv!,y₀,tel,params)

    @time begin
        # Benchmarks needed for: abstol, reltol, Algos: Rosenbrock23, RadauIIA5, Rodas4
        out = solve(prob,Rosenbrock23(autodiff=false),progress=true,callback=cb,abstol=absTol,reltol=relTol)
    end

    output = create_output_dict(out,gas)

end

function create_output_dict(ode_output::ODESolution,gas::PyObject)

    output = Dict() # could probably be optimized, but for now for convinience

    output["time"] = ode_output.t
    output["T"] = [u[1] for u in ode_output.u]
    output["speciesY"] = [u[2:end] for u in ode_output.u]

    # Initialize additional output matrices where needed
    b = length(output["time"])
    output["P"] = zeros(0)
    output["dTdt"] = zeros(0)
    output["speciesX"] = Vector{Float64}[]
    output["ind_time"] = 0
    output["ind_time_90"] = 0
    output["ind_time_10"] = 0
    output["exo_time"] = 0
    temp_grad = zeros(0)

    #############################################################################
    # Extract PRESSSURE and TEMPERATURE GRADIENT
    #############################################################################

    # Have to loop for operations involving the working gas object

    for (T,Y) in zip(output["T"],output["speciesY"])
        gas.TDY = T,ρ₁,Y
        M̄ = gas.mean_molecular_weight
        dTdt = cv!(zeros(length(Y)+1),vcat(T,Y),[gas,ρ₁,Mᵢ],0.)[1]::Float64

        push!(output["P"],gas.P)
        push!(output["speciesX"],gas.X)
        push!(output["dTdt"],dTdt)
    end

    n = argmax(output["dTdt"])

    if n == b
        println("Error: Maximum temperature gradient occurs at the end of the reaction zone")
        println("       Your final integration length may be too short,")
        println("       your mixture may be too rich/lean, or something else may be wrong")
        println(" ")
        output["ind_time"] = output["time"][b]
        output["ind_time_10"] = output["time"][b]
        output["ind_time_90"] = output["time"][b]
        output["exo_time"] = 0
        println("Induction Time: $(output["ind_time"])")
        println("Exothermic Pulse Time: $(output["exo_time"])")
        return output
    elseif n == 0
        println("Error: Maximum temperature gradient occurs at the beginning of the reaction zone")
        println("       Your final integration length may be too short,")
        println("       your mixture may be too rich/lean, or something else may be wrong")
        println(" ")
        output["ind_time"] = output["time"][1]
        output["ind_time_10"] = output["time"][1]
        output["ind_time_90"] = output["time"][1]
        output["exo_time"] = 0
        println("Induction Time: $(output["ind_time"])")
        println("Exothermic Pulse Time: $(output["exo_time"])")
        return output
    else
        output["ind_time"] = output["time"][n]

        k = 1
        MAX10 = 0.1*maximum(output["dTdt"])
        d = output["dTdt"][1]
        while d < MAX10 && k < n
            k = k + 1
            d = output["dTdt"][k]
        end
        output["ind_time_10"] = output["time"][k]

        k = 1
        MAX90 = 0.9*maximum(output["dTdt"])
        d = output["dTdt"][1]
        while d < MAX90 && k < n
            k = k + 1
            d = output["dTdt"][k]
        end
        output["ind_time_90"] = output["time"][k]

        # find exothermic time
        tstep1 = 1
        tstep2 = 1
        half_T_flag1 = 0
        half_T_flag2 = 0
        # Go into a loop to find two times when temperature is half its maximum
        for (j,tgrad) in enumerate(output["dTdt"])
            if half_T_flag1 == 0
                if tgrad > 0.5*maximum(output["dTdt"])
                    half_T_flag1 = 1
                    tstep1 = j
                end

            elseif half_T_flag2 == 0
                if tgrad < 0.5*maximum(output["dTdt"])
                    half_T_flag2 = 1
                    tstep2 = j
                else
                    tstep2 = 1
                end
            end
        end
    end

    # Exothermic time for CV explosion
    if tstep2 == 1
        println("Error: No pulse in the temperature gradient")
        println("       Your final integration length may be too short,")
        println("       your mixture may be too rich/lean, or something else may be wrong")
        output["exo_time"] = 0
    else
        output["exo_time"] = output["time"][tstep2] - output["time"][tstep1]
    end

    output["gas"] = gas
    return output
end

end #module


text:
export conv, plan_conv, conv_psf, plan_conv_psf

"""
    conv(u, v[, dims])

Convolve `u` with `v` over `dims` dimensions with an FFT based method.
Note, that this method introduces wrap-around artifacts without
proper padding/windowing.

# Arguments
* `u` is an array in real space.
* `v` is the array to be convolved in real space as well.
* Per default `ntuple(+, min(N, M)))` means that we perform the convolution
    over all dimensions of that array which has less dimensions.
    If `dims` is an array with integers, we perform convolution
    only over these dimensions. Eg. `dims=[1,3]` would perform the convolution
    over the first and third dimension. Second dimension is not convolved.

If `u` and `v` are both a real valued array we use `rfft` and hence
the output is real as well.
If either `u` or `v` is complex we use `fft` and output is hence complex.

# Examples
1D with FFT over all dimensions. We choose `v` to be a delta peak.
Therefore convolution should act as identity.
```jldoctest
julia> u = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5
julia> v = [0 0 1 0 0]
1×5 Array{Int64,2}:
 0  0  1  0  0
julia> conv(u, v)
1×5 Matrix{Float64}:
 4.0  5.0  1.0  2.0  3.0
```
2D with FFT with different `dims` arguments.
```jldoctest
julia> u = 1im .* [1 2 3; 4 5 6]
2×3 Matrix{Complex{Int64}}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im
julia> v = [1im 0 0; 1im 0 0]
2×3 Matrix{Complex{Int64}}:
 0+1im  0+0im  0+0im
 0+1im  0+0im  0+0im
julia> conv(u, v)
2×3 Matrix{ComplexF64}:
 -5.0+0.0im  -7.0+0.0im  -9.0+0.0im
 -5.0+0.0im  -7.0+0.0im  -9.0+0.0im
```
"""
function conv(u::AbstractArray{T, N}, v::AbstractArray{D, M}, dims=ntuple(+, min(N, M))) where {T, D, N, M}
    return ifft(fft(u, dims) .* fft(v, dims), dims)
end

function conv(u::AbstractArray{<:Real, N}, v::AbstractArray{<:Real, M}, dims=ntuple(+, min(N, M))) where {N, M}
    return irfft(rfft(u, dims) .* rfft(v, dims), size(u, dims[1]), dims)
end

"""
    conv_psf(u, psf[, dims])
`conv_psf` is a shorthand for `conv(u,ifftshift(psf))`. For examples see `conv`.
"""
function conv_psf(u::AbstractArray{T, N}, psf::AbstractArray{D, M}, dims=ntuple(+, min(N, M))) where {T, D, N, M}
    return conv(u, ifftshift(psf, dims), dims)
end

 # define custom adjoint for conv
 # so far only defined for the derivative regarding the first component
function ChainRulesCore.rrule(::typeof(conv), u::AbstractArray{T, N}, v::AbstractArray{D, M},
                              dims=ntuple(+, min(N, M))) where {T, D, N, M}
    Y = conv(u, v, dims)
    function conv_pullback(barx)
        return NoTangent(), conv(barx, conj(v), dims), NoTangent(), NoTangent()
    end
    return Y, conv_pullback
end

"""
    plan_conv(u, v [, dims])

Pre-plan an optimized convolution for arrays shaped like `u` and `v` (based on pre-plan FFT)
along the given dimenions `dims`.
`dims = 1:ndims(u)` per default.
The 0 frequency of `u` must be located at the first entry.
We return two arguments:
The first one is `v_ft` (obtained by `fft(v)` or `rfft(v)`).
The second return is the convolution function `pconv`.
`pconv` itself has two arguments. `pconv(u, v_ft=v_ft)` where `u` is the object and `v_ft` the v_ft.
This function achieves faster convolution than `conv(u, u)`.
Depending whether `u` is real or complex we do `fft`s or `rfft`s

# Warning
The resulting output of the `pconv` function is a reference to an internal, allocated array.
If you use the `pconv` function for different tasks,
a new call to `pconv` will change the previous result (since the previous result was only a reference, not a new array).

# Examples
```jldoctest
julia> u = [1 2 3 4 5]
1×5 Matrix{Int64}:
 1  2  3  4  5
julia> v = [1 0 0 0 0]
1×5 Matrix{Int64}:
 1  0  0  0  0
julia> v_ft, pconv = plan_conv(u, v);
julia> pconv(u, v_ft)
1×5 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0
julia> pconv(u)
1×5 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0
```
"""
function plan_conv(u::AbstractArray{T, N}, v::AbstractArray{T, M}, dims=ntuple(+, N)) where {T, N, M}
    plan = get_plan(T)
    # do the preplanning step
    P = plan(u, dims)
    u_ft_stor = P * u
    P_inv = inv(P)

    v_ft = fft_or_rfft(T)(v, dims)
    out = similar(u)
    # construct the efficient conv function
    # P and P_inv can be understood like matrices
    # but their computation is fast
    conv(u, v_ft=v_ft) = p_conv_aux!(P, P_inv, u, v_ft, u_ft_stor, out)
    return v_ft, conv
end

"""
    plan_conv_psf(u, psf [, dims]) where {T, N}

`plan_conv_psf` is a shorthand for `plan_conv(u, ifftshift(psf))`. For examples see `plan_conv`.
"""
function plan_conv_psf(u::AbstractArray{T, N}, psf::AbstractArray{T, M}, dims=ntuple(+, N)) where {T, N, M}
    return plan_conv(u, ifftshift(psf, dims), dims)
end

function p_conv_aux!(P, P_inv, u, v_ft, u_ft_stor, out)
    #return P_inv.scale .* (P_inv.p * ((P * u) .* v_ft))
    mul!(u_ft_stor, P, u)
    u_ft_stor .*= v_ft
    mul!(out, P_inv.p, u_ft_stor)
    #out2 = out .* P_inv.scale
    out .*= P_inv.scale
    return out
end

function ChainRulesCore.rrule(::typeof(p_conv_aux!), P, P_inv, u, v, u_ft_stor, out)
    Y = p_conv_aux!(P, P_inv, u, v, u_ft_stor, out)
    function conv_pullback(barx)
        conj_v = let
            if eltype(v) <: Real
                v
            else
                conj(v)
            end
        end
        barx = let
            if typeof(barx) <: FillArrays.Fill
                collect(eltype(u).(barx))
            else
                barx
            end
        end
        ∇ = p_conv_aux!(P, P_inv, barx, conj_v, u_ft_stor, copy(out))
        return NoTangent(), NoTangent(), NoTangent(), ∇, NoTangent(), NoTangent(), NoTangent()
    end
    return Y, conv_pullback
end

"""
    fft_or_rfft(T)
Small helper function to decide whether a real
or a complex valued FFT is appropriate.
"""
function fft_or_rfft(::Type{<:Real})
    return rfft
end

function fft_or_rfft(::Type{T}) where T
    return fft
end

"""
    get_plan(T)

Small helper function to decide whether a real
or a complex valued FFT plan is appropriate.
"""
function get_plan(::Type{<:Real})
    return plan_rfft
end

function get_plan(::Type{T}) where T
    return plan_fft
end


text:
# Base.eltype
eltype(::Type{BondGraph}) = AbstractNode
eltype(bg::BondGraph) = AbstractNode

# edgetype
g.edgetype(::Type{BondGraph}) = g.AbstractSimpleEdge{Integer}
g.edgetype(bg::BondGraph) = g.AbstractSimpleEdge{Integer}

# edges
g.edges(bg::BondGraph) = bg.bonds
g.ne(bg::BondGraph) = length(bg.bonds)
g.has_edge(bg::BondGraph, bond::Bond) =  any(b -> b === bond, bg.bonds) # strong equality
g.has_edge(bg::BondGraph, n1::AbstractNode, n2::AbstractNode) = g.has_edge(bg, vertex(n1), vertex(n2))
g.has_edge(bg::BondGraph, s::Int, d::Int) =
    any(b -> g.src(b) === s && g.dst(b) === d, bg.bonds)

# vertices
g.vertices(bg::BondGraph) = vertex.(bg.nodes)
g.nv(bg::BondGraph) = length(bg.nodes)
g.has_vertex(bg::BondGraph, node::AbstractNode) = any(n -> n === node, bg.nodes) # strong equality
g.has_vertex(bg::BondGraph, v::Int) = 1 <= v <= length(bg.nodes)

# inneighbors, outneighbors
g.inneighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[g.inneighbors(bg, vertex(n))]
g.inneighbors(bg::BondGraph, v::Int) = [g.src(b) for b in bg.bonds if g.dst(b) == v]
g.outneighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[g.outneighbors(bg, vertex(n))]
g.outneighbors(bg::BondGraph, v::Int) = [g.dst(b) for b in bg.bonds if g.src(b) == v]
g.all_neighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[g.all_neighbors(bg, vertex(n))]

# is_directed
g.is_directed(::Type{BondGraph}) = true
g.is_directed(bg::BondGraph) = true

# zero
g.zero(::Type{BondGraph}) = BondGraph()
g.zero(bg::BondGraph) = BondGraph()

# src, dst
g.src(b::Bond) = vertex(srcnode(b))
g.dst(b::Bond) = vertex(dstnode(b))

# weights
# TODO

# Mutations
function g.add_vertex!(bg::BondGraph, node::AbstractNode)
    g.has_vertex(bg, node) && return false
    push!(bg.nodes, node)
    set_vertex!(node, g.nv(bg))
    return true
end

function g.rem_vertex!(bg::BondGraph, node::AbstractNode)
    g.has_vertex(bg, node) || return false
    index = vertex(node)
    deleteat!(bg.nodes, index)
    for n in bg.nodes[index:end]
        n.vertex[] -= 1
    end
    return true
end

# add_edge is done in terms of ports so that a port index can be specified higher up
function g.add_edge!(bg::BondGraph, srcport::Port, dstport::Port)
    srcnode = srcport.node
    dstnode = dstport.node
    g.has_edge(bg, srcnode, dstnode) && error("Bond already exists between $srcnode and $dstnode")
    g.has_edge(bg, dstnode, srcnode) && error("Bond already exists between $dstnode and $srcnode")
    new_bond = Bond(srcport, dstport)
    push!(bg.bonds, new_bond)
    srcportidx = new_bond.srcport.index
    dstportidx = new_bond.dstport.index
    updateport!(srcnode, srcportidx)
    updateport!(dstnode, dstportidx)
    srcnode isa Junction && set_weight!(srcnode, srcportidx, -1)
    dstnode isa Junction && set_weight!(dstnode, dstportidx, +1)
    return new_bond
end
function g.add_edge!(bg::BondGraph, srcnode::AbstractNode, dstnode::AbstractNode)
    g.add_edge!(bg, Port(srcnode), Port(dstnode))
end

function g.rem_edge!(bg::BondGraph, node1::AbstractNode, node2::AbstractNode)
    g.has_edge(bg, node1, node2) || g.has_edge(bg, node2, node1) || return
    index = findfirst(b -> node1 in b && node2 in b, bg.bonds)
    deleted_bond = bg.bonds[index]
    deleteat!(bg.bonds, index)

    src_node = srcnode(deleted_bond)
    dst_node = dstnode(deleted_bond)
    srcportidx = deleted_bond.srcport.index
    dstportidx = deleted_bond.dstport.index

    updateport!(src_node, srcportidx)
    updateport!(dst_node, dstportidx)
    src_node isa Junction && set_weight!(src_node, srcportidx, 0)
    dst_node isa Junction && set_weight!(dst_node, dstportidx, 0)
    return deleted_bond
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
using DiffEqBase, OrdinaryDiffEq, DiffEqJump, DiffEqBiological
using StaticArrays, DataStructures
using DiffEqJumpExtensions
using Plots; pyplot()
using BenchmarkTools
using BioSimulator

doPlot      = false
doBenchMark = true

# DNA repression model DiffEqBiological
rs = @reaction_network dtype begin
    k1, DNA --> mRNA + DNA
    k2, mRNA --> mRNA + P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
tf        = 20000.0
params    =  (.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.)
prob      = DiscreteProblem([1,0,0,0], (0.0, tf), params)

# Biosimlator model of the same network
model = Network("Gene Expression with Neg Feedback")
model <= Species("DNA", 1)
model <= Species("mRNA", 0)
model <= Species("P", 0)
model <= Species("DNAR", 0)
model <= Reaction("transcription", .5, "DNA --> mRNA + DNA")
model <= Reaction("translation", 20*log(2.)/120., "mRNA --> mRNA + P")
model <= Reaction("mRNA degradation", (log(2.)/120.), "mRNA --> 0")
model <= Reaction("protein degradation", (log(2.)/600.), "P --> 0")
model <= Reaction("DNA repression", .025, "DNA + P --> DNAR")
model <= Reaction("DNA freedom", 1., "DNAR --> DNA + P")

spec = OrderedDict( :DNA => 1, :mRNA => 2, :P => 3, :DNAR => 4)
scaled_rates = [rate for rate in params]
# reactant_stoch = [
#    [spec[:DNA] 1]',
#    [spec[:mRNA] 1]',
#    [spec[:mRNA] 1]',
#    [spec[:P] 1]',
#    [spec[:DNA] spec[:P]; 1 1],
#    [spec[:DNAR] 1]'
#    ]
# net_stoch = [
#     [spec[:mRNA] 1]',
#     [spec[:P] 1]',
#     [spec[:mRNA] -1]',
#     [spec[:P] -1]',
#     [spec[:DNA] spec[:P] spec[:DNAR]; -1 -1 1],
#     [spec[:DNA] spec[:P] spec[:DNAR]; 1 1 -1]
#     ]
    reactant_stoch = [
        [spec[:DNA] => 1],
        [spec[:mRNA] => 1],
        [spec[:mRNA] => 1],
        [spec[:P] => 1],
        [spec[:DNA] => 1, spec[:P] => 1],
        [spec[:DNAR] => 1]
        ]
     net_stoch = [
         [spec[:mRNA] => 1],
         [spec[:P] => 1],
         [spec[:mRNA] => -1],
         [spec[:P] => -1],
         [spec[:DNA] => -1, spec[:P] => -1, spec[:DNAR] => 1],
         [spec[:DNA] => 1, spec[:P] => 1, spec[:DNAR] => -1]
         ]

# vectors of SMatrices seem to be really slow here...
# reactant_stoch = [
#     SMatrix{2,1}([spec[:DNA] 1]),
#     SMatrix{2,1}([spec[:mRNA] 1]),
#     SMatrix{2,1}([spec[:mRNA] 1]),
#     SMatrix{2,1}([spec[:P] 1]),
#     SMatrix{2,2}([spec[:DNA] spec[:P]; 1 1]),
#     SMatrix{2,1}([spec[:DNAR] 1])
#     ]
# net_stoch = [
#        SMatrix{2,1}([spec[:mRNA] 1]),
#        SMatrix{2,1}([spec[:P] 1]),
#        SMatrix{2,1}([spec[:mRNA] -1]),
#        SMatrix{2,1}([spec[:P] -1]),
#        SMatrix{2,3}([spec[:DNA] spec[:P] spec[:DNAR]; -1 -1 1]),
#        SMatrix{2,3}([spec[:DNA] spec[:P] spec[:DNAR]; 1 1 -1])
#        ]

# DiffEqJump methods to benchmark
#methods = (Direct2(), DirectVEC(), DiffEqJumpExtensions.FRM(), FRMVEC())
methods = (DirectVEC(), DirectMA(scaled_rates, reactant_stoch, net_stoch))

# plotting solutions
if doPlot
    result  = simulate(model, algorithm=BioSimulator.SSA, time=tf, epochs=25000, trials=1)
    gui( plot(result.t_index, result[1,:,1], reuse=false) )

    for method in methods
        jump_prob = JumpProblem(prob, method, rs)
        sol = solve(jump_prob, SSAStepper())
        gui(plot!(sol.t, sol[3,:]))
    end

end

if doBenchMark
    # number of save times to look at
    num_save_vec = 25000 #10.^(3:6)

    # exact methods
    for method in methods
        # println("Solving with method: ", typeof(method), ", using FunctionMap")
        # jump_prob = JumpProblem(prob, method, rs)
        # @btime solve($jump_prob, FunctionMap())

        println("Solving with method: ", typeof(method), ", using SSAStepper")
        jump_prob = JumpProblem(prob, method, rs)
        @btime solve($jump_prob, SSAStepper())
    end
    println()

    # methods with different numbers of points
    for num_save in num_save_vec
        dtsave = tf / num_save

        for method in methods
            # println("Solving with method: ", typeof(method), ", using FunctionMap and ", num_save, " points")
            # jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
            # @btime solve($jump_prob, FunctionMap(), saveat=$dtsave)

            println("Solving with method: ", typeof(method), ", using SSAStepper and ", num_save, " points")
            jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
            @btime solve($jump_prob, SSAStepper(), saveat=$dtsave)
        end

        # BioSimulator:
        println("BioSimulator, Direct, saving at ", num_save, " points:")
        @btime simulate($model, BioSimulator.SSA, time=$tf, epochs=$num_save, trials=1);

        println("BioSimulator, FRM, saving at ", num_save, " points:")
        @btime simulate($model, BioSimulator.FRM, time=$tf, epochs=$num_save, trials=1);

        println()
    end

end

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

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
import Base: peek

const whitespace = " \t"

function skipwhitespace(io::IO; newlines = true)
  while !eof(io) && (peek(io) in whitespace || (newlines && peek(io) == '\n'))
    read(io, Char)
  end
  return io
end

function startswith(stream::IO, s::AbstractString; eat = true, padding = false)
  start = position(stream)
  padding && skipwhitespace(stream)
  result = true
  for char in s
    !eof(stream) && read(stream, Char) == char ||
      (result = false; break)
  end
  !(result && eat) && seek(stream, start)
  return result
end

function startswith(stream::IO, c::Char; eat = true)
  if peek(stream) == c
    eat && read(stream, Char)
    return true
  else
    return false
  end
end

function startswith{T<:AbstractString}(stream::IO, ss::Vector{T}; eat = true)
  for s in ss
    startswith(stream, s, eat = eat) && return true
  end
  return false
end

function startswith(stream::IO, r::Regex; eat = true, padding = false)
  @assert beginswith(r.pattern, "^")
  start = position(stream)
  padding && skipwhitespace(stream)
  line = chomp(readline(stream))
  seek(stream, start)
  m = match(r, line)
  m == nothing && return ""
  eat && @dotimes length(m.match) read(stream, Char)
  return m.match
end

function peekbehind(stream::IO, offset = 0)
  c = '\0'
  skip(stream, offset)
  if position(stream) > 0
    skip(stream, -1)
    c = read(stream, Char)
  end
  skip(stream, -offset)
  return c
end


text:
using BigArrays
using BigArrays.BinDicts
using Test
using RealNeuralNetworks
using RealNeuralNetworks.Manifests
using RealNeuralNetworks.NodeNets
using RealNeuralNetworks.SWCs
using JSON

const cellId = 77497
const VOXEL_SIZE = (80,80,45)

@testset "test basic data structure" begin
    h = BinDict("/neuroglancer/zfish_v1/consensus-20180123/mesh_mip_4")
    manifest = String(h["$(cellId):0"])
    manifest = JSON.parse(manifest; dicttype=Dict{Symbol, Any})

    str = split(manifest[:fragments][1], ":")[end]
    @show str
    range = BigArrays.Indexes.string2unit_range( str )
    @test range == [2865:3376, 2017:2528, 16401:16912]
end

@testset "test manifest iteration" begin
    manifest = Manifest("/neuroglancer/zfish_v1/consensus-20180123/mesh_mip_4",
                        "$(cellId):0", "/neuroglancer/zfish_v1/consensus-20180123/40_40_45")
    println("trace a cell...")
    @time nodeNet = Manifests.trace(manifest, cellId)
    swc = SWC( nodeNet )
    SWCs.stretch_coordinates!(swc, VOXEL_SIZE)
    SWCs.save(swc, "/tmp/$(cellId).swc")
end


text:
@testset "Rosenbrock" begin

    function test_result(result::Evolutionary.EvolutionaryOptimizationResults, N::Int, tol::Float64)
        fitness = minimum(result)
        extremum = Evolutionary.minimizer(result)
        @test sum(abs2,extremum.-ones(N)) ≈ 0.0 atol=tol
        @test fitness ≈ 0.0 atol=tol
    end

    # Objective function
    rosenbrock(x::AbstractVector) = (1.0 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

    # Parameters
    N = 2

    # Testing: (15,3/(+,)100)-ES
    settings = [
        :isotropic=>(gaussian, gaussian, IsotropicStrategy(N)),
        :isotropic=>(cauchy, gaussian, IsotropicStrategy(N)),
        :anisotropic=>(gaussian, gaussian, AnisotropicStrategy(N))
    ]
    selections = [:plus, :comma]
    @testset "ES settings" for (sn,ss) in settings, sel in selections
        result = Evolutionary.optimize( rosenbrock, (() -> rand(N)),
            ES(
                initStrategy = ss[3],
                recombination = average, srecombination = average,
                mutation = ss[1], smutation = ss[2],
                μ = 10, ρ = 3, λ = 100, selection=sel
            ), Evolutionary.Options(iterations=1000, successive_f_tol=25)
        )
        println("(15/3$(sel == :plus ? "+" : ",")100)-ES:$(sn) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
        test_result(result, N, sel == :plus ? 0.1 : 0.5)
    end

    # Testing: CMA-ES
    result = Evolutionary.optimize(rosenbrock, (() -> rand(Float32,N)), CMAES(μ = 5, λ = 100))
    println("(5/5,100)-CMA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-2)

    # Testing: GA
    result = Evolutionary.optimize(rosenbrock, (() -> rand(N)),
        GA(
            populationSize = 100,
            ɛ = 0.1,
            selection = rouletteinv,
            crossover = intermediate(0.25),
            mutation = domainrange(fill(0.5,N))
        ))
    println("GA(p=100,x=0.8,μ=0.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)

    # Testing: DE
    result = Evolutionary.optimize(rosenbrock, (() -> rand(N)), DE(populationSize = 100))
    println("DE/rand/1/bin(F=1.0,Cr=0.5) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)

end


text:
function optimize(solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
                    transit_map::InputGraph, planarity_constraints::Integer,
                    elen_factor::Integer = 3,
                    bc_factor::Integer = 3,
                    rpos_factor::Integer = 1)
    input_graph = transit_map
    station_list = stations(input_graph)
    edge_list = edges(input_graph)
    faces = edge_faces(input_graph)
    model = Model(solver = solver)
    variables = build_model!(model, input_graph, faces, planarity_constraints,
                                    elen_factor, elen_factor, rpos_factor)
    status = solve(model)
    x_val = getvalue(variables.x)
    y_val = getvalue(variables.y)
    function map_station(i)
        station = station_list[i]
        coord = EuclideanCoordinate(round(x_val[i], 6), round(y_val[i], 6))
        EuclideanStation(station.id, coord, station.label, station.is_dummy)
    end
    euclidian_stations = map(map_station, 1:nstations(input_graph))
    function map_edge(edge)
        idx_from = first(indexin([edge.from], station_list))
        idx_to = first(indexin([edge.to], station_list))
        Edge(euclidian_stations[idx_from],
             euclidian_stations[idx_to], edge.line)
    end
    euclidian_edges = map(map_edge, edges(input_graph))
    TransitMapLayout(euclidian_stations,
        euclidian_edges,
        input_graph.lines,
        faces)
end

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
"""
File defining Trajectory, SpatialTrajectory and related functions
"""

"""
A struct to store all data needed to iterate. principal_vectors have to be orthonormal
"""
struct Trajectory
    data_set
    n :: Int
    k :: Int
    error :: Float64
    indices :: Array{Int,1}
    principal_vectors :: Array{Float64,2}
end

    function show(io :: IO, traj :: Trajectory)
        return print(io, "error: $(traj.error)\nindices: $(traj.indices)")
    end

    "start with indices"
    function Trajectory(data_set, n, k, indices :: AbstractArray{I}) where I <: Integer
        @assert length(indices) == n
        error, svd = truncated_svd(hcat(data_set[indices]...), k)
        Trajectory(data_set, n, k, error, sort(indices), svd.U)
    end

    "pick the points closest to the span of vectors"
    function Trajectory(data_set, n, k, vectors :: Array{T, 2}) where T <: Number
        vectors = Matrix(qr(vectors).Q) # to orthonormalize principal_vectors
        indices, S = find_n_closest_pts(data_set, vectors, n)
        error, svd = truncated_svd(hcat(S...), k)
        Trajectory(data_set, n, k, error, sort(indices), svd.U)
    end

    "pick k points from the dataset at random and pick the rest n-k points to be the closest to the span of the k random points"
    function Trajectory(data_set, n, k)
        principal_vectors = hcat(sample(data_set, k)...)
        Trajectory(data_set, n, k, principal_vectors)
    end

    "perform one step: find n closest points and compute principal vectors. Return a new Trajectory object"
    function L2_step(trajectory)
        @unpack data_set, n, k, principal_vectors = trajectory
        indices, S = find_n_closest_pts(data_set, principal_vectors, n)
        error, svd = truncated_svd(hcat(S...), k)
        Trajectory(data_set, n, k, error, sort(indices), svd.U)
    end

    "perform one step: find n closest points and compute principal vectors. Return a new Trajectory object"
    function L2_step_local(trajectory)
        @unpack data_set, n, k, principal_vectors, indices = trajectory
        worst_index = maximum(i -> (distance_squared(data_set[i], principal_vectors), i), indices)[2]
        indices_minus_worst = setdiff(indices, worst_index)
        indices_to_search = setdiff(1:length(data_set), indices_minus_worst)
        new_ind = minimum(i -> (distance_squared(data_set[i], principal_vectors), i), indices_to_search)[2]
        indices = union(indices_minus_worst, new_ind)
        S = data_set[indices]
        error, svd = truncated_svd(hcat(S...), k)
        d1 = distance_squared( data_set[worst_index], principal_vectors)
        d2 = distance_squared( data_set[new_ind], principal_vectors)
        Trajectory(data_set, n, k, error, sort(indices), svd.U)
    end

"""
A struct to store all data needed to iterate, including a graph. principal_vectors have to be orthonormal
"""
struct SpatialTrajectory
    graph
    data_set
    n :: Int
    k :: Int
    λ :: Float64
    error :: Float64
    indices :: Array{Int,1}
    principal_vectors :: Array{Float64,2}
end

    function show(io :: IO, traj :: SpatialTrajectory)
        return print(io, "error: $(traj.error)\nindices: $(traj.indices)\nprincipal vectors: $(traj.principal_vectors)")
    end


text:
using HDF5
using Plots

mode = 2
folder_dir = "./Data/default_c=4_5/"

if mode==1
    for i in 0.0:0.2:1.0
        path = folder_dir*"break_rules="*string(i)*".h5"
        h5open(path, "r") do file
            data = read(file, "data")
            #plot(data[1,:], [data[2,:], data[3,:], data[4,:], data[5,:], data[6,:]])
            plot(data[1,:], data[2,:])
            savefig(folder_dir*"break_rules="*string(i)*".svg")
        end
    end
end

if mode==2
    path = folder_dir
    data = Any[[], [], [], [], [], []]
    for i in 0.0:0.2:1.0
        h5open(path*"break_rules="*string(i)*".h5", "r") do file
            global data
            data[Int(i*5+1)] = read(file, "data")
        end
    end
    plot(data[1][1,:], data[1][2,:])
    for i in 2:6
        plot!(data[i][1,:], data[i][2,:])
    end
    savefig(path*"summary.svg")
end


text:
module PackVec

export @pack!, @unpack!

using Requires

if VERSION < v"0.7.0-"
    using Compat
end

"""
    @pack!(vec, A, B, C, ...)

Pack the objects `A`, `B`, `C` etc. into `vec`.

Each object requires a `pack!` method. It returns the number of objects inserted.
"""
macro pack!(vec, args...)
    blk = Expr(:block)
    push!(blk.args, :(offset = 0))
    for X in args
        push!(blk.args, :(offset = pack!($(esc(vec)), offset, $(esc(X)))))
    end
    blk
end

"""
    @unpack!(vec, A[::spec], B[::spec], ...)

Unpack the vector into objects named `A`, `B`, etc.

If a `spec` is defined, it creates a new object of that type by calling
`unpack`. Otherwise it requires that the object already exist and mutates it by calling
`unpack!`.
"""
macro unpack!(vec, args...)
    blk = Expr(:block)
    push!(blk.args, :(offset = 0))
    for ex in args
        if ex isa Expr && ex.head == :(::)
            X = ex.args[1]
            spec = ex.args[2]
            push!(blk.args, :(($(esc(X)), offset) = unpack($(esc(vec)), offset, $(esc(spec)))))
        else
            X = ex
            push!(blk.args, :(offset = unpack!($(esc(vec)), offset, $(esc(X)))))
        end
    end
    push!(blk.args, :(offset))
    blk
end

"""
    nextoffset = pack!(vec, offset, obj)

This is called by `@pack!`. It should pack the object into `vec`, starting at `offset+1`,
returning the index of the final insertion.
"""
function pack! end

"""
    nextoffset = unpack!(vec, offset, obj)

This is called by `@unpack!`. It should unpack `vec` (starting at `offset+1`) into `obj`,
returning the next offset.

See also `unpack`.
"""
function unpack! end

"""
    obj, nextoffset = unpack(vec, offset, T)

This is called by `@unpack!`. It should unpack `vec` (starting at `offset+1`) into a new
object of type `T`, returning a tuple of the object and next offset.

See also `unpack`.
"""
function unpack end

function pack!(vec, offset, obj::Real)
    vec[offset += 1] = obj
    return offset
end
function unpack(vec, offset, ::Type{T}) where {T<:Real}
    obj = T(vec[offset += 1])
    return obj, offset
end

function pack!(vec, offset, obj::AbstractArray)
    n = length(obj)
    copyto!(vec, offset+1, obj, 1, n)
    return offset+n
end
function unpack!(vec, offset, obj::AbstractArray{T}) where {T<:Real}
    n = length(obj)
    copyto!(obj, 1, vec, offset+1, n)
    return offset+n
end

pack!(vec, offset, obj::Diagonal) = pack!(vec, offset, obj.diag)
unpack!(vec, offset, obj::Diagonal) = unpack!(vec, offset, obj.diag)

function pack!(vec, offset, obj::UpperTriangular)
    n = size(obj,1)
    for i = 1:n
        copyto!(vec, offset+1, obj.data, (i-1)*n+1, i)
        offset += i
    end
    return offset
end
function unpack!(vec, offset, obj::UpperTriangular)
    n = size(obj,1)
    for i = 1:n
        copyto!(obj.data, (i-1)*n+1, vec, offset+1, i)
        offset += i
    end
    return offset
end

function pack!(vec, offset, obj::LowerTriangular)
    n = size(obj,1)
    for i = 1:n
        ii = n-i+1
        copyto!(vec, offset+1, obj, (i-1)*n+i, ii)
        offset += ii
    end
    return offset
end
function unpack!(vec, offset, obj::LowerTriangular)
    n = size(obj,1)
    for i = 1:n
        ii = n-i+1
        copyto!(obj, (i-1)*n+i, vec, offset+1, ii)
        offset += ii
    end
    return offset
end

function pack!(vec, offset, obj::Symmetric)
    if obj.uplo == 'U'
        pack!(vec, offset, UpperTriangular(obj.data))
    else
        pack!(vec, offset, LowerTriangular(obj.data))
    end
end
function unpack!(vec, offset, obj::Symmetric)
    if obj.uplo == 'U'
        unpack!(vec, offset, UpperTriangular(obj.data))
    else
        unpack!(vec, offset, LowerTriangular(obj.data))
    end
end

function pack!(vec, offset, obj::LinAlg.Cholesky)
    if obj.uplo == 'U'
        pack!(vec, offset, UpperTriangular(obj.factors))
    else
        pack!(vec, offset, LowerTriangular(obj.factors))
    end
end
function unpack!(vec, offset, obj::LinAlg.Cholesky)
    if obj.uplo == 'U'
        unpack!(vec, offset, UpperTriangular(obj.factors))
    else
        unpack!(vec, offset, LowerTriangular(obj.factors))
    end
end

@require StaticArrays begin
    using StaticArrays
    function unpack(vec, offset, ::Type{SVector{N}}) where {N}
        obj = SVector{N}(@view vec[offset+1:offset+N])
        return obj, offset+N
    end
    function unpack(vec, offset, ::Type{SMatrix{M,N}}) where {M,N}
        obj = SMatrix{M,N}(@view vec[offset+1:offset+M*N])
        return obj, offset+M*N
    end
end

end # module

max_stars_repo_name: simonbyrne/PackVec


text:
#

using Gadfly, Colors
using Distributions
using FFTW
using DSP

p1 = Normal(-pi/2, 0.5)
p2 = Normal(pi/2, 0.5)
p12_t = Normal(0.0, 1.0/sqrt(8.0))

x = range(-pi, stop=pi, length=1024)

g1 = (x)->pdf(p1,x)
g2 = (x)->pdf(p2,x)
g12_t = (x)->pdf(p12_t,x)

y1 = g1.(x)
y2 = g2.(x)
y12_t = g12_t.(x)

y12_t .*= 0.5
temp = fftshift(y12_t)
y12_t .+= temp

plot(y=y12_t, Geom.line)

Y12_t = fftshift(fft(y12_t))

pl = Gadfly.plot(
layer(y=abs.(Y1), Geom.line, Theme(default_color=colorant"red")),
layer(y=abs.(Y2), Geom.line, Theme(default_color=colorant"blue")),
layer(y=abs.(Y12_t), Geom.line, Theme(default_color=colorant"magenta"))
)

pl |> SVG("/tmp/test.svg",100cm,100cm)

y12_tr = ifft(ifftshift(Y12_t))

plot(y=abs.(y12_tr), Geom.line)

plot(
layer(x=x, y=y1, Geom.line, Theme(default_color=colorant"red")),
layer(x=x, y=y2, Geom.line, Theme(default_color=colorant"blue"))
)

using Makie

scene = lines(x, abs.(y12_tr), color = :blue)
# scatter!(scene, x, y1, color = :red, markersize = 0.1)

Y1 = fftshift(fft(y1))
Y2 = fftshift(fft(y2))
Y12_t = fft(y12_t)

plot(y=abs.(Y1), Geom.line)
plot(y=abs.(Y2), Geom.line)

Y12 = conv((Y1), (Y2))

st = 1
off = 2046
y12 = ifft(ifftshift(Y12[st:(st+off)]))

plot(y=abs.(y12), Geom.line)

# plot(x=real.(Y1), y=imag.(Y1), Geom.point)
# plot(x=real.(Y2), y=imag.(Y2), Geom.point)

plot(y=abs.(Y12), Geom.line)

##

#YY1 = [Y1; Y1; Y1]
YY2 = [ifftshift(Y2); ifftshift(Y2); ifftshift(Y2)]

plot(y=abs.(YY1), Geom.line)

YY12 = conv(ifftshift(Y1), (YY2))

plot(y=abs.(YY12), Geom.line)

st = 1024
off = 2047
yy12 = ifft(ifftshift(YY12[st:(st+off)]))

pl = plot(y=abs.(yy12), Geom.line)

##

using Cairo, Fontconfig
using Gadfly

pl |> PDF("/tmp/test.pdf", 30cm, 20cm)

run(`evince /tmp/test.pdf`)

##

plot(x=real.(Y12_t), y=imag.(Y12_t), Geom.point)

# brute force product
y12_d = y1 .* y2
y12_d ./= sum(y12_d)*(x[2]-x[1])

plot(
layer(x=x, y=y1, Geom.line, Theme(default_color=colorant"red")),
layer(x=x, y=y2, Geom.line, Theme(default_color=colorant"blue")),
layer(x=x, y=y12_d, Geom.line, Theme(default_color=colorant"magenta"))
)

# calculating the square using convolution

y1sqrt = sqrt.(y1)

Y1sqrt = fft(y1sqrt)

Y11 = conv(Y1sqrt, conj(Y1sqrt))

##

st = 1
off = 512
y11 = ifft(Y11[st:(st+off)])

plot(
layer(y=y1, Geom.line, Theme(default_color=colorant"red")),
layer(y=abs.(y11)./(sum(abs.(y11))*(x[2]-x[1])), Geom.line, Theme(default_color=colorant"blue"))
# layer(x=x, y=y12_d, Geom.line, Theme(default_color=colorant"magenta"))
)

## Linear convolution in Euclid space

xx = 0:127;

Y = (x::Real; mu::Real=0.0) -> pdf(Normal(mu, 1.0), x)

y1 = Y.(xx, mu=10)
y2 = Y.(xx, mu=20)

fY1 = fft(y1)
fY2 = fft(y2)

fYY = fY1.*fY2
yy = abs.(ifft(fYY))

plot(
  layer(
    x=xx,
    y=y1,
    Geom.line
  ),
  layer(
    x=xx,
    y=y2,
    Geom.line,
    Theme(default_color=colorant"deepskyblue")
  ),
  layer(
    x=xx,
    y=yy,
    Geom.line,
    Theme(default_color=colorant"magenta")
  )
)

## Test with Sine function

x = range(-7.5, stop=7.5, length=2048)
y = sin.(10.0*x) + sin.(100.0*x)

Y = fft(y)

Ys = fftshift(Y)

plot(y=abs.(Ys), Geom.line)

#


text:
#=
root_finding:
- Author: Berk
- Date: 2020-08-07
Root-finding methods sampling nonlinear functions
=#
"""
    normalized_data(bbl::BlackBoxLearner)
    normalized_data(X::DataFrame)

Normalizes and returns data (0-1) by lower and upper bounds.
If used on DataFrame, returns lbs and ubs as well.
"""
function normalized_data(bbl::BlackBoxLearner)
    bounds = get_bounds(bbl.vars)
    vks = string.(bbl.vars)
    lbs = [minimum(bounds[key]) for key in bbl.vars]
    ubs = [maximum(bounds[key]) for key in bbl.vars]
    divisor = ones(length(lbs))
    for i = 1:length(lbs)
        if ubs[i] - lbs[i] != 0
            divisor[i] =  ubs[i] - lbs[i]
        end
    end
    normalized_X = reduce(hcat,[(bbl.X[:, vks[i]] .- lbs[i]) ./divisor[i] for i=1:length(lbs)]);
    return normalized_X
end

function normalized_data(X::DataFrame)
    lbs = [minimum(col) for col in eachcol(X)]
    ubs = [maximum(col) for col in eachcol(X)]
    normalized_X = reduce(hcat,[(X[:, i] .- lbs[i]) ./ (ubs[i] - lbs[i]) for i=1:length(lbs)]);
    return normalized_X, lbs, ubs
end

function build_knn_tree(bbl::BlackBoxLearner)
    """ Builds an efficient KNN tree over existing data."""
    X = normalized_data(bbl); #TODO: optimize calls to normalize data.
    kdtree = KDTree(X', reorder = false);
    bbl.knn_tree = kdtree;
    return
end

function find_knn(bbl::BlackBoxLearner; k::Int64 = 10)
    """ Returns idxs and Euclidian distances of KNNs of each data point. """
    idxs, dists = knn(bbl.knn_tree, bbl.knn_tree.data, k+1);
    return idxs, dists
end

""" Classifies KNN domains by feasibility. """
function classify_patches(bbc::BlackBoxLearner, idxs::Array, threshold = 0)
    arr = []
    for idx in idxs
        signs = [bbc.Y[i] for i in idx];
        if all(signs .>= threshold)
            push!(arr, "feas")
        elseif all(signs .< threshold)
            push!(arr, "infeas")
        else
            push!(arr, "mixed")
        end
    end
    return arr
end

function secant_method(X::DataFrame, Y::Array)
    """ Generates samples estimated to be near a zero from X and Y data with mixed feasibility. """
    feas = Y .>= 0;
    np = DataFrame(propertynames(X) .=> [Float64[] for i in propertynames(X)])
    for i = 1:length(feas)
        for j = i+1:length(feas)
            if feas[i] + feas[j] == 1
                 if !isinf(Y[j]) && !isinf(Y[i]) # For non-infinite outputs
                    push!(np, Array(X[j,:]) - (Array(X[j,:]) - Array(X[i,:])) ./ (Y[j] - Y[i]) * Y[j])
                 else # For Infs, try midpoint.
                    push!(np, Array(X[j,:]) - (Array(X[j,:]) - Array(X[i,:])) ./ 2.)
                 end
            end
        end
    end
    return np
end

"""
    classify_curvature(bbl::BlackBoxLearner, idxs = collect(1:size(bbl.X, 1)))

Classify curvature of KNN patches, over given point indices.
"""
function classify_curvature(bbl::BlackBoxLearner, idxs = collect(1:size(bbl.X, 1)))
    size(bbl.X, 1) > 0 || throw(OCTHaGOnException("Must sample Learner $(bbl.name) first."))
    build_knn_tree(bbl);
    knn_idxs, dists = find_knn(bbl, k=length(bbl.vars) + 1);
    if any(ismissing.(bbl.gradients[idxs,1])) # Only need to check the first element
        @info("Missing gradient information was updated.")
        miss_idxs = [idx for idx in idxs if ismissing(bbl.gradients[idx,1])]
        update_gradients(bbl, miss_idxs)
    end
    thresh = 1e-10 * (bbl.max_Y - bbl.min_Y)
    if bbl isa BlackBoxRegressor
        for i in idxs
            diffs = [Array(bbl.X[i, :]) - Array(bbl.X[j, :]) for j in knn_idxs[i]]
            center_grad = Array(bbl.gradients[i,:])
            under_offsets = [-dot(center_grad, differ) for differ in diffs]
            actual_offsets = bbl.Y[knn_idxs[i]] .- bbl.Y[i]
            if all(actual_offsets - under_offsets .>= -thresh)
                bbl.curvatures[i] = 1
            elseif all(actual_offsets - under_offsets .<= thresh)
                bbl.curvatures[i] = -1
            else
                bbl.curvatures[i] = 0
            end
        end
    else
        for i in idxs
            diffs = [Array(bbl.X[i, :]) - Array(bbl.X[j, :]) for j in knn_idxs[i]]
            center_grad = Array(bbl.gradients[i,:])
            under_offsets = [-dot(center_grad, differ) for differ in diffs]
            actual_offsets = bbl.Y[knn_idxs[i]] .- bbl.Y[i]
            if all(actual_offsets - under_offsets .<= thresh)
                bbl.curvatures[i] = 1
            elseif all(actual_offsets - under_offsets .>= -thresh)
                bbl.curvatures[i] = -1
            else
                bbl.curvatures[i] = 0
            end
        end
    end
end