text:
using SimpleWeightedGraphs

@testset "SimpleWeightedGraphs" begin
    @info("Ignore warnings relating to adding and removing vertices and edges")
    adjmx1 = [0 1 0; 1 1 1; 0 1 0] # SimpleWeightedGraph
    adjmx2 = [0 1 0; 1 1 1; 1 1 0] # SimpleWeightedDiGraph
    # specific concrete generators - no need for loop
    @test @inferred(eltype(SimpleWeightedGraph())) == Int
    @test @inferred(eltype(SimpleWeightedGraph(adjmx1))) == Int
    @test_throws ErrorException SimpleWeightedGraph(adjmx2)

    @test @inferred(ne(SimpleWeightedGraph(path_digraph(5)))) == 4
    @test @inferred(ne(SimpleWeightedGraph(adjmx1))) == 3
    @test @inferred(!is_directed(SimpleWeightedGraph))

    @test @inferred(eltype(SimpleWeightedDiGraph())) == Int
    @test @inferred(eltype(SimpleWeightedDiGraph(adjmx2))) == Int
    @test @inferred(ne(SimpleWeightedDiGraph(path_graph(5)))) == 8
    @test @inferred(ne(SimpleWeightedDiGraph(adjmx2))) == 6
    @test @inferred(is_directed(SimpleWeightedDiGraph))

    for gbig in [SimpleWeightedGraph(0xff), SimpleWeightedDiGraph(0xff)]
        @test @inferred(!add_vertex!(gbig))    # overflow
        @test @inferred(!add_vertices!(gbig, 10))
    end

    gdx = SimpleWeightedDiGraph(path_digraph(4))
    gx = SimpleWeightedGraph()
    for g in testgraphs(gx)
        T = eltype(g)
        U = weighttype(g)
        @test sprint(show, g) == "{0, 0} undirected simple $T graph with $U weights"
        @inferred(add_vertices!(g, 5))
        @test sprint(show, g) == "{5, 0} undirected simple $T graph with $U weights"
    end
    gx = SimpleWeightedDiGraph()
    for g in testdigraphs(gx)
        T = eltype(g)
        U = weighttype(g)
        @test sprint(show, g) == "{0, 0} directed simple $T graph with $U weights"
        @inferred(add_vertices!(g, 5))
        @test sprint(show, g) == "{5, 0} directed simple $T graph with $U weights"
    end

    gx = SimpleWeightedGraph(path_graph(4))

    gc = copy(gx)
    @test_logs (:warn, "Note: adding edges with a zero weight to this graph type has no effect.") add_edge!(gc, 4, 1, 0.0)
    @test !(add_edge!(gc, 4, 1, 0.0))

    for g in testgraphs(gx)
        @test @inferred(vertices(g)) == 1:4
        @test SimpleWeightedEdge(2,3) in edges(g)
        @test @inferred(nv(g)) == 4
        @test @inferred(outneighbors(g,2)) == inneighbors(g,2) == neighbors(g,2)
        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(has_edge(g, 3, 2))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4=>1)) && gc == SimpleWeightedGraph(cycle_graph(4))
        @test @inferred(has_edge(gc, 4=>1)) && has_edge(gc, 0x04=>0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4,1))) && gc == SimpleWeightedGraph(cycle_graph(4))
        @test @inferred(has_edge(gc, (4,1))) && has_edge(gc, (0x04, 0x01))
        gc = copy(g)
        @test add_edge!(gc, 4, 1) && gc == SimpleWeightedGraph(cycle_graph(4))

        @test @inferred(inneighbors(g, 2)) == @inferred(outneighbors(g, 2)) == @inferred(neighbors(g,2)) == [1,3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleWeightedGraph{eltype(g), weighttype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(outneighbors(g,1)) == eltype(nv(g))
        T = @inferred(eltype(g))
        U = @inferred(weighttype(g))
        @test @inferred(nv(SimpleWeightedGraph{T, U}(6))) == 6

        # @test @inferred(eltype(SimpleWeightedGraph(T))) == T
        @test @inferred(eltype(SimpleWeightedGraph{T, U}(adjmx1))) == T

        ga = SimpleWeightedGraph(10)
        @test @inferred(eltype(SimpleWeightedGraph{T, U}(ga))) == T

        for gd in testdigraphs(gdx)
            T2 = eltype(gd)
            @test @inferred(eltype(SimpleWeightedGraph(gd))) == T2
        end

        @test @inferred(edgetype(g)) == SimpleWeightedGraphEdge{T, U}
        @test @inferred(copy(g)) == g
        @test @inferred(!is_directed(g))

        e = first(edges(g))
        @test @inferred(has_edge(g, e))
    end

    gdx = SimpleWeightedDiGraph(path_digraph(4))

    gc = copy(gdx)
    @test_logs (:warn, "Note: adding edges with a zero weight to this graph type has no effect.") add_edge!(gc, 4, 1, 0.0)
    @test !(add_edge!(gc, 4, 1, 0.0))

    for g in testdigraphs(gdx)
        @test @inferred(vertices(g)) == 1:4
        @test SimpleWeightedEdge(2,3) in edges(g)
        @test !(SimpleWeightedEdge(3,2) in edges(g))
        @test @inferred(nv(g)) == 4
        @test outneighbors(g, 2) == [3]
        @test inneighbors(g, 2) == [1]

        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(!has_edge(g, 3, 2))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4=>1)) && gc == SimpleWeightedDiGraph(cycle_digraph(4))
        @test @inferred(has_edge(gc, 4=>1)) && has_edge(gc, 0x04=>0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4,1))) && gc == SimpleWeightedDiGraph(cycle_digraph(4))
        @test @inferred(has_edge(gc, (4,1))) && has_edge(gc, (0x04, 0x01))
        gc = @inferred(copy(g))
        @test @inferred(add_edge!(gc, 4, 1)) && gc == SimpleWeightedDiGraph(cycle_digraph(4))

        @test @inferred(inneighbors(g, 2)) == [1]
        @test @inferred(outneighbors(g, 2)) == @inferred(neighbors(g,2)) == [3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(!rem_edge!(gc, 2, 1))
        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleWeightedDiGraph{eltype(g), weighttype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(@inferred(outneighbors(g,1))) == eltype(nv(g))
        T = @inferred(eltype(g))
        U = @inferred(weighttype(g))
        @test @inferred(nv(SimpleWeightedDiGraph{T, U}(6))) == 6

        # @test @inferred(eltype(SimpleWeightedDiGraph(T))) == T
        @test @inferred(eltype(SimpleWeightedDiGraph{T, U}(adjmx2))) == T

        ga = SimpleWeightedDiGraph(10)
        @test @inferred(eltype(SimpleWeightedDiGraph{T, U}(ga))) == T

        for gu in testgraphs(gx)
            T2 = @inferred(eltype(gu))
            @test @inferred(eltype(SimpleWeightedDiGraph(gu))) == T2
        end

        @test @inferred(edgetype(g)) == SimpleWeightedDiGraphEdge{T, U}
        @test @inferred(copy(g)) == g
        @test @inferred(is_directed(g))

        e = first(@inferred(edges(g)))
        @test @inferred(has_edge(g, e))
    end

    gdx = SimpleWeightedDiGraph(complete_digraph(4))
    for g in testdigraphs(gdx)
        @test rem_vertex!(g, 2)
        @test nv(g) == 3 && ne(g) == 6
    end
    g = SimpleWeightedGraph(complete_graph(3), 3)
    @test sum(weights(g)) == 2 * ne(g) * 3
    @test @inferred(get_weight(g, 1, 2)) == 3

    g = SimpleWeightedDiGraph(path_graph(5), 4.0)
    @test sum(weights(g)) == ne(g) * 4.0

    gx = Graph(4,3)
    for g in testsimplegraphs(gx)
        @test eltype(SimpleWeightedGraph(g)) == eltype(g)
    end

    gx = DiGraph(4,3)
    for g in testsimpledigraphs(gx)
        @test eltype(SimpleWeightedGraph(g)) == eltype(g)
    end

    s = SimpleWeightedGraph(path_graph(5), 2)
    s2 = SimpleWeightedGraph([1,2,3,4], [2,3,4,5], [2,2,2,2])
    @test s == s2

    s = SimpleWeightedDiGraph(path_digraph(5), 2)
    s2 = SimpleWeightedDiGraph([1,2,3,4], [2,3,4,5], [2,2,2,2])
    @test s == s2

    s = SimpleWeightedGraph([10,20,30,40], [20,30,40,50], [2,2,2,2])
    @test size(s.weights) == (50, 50)

    s = SimpleWeightedDiGraph([10,20,30,40], [20,30,40,50], [2,2,2,2])
    @test size(s.weights) == (50, 50)

    s = SimpleWeightedGraph([1,2,1], [2,1,2], [1,1,1]; combine = +)
    @test s.weights[2,1] == s.weights[1,2] == 3

    s = SimpleWeightedDiGraph([1,2,1], [2,1,2], [1,1,1]; combine = +)
    @test s.weights[1,2] == 1
    @test s.weights[2,1] == 2

    s = SimpleWeightedDiGraph([1,2,1], [2,1,2], [1,1,2]; combine = (x,y) -> y)
    @test s.weights[1,2] == 1
    @test s.weights[2,1] == 2

    @test SimpleDiGraph(SimpleWeightedDiGraph(cycle_graph(4))) == SimpleDiGraph(cycle_graph(4))
    @test SimpleGraph(SimpleWeightedGraph(path_graph(5))) == path_graph(5)

    @test SimpleWeightedGraph(cycle_graph(4)) == SimpleWeightedGraph(SimpleWeightedGraph(cycle_graph(4)))
    @test SimpleWeightedDiGraph(cycle_digraph(4)) == SimpleWeightedDiGraph(SimpleWeightedDiGraph(cycle_digraph(4)))

    @test SimpleWeightedDiGraph(Matrix(adjacency_matrix(cycle_digraph(4)))) == SimpleWeightedDiGraph(cycle_digraph(4))
    @test SimpleWeightedDiGraph{Int32}(Matrix(adjacency_matrix(cycle_digraph(4)))) == SimpleWeightedDiGraph{Int32, Float64}(SimpleWeightedDiGraph(cycle_digraph(4)))

    @test SimpleWeightedGraph{Int32}(Matrix(adjacency_matrix(cycle_graph(4)))) == SimpleWeightedGraph{Int32, Float64}(SimpleWeightedGraph(cycle_graph(4)))
    @test SimpleWeightedGraph{Int32}(adjacency_matrix(cycle_graph(4))) == SimpleWeightedGraph{Int32, Float64}(SimpleWeightedGraph(cycle_graph(4)))

    @testset "Typed constructors $T" for T in (UInt8, Int32)
        g = SimpleWeightedGraph(T)
        @test g isa AbstractGraph{T}
        @test g isa SimpleWeightedGraph{T, Float64}
        dg = SimpleWeightedDiGraph(T)
        @test dg isa AbstractGraph{T}
        @test dg isa SimpleWeightedDiGraph{T, Float64}
        for U in (Float16, Float32)
            g = SimpleWeightedGraph(T, U)
            @test g isa AbstractGraph{T}
            @test g isa SimpleWeightedGraph{T, U}
            dg = SimpleWeightedDiGraph(T, U)
            @test dg isa AbstractGraph{T}
            @test dg isa SimpleWeightedDiGraph{T, U}
        end
    end

    @testset "Getting weights" begin
        @testset "Testing $G" for G in (SimpleWeightedGraph, SimpleWeightedDiGraph)
            g = G(3)
            @test g[1, 2, Val{:weight}()] ≈ 0
            @test g[1, 3, Val{:weight}()] ≈ 0
            for e in edges(g)
                @test g[e, Val{:weight}] ≈ 0
                esimple = SimpleEdge(src(e), dst(e))
                @test g[esimple, Val{:weight}] ≈ 0
            end
            @test_throws BoundsError g[3, 4, Val{:weight}()]
            @test_throws MethodError g[1, 2, Val{:wight}()]
            add_edge!(g, 1, 2, 5.0)

            @test g[1, 2, Val{:weight}()] ≈ 5
            if is_directed(G)
                @test g[2, 1, Val{:weight}()] ≈ 0
            else
                @test g[2, 1, Val{:weight}()] ≈ 5
            end
            m = adjacency_matrix(g)
            @test g[2, 1, Val{:weight}()] ≈ g.weights[1, 2]
        end
    end

    @testset "Copying graph copies matrix" begin
        dg = SimpleWeightedDiGraph(cycle_digraph(4))
        dg2 = SimpleWeightedDiGraph(dg)
        @test dg[1, 3, Val{:weight}()] ≈ 0
        @test dg2[1, 3, Val{:weight}()] ≈ 0
        @test add_edge!(dg, 1, 3, 2.5)
        @test dg[1, 3, Val{:weight}()] ≈ 2.5
        @test dg2[1, 3, Val{:weight}()] ≈ 0
        dg3 = SimpleWeightedDiGraph{Int, Float64}(dg)
        @test add_edge!(dg, 1, 4, 3.5)
        @test dg[1, 4, Val{:weight}()] ≈ 3.5
        @test dg3[1, 4, Val{:weight}()] ≈ 0
        g1 = SimpleWeightedGraph{Int, Float64}(dg)
        g2 = SimpleWeightedGraph(dg)
        @test g1 == g2
        @test ne(g1) == 5 # 1-2 1-3 2-3 3-4 4-1
        @test g1[1, 3, Val{:weight}()] ≈ 2.5

        g = SimpleWeightedGraph(cycle_graph(5))
        g2 = SimpleWeightedGraph(g)
        @test g[1, 3, Val{:weight}()] ≈ 0
        @test g2[1, 3, Val{:weight}()] ≈ 0
        @test add_edge!(g, 1, 3, 2.5)
        @test g[1, 3, Val{:weight}()] ≈ 2.5
        @test g2[1, 3, Val{:weight}()] ≈ 0.0
        g3 = SimpleWeightedDiGraph{Int, Float64}(g)
        @test add_edge!(g, 1, 4, 3.5)
        @test g[1, 4, Val{:weight}()] ≈ 3.5
        @test g3[1, 4, Val{:weight}()] ≈ 0

        # copy from undirected to directed
        g = SimpleWeightedGraph(cycle_graph(4), 0.5)
        dg = SimpleWeightedDiGraph(g)
        @test g[1,3,Val{:weight}()] ≈ 0
        add_edge!(g, 1, 3, 6.5)
        @test g[1,3,Val{:weight}()] ≈ 6.5
        @test dg[1,3,Val{:weight}()] ≈ 0

        dg = SimpleWeightedDiGraph(cycle_digraph(4), 0.5)
        @test dg[2, 1, Val{:weight}()] ≈ 0
        add_edge!(dg, 2, 1, 0.6)
        g = SimpleWeightedGraph(dg)
        @test g[1, 2, Val{:weight}()] ≈ 1.1
        @test g[1, 3, Val{:weight}()] ≈ 0
        @test g[2, 3, Val{:weight}()] ≈ 0.5
    end
end


text:
"""
$(TYPEDEF)
Context is a container for storing any contextual information.
It uses Memento pattern to keep track of prior history.  While
the context can be changed, explict save/restore can be used
to save the current context and restore the most recently saved
one.

# Fields
$(TYPEDFIELDS)
"""
struct Context{T}
    "ID of the context"
    id::UInt
    "History of context data"
    history::Stack{T}
    "Call path"
    path::Vector{Symbol}
end

ext: jl


text:
@testset "create_main_index()" begin

    # create temp csvs to test on
    temp_csv1 = "./temp1.tsv"
    temp_csv2 = "./temp2.tsv"
    open(temp_csv1, "w") do io
        write(io, "col1|col1\n1|2")
    end
    open(temp_csv2, "w") do io
        write(io, "col1|col1\n3|4")
    end

    # run test
    main_file = "test_main.tsv"
    create_main_index("./", main_file)
    @test read(main_file, String) == "col1|col1_1\n1|2\n3|4\n"

    # remove temp files
    rm(temp_csv1)
    rm(temp_csv2)
    rm(main_file)

end


text:
fp = dirname(@__FILE__)
@testset "DSGE impulse responses to a pre-specified impact matrix" begin
    matdata = load(joinpath(fp, "../../../reference/test_irfdsge.jld2"))
    irfout = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                               zeros(size(matdata["zz"], 1)), matdata["mmm"],
                               matdata["impact"], 1; accumulate = true, cum_inds = 1)
    nobs = size(matdata["zz"], 1)
    for i = 1:nobs
        @test vec(matdata["aairf"][:, nobs * (i - 1) + 1:nobs * i]) ≈
            vec(irfout[:, :, i])
    end

    irfout_accum1 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = true, cum_inds = 1)
    irfout_nocum1 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = false, cum_inds = 1)
    irfout_accum2 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = true,
                                      cum_inds = 1:2)
    irfout_nocum2 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = false,
                                      cum_inds = 1:2)
    irfout_accum3 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = true,
                                      cum_inds = [1,3])
    irfout_nocum3 = impulse_responses(matdata["TTT"], matdata["RRR"], matdata["zz"],
                                      zeros(size(matdata["zz"], 1)), matdata["mmm"],
                                      matdata["impact"], 2; accumulate = false,
                                      cum_inds = [1,3])
    for i = 1:size(irfout_accum1, 3)
        @test cumsum(irfout_nocum1[1, :, i]) ≈ irfout_accum1[1, :, i]
        @test irfout_nocum1[2:end, :, i] ≈ irfout_accum1[2:end, :, i]
        @test cumsum(irfout_nocum2[1:2, :, i]; dims = 2) ≈ irfout_accum2[1:2, :, i]
        @test irfout_nocum2[end, :, i] ≈ irfout_accum2[end, :, i]
        @test cumsum(irfout_nocum3[[1,3], :, i]; dims = 2) ≈ irfout_accum3[[1,3], :, i]
        @test irfout_nocum3[2, :, i] ≈ irfout_accum3[2, :, i]
    end
end

@testset "Impulse responses to structural shocks identified by a DSGE-VAR" begin
    matdata1 = load(joinpath(fp, "../../../reference/test_irfdsge.jld2"))
    matdata2 = load(joinpath(fp, "../../../reference/test_rotations.jld2"))
    ŷ = DSGE.impulse_responses(matdata1["TTT"], matdata1["RRR"], matdata1["zz"],
                               zeros(size(matdata1["zz"], 1)), matdata1["mmm"],
                               matdata1["impact"].^2, Int(matdata2["k"]),
                               convert(Matrix{Float64}, matdata2["cct_sim"]'),
                               matdata2["sig_sim"], Int(matdata2["qahead"]),
                               vec(matdata2["XXpred"]);normalize_rotation = false,
                               test_shocks =
                               convert(Matrix{Float64}, matdata2["Shocks"]'))

    @test ŷ ≈ matdata2["yypred"]'

    Random.seed!(1793) # drawing shocks, so this shouldn't be the same
    ŷ1 = DSGE.impulse_responses(matdata1["TTT"], matdata1["RRR"], matdata1["zz"],
                                zeros(size(matdata1["zz"], 1)), matdata1["mmm"],
                                matdata1["impact"].^2, Int(matdata2["k"]),
                                convert(Matrix{Float64}, matdata2["cct_sim"]'),
                                matdata2["sig_sim"], Int(matdata2["qahead"]),
                                vec(matdata2["XXpred"]), draw_shocks = true,
                                normalize_rotation = false)
    Random.seed!(1793) # re-seeding should work
    ŷ2 = DSGE.impulse_responses(matdata1["TTT"], matdata1["RRR"], matdata1["zz"],
                                zeros(size(matdata1["zz"], 1)), matdata1["mmm"],
                                matdata1["impact"].^2, Int(matdata2["k"]),
                                convert(Matrix{Float64}, matdata2["cct_sim"]'),
                                matdata2["sig_sim"], Int(matdata2["qahead"]),
                                vec(matdata2["XXpred"]), draw_shocks = true,
                                normalize_rotation = false)
    Random.seed!(1793) # flipping shocks shouldn't yield anything different
    @info "The following warning is expected."
    ŷ3 = DSGE.impulse_responses(matdata1["TTT"], matdata1["RRR"], matdata1["zz"],
                                zeros(size(matdata1["zz"], 1)), matdata1["mmm"],
                                matdata1["impact"].^2, Int(matdata2["k"]),
                                convert(Matrix{Float64}, matdata2["cct_sim"]'),
                                matdata2["sig_sim"], Int(matdata2["qahead"]),
                                vec(matdata2["XXpred"]), draw_shocks = true,
                                flip_shocks = true, normalize_rotation = false)
    @test !(ŷ1 ≈ matdata2["yypred"]')
    @test ŷ1 ≈ ŷ2
    @test ŷ1 ≈ ŷ3
end

@testset "Impulse responses of a VAR using a DSGE as a prior" begin
    jlddata = load(joinpath(fp, "../../../reference/test_dsgevar_lambda_irfs.jld2"))
    m = Model1002("ss10", custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    m <= Setting(:impulse_response_horizons, 10)
    dsgevar = DSGEVAR(m, collect(keys(m.exogenous_shocks)), "ss11")
    DSGE.update!(dsgevar, λ = 1.)
    DSGE.update!(dsgevar, jlddata["modal_param"])

    Random.seed!(1793)
    out = impulse_responses(dsgevar, jlddata["data"], :cholesky, 1)
    out_lr = impulse_responses(dsgevar, jlddata["data"], :choleskyLR, 1)
    out_maxbc = impulse_responses(dsgevar, jlddata["data"], :maxBC, 1)

    Random.seed!(1793)
    _ = impulse_responses(dsgevar, jlddata["data"], :cholesky, 1)
    out_lr2 = impulse_responses(dsgevar, jlddata["data"], :cholesky_long_run, 1)
    out_maxbc2 = impulse_responses(dsgevar, jlddata["data"], :maximum_business_cycle_variance,
                                   1)

    Random.seed!(1793)
    out_flip = impulse_responses(dsgevar, jlddata["data"], :cholesky, 1,
                                 flip_shocks = true,)
    out_lr_flip = impulse_responses(dsgevar, jlddata["data"], :choleskyLR, 1,
                                    flip_shocks = true)
    out_maxbc_flip = impulse_responses(dsgevar, jlddata["data"], :maxBC, 1,
                                       flip_shocks = true)

    Random.seed!(1793)
    out_h = impulse_responses(dsgevar, jlddata["data"], :cholesky, 1,
                              horizon = impulse_response_horizons(dsgevar))
    out_lr_h = impulse_responses(dsgevar, jlddata["data"], :choleskyLR, 1,
                              horizon = impulse_response_horizons(dsgevar))
    out_maxbc_h = impulse_responses(dsgevar, jlddata["data"], :maxBC, 1,
                              horizon = impulse_response_horizons(dsgevar))

    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] out
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] out_lr

    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] out_lr2

    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] -out_flip
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] -out_lr_flip

    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] out_h
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] out_lr_h

    # Test maxBC separately b/c these have a slightly different error bound that leads to errors on Julia 1.0 but not Julia 1.1

    if VERSION >= v"1.1"
        @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] out_maxbc
        @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] out_maxbc2
        @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] out_maxbc_h
        @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] -out_maxbc_flip
    else
        @test maximum(abs.(jlddata["exp_modal_maxBC_irf"] - out_maxbc)) < 6e-6
        @test maximum(abs.(jlddata["exp_modal_maxBC_irf"] - out_maxbc2)) < 6e-6
        @test maximum(abs.(jlddata["exp_modal_maxBC_irf"] - out_maxbc_h)) < 6e-6
        @test maximum(abs.(jlddata["exp_modal_maxBC_irf"] + out_maxbc_flip)) < 6e-6
    end
end

@testset "Impulse responses of VAR by using a DSGE as prior and to identify the rotation matrix" begin
    jlddata = load(joinpath(fp, "../../../reference/test_dsgevar_lambda_irfs.jld2"))
    m = Model1002("ss10", custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    m <= Setting(:impulse_response_horizons, 10)
    dsgevar = DSGEVAR(m, collect(keys(m.exogenous_shocks)), "ss11")
    DSGE.update!(dsgevar, λ = 1.)
    DSGE.update!(dsgevar, jlddata["modal_param"])

    Random.seed!(1793)
    out = impulse_responses(dsgevar, jlddata["data"], normalize_rotation = false)
    Random.seed!(1793)
    out_flip = impulse_responses(dsgevar, jlddata["data"]; flip_shocks = true, normalize_rotation = false)
    nobs = size(jlddata["data"], 1)
    lags = DSGE.get_lags(dsgevar)
    k = nobs * lags + 1
    XX = DSGE.lag_data(jlddata["data"], lags; use_intercept = true)
    X̂ = vcat(1, jlddata["data"][:, end], XX[end, 1+1:k - nobs])
    Random.seed!(1793)
    out_X̂ = impulse_responses(dsgevar, jlddata["data"], X̂, normalize_rotation = false)
    Random.seed!(1793)
    out_MM1 = impulse_responses(dsgevar, jlddata["data"]; MM = zeros(DSGE.n_observables(dsgevar),
                                                                     DSGE.n_shocks(dsgevar)), normalize_rotation = false)
    Random.seed!(1793)
    out_MM2 = impulse_responses(dsgevar, jlddata["data"]; MM = rand(DSGE.n_observables(dsgevar),
                                                                    DSGE.n_shocks(dsgevar)), normalize_rotation = false)
    Random.seed!(1793)
    out_draw = impulse_responses(dsgevar, jlddata["data"]; draw_shocks = true, normalize_rotation = false)

    Random.seed!(1793) # now use deviations
    out_dev = impulse_responses(dsgevar, jlddata["data"], deviations = true, normalize_rotation = false)

    Random.seed!(1793) # flip shock
    out_dev_flip = impulse_responses(dsgevar, jlddata["data"], deviations = true,
                                     flip_shocks = true, normalize_rotation = false)

    Random.seed!(1793) # draw shock
    out_dev_draw = impulse_responses(dsgevar, jlddata["data"], deviations = true,
                                     draw_shocks = true, normalize_rotation = false)

    @test @test_matrix_approx_eq jlddata["rotation_irf_by_shock"] out
    @test @test_matrix_approx_eq jlddata["flip_rotation_irf_by_shock"] out_flip
    @test @test_matrix_approx_eq out out_MM1
    @test !(out ≈ out_MM2)
    @test @test_matrix_approx_eq out out_X̂
    @test @test_matrix_approx_eq jlddata["rotation_irf_draw_shock"] out_draw
    @test @test_matrix_approx_eq out_dev jlddata["deviations_rotation_irf_by_shock"]
    @test @test_matrix_approx_eq out_dev_draw jlddata["deviations_rotation_irf_draw_shock"]
    @test @test_matrix_approx_eq out_dev -out_dev_flip
end

@testset "Impulse responses of a VAR approximation to a DSGE (or λ = ∞)" begin
    m = AnSchorfheide()
    m <= Setting(:impulse_response_horizons, 10)
    Random.seed!(1793)
    observables = [:obs_gdp, :obs_nominalrate, :z_t]

    shocks = collect(keys(m.exogenous_shocks))
    fp = dirname(@__FILE__)
    jlddata = load(joinpath(fp, "../../../reference/var_approx_dsge_irfs.jld2"))
    dsgevar = DSGEVAR(m, shocks, "ss0")
    DSGE.update!(dsgevar, lags = 4, observables = observables, λ = 1.)

    out1 = impulse_responses(dsgevar, :cholesky, 1)
    out2 = impulse_responses(dsgevar, :choleskyLR, 1)
    out3 = impulse_responses(dsgevar, :maximum_business_cycle_variance, 1)
    out4 = impulse_responses(dsgevar, :cholesky_long_run, 1)
    out5 = impulse_responses(dsgevar, :maxBC, 1)
    out6 = impulse_responses(dsgevar, :cholesky, 1,
                             flip_shocks = true)
    out7 = impulse_responses(dsgevar, :choleskyLR, 1,
                             flip_shocks = true)
    out8 = impulse_responses(dsgevar, :maxBC, 1,
                             flip_shocks = true)

    out9 = impulse_responses(dsgevar, :cholesky, 1,
                             use_intercept = true)
    out10 = impulse_responses(dsgevar, :choleskyLR, 1,
                              use_intercept = true)
    out11 = impulse_responses(dsgevar, :maxBC, 1,
                              use_intercept = true)
    out12 = impulse_responses(dsgevar, :cholesky, 1,
                              use_intercept = true, flip_shocks = true)
    out13 = impulse_responses(dsgevar, :choleskyLR, 1,
                              use_intercept = true, flip_shocks = true)
    out14 = impulse_responses(dsgevar, :maxBC, 1,
                              use_intercept = true, flip_shocks = true)

    @test @test_matrix_approx_eq jlddata["exp_cholesky"][:, :, 1] out1
    @test @test_matrix_approx_eq jlddata["exp_choleskyLR"][:, :, 1] out2
    @test @test_matrix_approx_eq jlddata["exp_maxBC"][:, :, 1] out3
    @test @test_matrix_approx_eq jlddata["exp_choleskyLR"][:, :, 1] out4
    @test @test_matrix_approx_eq jlddata["exp_maxBC"][:, :, 1] out5
    @test @test_matrix_approx_eq jlddata["exp_cholesky"][:, :, 1] -out6
    @test @test_matrix_approx_eq jlddata["exp_choleskyLR"][:, :, 1] -out7
    @test @test_matrix_approx_eq jlddata["exp_maxBC"][:, :, 1] -out8

    @test @test_matrix_approx_eq jlddata["exp_cholesky_int"][:, :, 1] out9
    @test @test_matrix_approx_eq jlddata["exp_choleskyLR_int"][:, :, 1] out10
    @test @test_matrix_approx_eq jlddata["exp_maxBC_int"][:, :, 1] out11
    @test @test_matrix_approx_eq jlddata["exp_cholesky_int"][:, :, 1] -out12
    @test @test_matrix_approx_eq jlddata["exp_choleskyLR_int"][:, :, 1] -out13
    @test @test_matrix_approx_eq jlddata["exp_maxBC_int"][:, :, 1] -out14
end

nothing


text:
# Implementation of SPN2MILP algorithm
import AlgebraicDecisionDiagrams
# Aliases
const ADD = AlgebraicDecisionDiagrams
const MLExpr = ADD.MultilinearExpression

"""
    spn2milp(spn::SumProductNetwork)

Translates sum-product network `spn` into MAP-equivalent mixed-integer linear program.
Require that sum nodes have exactly two children.
"""
function spn2milp(io, spn::SumProductNetwork, ordering::Union{Nothing,Array{<:Integer}}=nothing)
    # obtain scope of every node
    scopes = scopes(spn)
    # Extract ADDs for each variable
    ## Colect ids of sum nodes
    sumnodes = filter(i -> issum(spn[i]), 1:length(spn))
    ## Create a bucket for each sum node / latent variable
    buckets = Dict{Int,Array{ADD.DecisionDiagram{MLExpr}}}( i => [] for i in sumnodes )
    # variable elimination sequence
    if isnothing(ordering) # if no ordering is given, obtain one
        ## Domain graph
        graph = Dict{Int,Set{Int}}( i => Set{Int}() for i in sumnodes )
    end
    ## First obtain ADDs for manifest variables
    offset = 0 # offset to apply to variable indices at ADD leaves
    vdims = vardims(spn) # var id => no. of values
    pool = ADD.DecisionDiagram{MLExpr}[]
    for var in sort(scopes[1])
        # Extract ADD for variable var
        α = ADD.reduce(extractADD!(Dict{Int,ADD.DecisionDiagram{MLExpr}}(), spn, 1, var, scopes, offset))
        offset += vdims[var] # update start index for next variable
        # Create corresponding optimization variables for leaves (interacts with Gurobi)
        map(t -> begin
            println(io, "binary ", ADD.value(t))
            end,
                Base.filter(n -> isa(n,ADD.Terminal), collect(α))
            )
        push!(pool, α)
        if isnothing(ordering)
            # update domain graph (connect variables in the ADD's scope)
            sc = map(ADD.index, Base.filter(n -> isa(n,ADD.Node), collect(α)))
            for i in sc
                union!(graph[i],sc)
            end
        end
    end
    ## Then build ADDs for sum nodes (latent variables)
    for id in sumnodes
        # construct ADD
        α = ADD.Node(id,MLExpr(spn[id].weights[1]),MLExpr(spn[id].weights[2]))
        # associate ADD to corresponding bucket
        push!(buckets[id],α)
    end
    # Find variable elimination sequence
    if isnothing(ordering)
        ordering = Int[]
        sizehint!(ordering, length(sumnodes))
        while !isempty(graph)
            # find minimum degree node -- break ties by depth/variable id
            deg, k = minimum( p -> (length(p[2]), -p[1]), graph )
            j = -k
            push!(ordering, j)
            # remove j and incident edges
            for k in graph[j]
                setdiff!(graph[k], j)
            end
            delete!(graph, j)
        end
    end
    @assert length(ordering) == length(sumnodes)
    vorder = Dict{Int,Int}() # ordering of elimination of each variable (inverse mapping)
    for i=1:length(ordering)
        vorder[ordering[i]] = i
    end
    # Add ADDs to appropriate buckets
    for α in pool
        # get index of first variable to be eliminated
        i, id = minimum(n -> (vorder[ADD.index(n)],ADD.index(n)), Base.filter(n -> isa(n,ADD.Node), collect(α)))
        push!(buckets[id], α)
    end
    # release pool of ADDs to be collected by garbage collector
    pool = nothing
    # To map each expression in a leaf into a fresh monomial
    cache = Dict{MLExpr,MLExpr}()
    bilinterms = Dict{ADD.Monomial,Int}()
    function renameleaves(e::MLExpr)
        # get!(cache,e,MLExpr(1.0,offset+length(cache)+1))
        # If cached, do nothing
        if haskey(cache, e)
            return cache[e]
        end
        # Generate corresponding variable and constraint (interacts with Gurobi)
        f = MLExpr(1.0,offset+1)
        offset += 1 # increase opt var counter
        # println("continuous ", f)
        idx = [offset] # indices of variables in constraint
        coeff = [-1.0] # coefficients in linear constraint
        terms = String[]
        for (m,c) in e
            if length(m.vars) == 1
                # push!(idx, m.vars[1])
                push!(terms, "$(c)$(m)")
            else # w = x*y
                @assert length(m.vars) == 2
                # Assumes both variables have domain [0,1]
                # Might lead to unfeasible program if this is violated
                if haskey(bilinterms,m)
                    id = bilinterms[m] # id of w
                else # add continuous variable w to problem with 0 ≤ w ≤ 1
                    offset += 1
                    id = offset # id of w
                    bilinterms[m] = id
                end
                push!(terms, "$(c)χ$(id)")
                push!(idx, id)
                # w - x ≤ 0
                println(io, "χ$(id) - χ$(m.vars[1]) <= 0.0")
                # w - y ≤ 0
                println(io, "χ$(id) - χ$(m.vars[2]) <= 0.0")
                # w - y - x ≥ -1
                println(io, "χ$(id) - χ$(m.vars[1]) - χ$(m.vars[1]) >= -1")
            end
            push!(coeff, c)
        end
        # println(idx, coeff)
        # for (i,c) in zip(idx, coeff)
        #     print("$c χ$(id) + ")
        # end
        println(io, "$f = ", join(terms, " + "))
        cache[e] = f
    end
    # Run variable elimination to generate constraints
    α = ADD.Terminal(MLExpr(1.0))
    for i = 1:length(ordering)
        var = ordering[i] # variable to eliminate
        # printstyled("Eliminate: ", var, "\n"; color = :red)
        α = α * reduce(*, buckets[var]; init = ADD.Terminal(MLExpr(1.0)))
        α = ADD.marginalize(α, var)
        # Obtain copy with modified leaves and generate constraints (interacts with JUMP / Gurobi)
        α = ADD.apply(renameleaves, α)
        # For standard bucket elimination (generates tree-decomposition)
        # scope = Base.filter(n -> isa(n,ADD.Node), collect(α))
        # if isempty(scope)
        #     id = ordering[end]
        # else
        #     _, id = minimum(n -> (vorder[ADD.index(n)],ADD.index(n)), scope)
        # end
        # push!(buckets[id], β)
        # printstyled("-> $id\n"; color = :green)
    end
    println(io, "Objective: ", ADD.value(α))
    ADD.value(α)
end

"""
    extractADD!(cache::Dict{Int,ADD.DecisionDiagram{MLExpr}},spn::SumProductNetwork,node::Integer,var::Integer,scopes,offset)

Extract algebraic decision diagram representing the distribution of a variable `var`, using a `cache` of ADDs and `scopes`.
"""
function extractADD!(cache::Dict{Int,ADD.DecisionDiagram{MLExpr}},spn::SumProductNetwork,node::Integer,var::Integer,scopes,offset)
    if haskey(cache, node) return cache[node] end
    if issum(spn[node])
        @assert length(spn[node].children) == 2
        low = extractADD!(cache,spn,spn[node].children[1],var,scopes,offset)
        high = extractADD!(cache,spn,spn[node].children[2],var,scopes,offset)
        γ = ADD.Node(Int(node),low,high)
        cache[node] = γ
        return γ
    elseif isprod(spn[node])
        for j in spn[node].children
            if var in scopes[j]
                γ = extractADD!(cache,spn,j,var,scopes,offset)
                cache[node] = γ
                return γ
            end
        end
        @error "Something went wrong. $var is not in scope of node $(node)."
    elseif isa(spn[node], IndicatorFunction) # leaf
        @assert spn[node].scope == var
        stride = convert(Int, offset + spn[node].value)
        γ = ADD.Terminal(MLExpr(1.0,stride))
        cache[node] = γ
        return γ
    elseif isa(spn[node], CategoricalDistribution) # leaf
        @assert spn[node].scope == var
        expr = mapreduce(i -> MLE(spn[node].values[i], offset + i), +, 1:length(spn[node].values) )
        γ = ADD.Terminal(expr)
        cache[node] = γ
        return γ
    else
        @error "Unsupported node type: $(typeof(spn[node]))."
    end
end


text:
PERFORMANCE = haskey(ENV, "PERFORMANCE")
CODESPEED = get(ENV, "CODESPEED", nothing)
if CODESPEED != nothing && !PERFORMANCE
    error("Cannot submit to Codespeed without enabling performance measurements")
end

if CODESPEED != nothing
    using JSON
    using HTTPClient.HTTPC

    # Setup codespeed data dict for submissions to codespeed's JSON endpoint.  These parameters
    # are constant across all benchmarks, so we'll just let them sit here for now
    pkgdir = dirname(Base.source_path())
    csdata = Dict()
    csdata["commitid"] = chomp(readall(`git -C $pkgdir rev-parse HEAD`))
    csdata["project"] = "CUDA.jl"
    csdata["branch"] = chomp(readall(`git -C $pkgdir symbolic-ref -q --short HEAD`))
    csdata["executable"] = CUDA_VENDOR
    csdata["environment"] = chomp(readall(`hostname`))
    csdata["result_date"] = chomp(readall(`date +'%Y-%m-%d %H:%M:%S'`))
    csdata["revision_date"] = join(split(
            chomp(readall(`git -C $pkgdir log --pretty=format:%cd -n 1 --date=iso`)
        ))[1:2], " " )    # ISO date format minus the timezone
end

# Takes in the raw array of values in vals, along with the benchmark name, description, unit and whether less is better
function submit_to_codespeed(vals,name,desc,unit,test_group,lessisbetter=true)
    csdata["benchmark"] = name
    csdata["description"] = desc
    csdata["result_value"] = mean(vals)
    csdata["std_dev"] = std(vals)
    csdata["min"] = minimum(vals)
    csdata["max"] = maximum(vals)
    csdata["units"] = unit
    csdata["units_title"] = test_group
    csdata["lessisbetter"] = lessisbetter

    println( "$name: $(mean(vals))" )
    ret = post( "http://$CODESPEED/result/add/json/", Dict("json" => json([csdata])) )
    println( json([csdata]) )
    if ret.http_code != 200 && ret.http_code != 202
        error("could not submit $name [HTTP code $(ret.http_code)]")
    end
end

function readable(d)
    if d > 60
        error("unimplemented")
    elseif d < 1
        t = ["m", "µ", "n", "p"]
        for i in 1:length(t)
            scale = 10.0^-3i
            if scale < d <= scale*1000
                return "$(signif(d/scale, 2)) $(t[i])s"
            end
        end
        error()
    else
        return "$(round(d, 2)) s"
    end
end

macro output_timings(t,name,desc,group)
    if CODESPEED == nothing
        ex = quote
            @printf "%-20s: %s ± %s\n" $name readable(mean($t)) readable(std($t))
        end
    else
        ex = quote
            # If we weren't given anything for the test group, infer off of file path!
            test_group = length($group) == 0 ? basename(dirname(Base.source_path())) : $group[1]

            submit_to_codespeed($t, $name, $desc, "seconds", test_group)
        end
    end

    ex
end

const mintrials = 5
const maxtime = 2.5       # in seconds

# TODO: wrap in let -- end?
# TODO: trials shouldn't be accessible afterwards
macro timeit(setup,ex,verification,teardown,name,desc,group...)
    quote
        trials = $mintrials
        t = zeros(trials)

        # warm up and verify
        let
            $(esc(setup))
            e = @elapsed $(esc(ex))
            $(esc(verification))
            $(esc(teardown))
        end
        e = @elapsed ()

        # benchmark
        # TODO: compile-time branch?
        if PERFORMANCE
            i = 1
            start = time()
            while i <= trials
                let
                    $(esc(setup))
                    gc_disable()
                    e = @elapsed $(esc(ex))
                    gc_enable()
                    $(esc(teardown))
                end

                t[i] = e
                if i == trials && (time()-start) < $maxtime
                    # check if accurate enough
                    uncertainty = std(t[1:i])/mean(t[1:i])
                    if uncertainty > .05
                        trials *= 2
                        resize!(t, trials)
                    end
                end

                i += 1
            end
            @output_timings t $name $desc $group
        end
    end
end

# seed rng for more consistent timings
srand(1776)


text:
using OpenCV

println("")
println("This is a simple exmample demonstrating the use of SimpleBlobDetector")
println("")
print("Path to image: ")
img_dir = readline()

img = OpenCV.imread(img_dir)
img_gray = OpenCV.cvtColor(img, OpenCV.COLOR_BGR2GRAY)

OpenCV.namedWindow("Img - Color")
OpenCV.namedWindow("Img - Gray")

OpenCV.imshow("Img - Color", img)
OpenCV.imshow("Img - Gray", img_gray)

OpenCV.waitKey(Int32(0))

OpenCV.destroyAllWindows()

detector = OpenCV.SimpleBlobDetector_create()
kps = OpenCV.detect(detector, img_gray)

println("Number of keypoints: ", size(kps))
for kp in kps
    println(kp.pt, "\t", kp.size)
end


text:
import Base:
    size, length,
    getindex, lastindex, setindex!,
    iterate, similar, *, +, -

"""
    CellField{A<:AbstractArray, G<:Grid} <: Field

A cell-centered field defined on a grid `G` whose values are stored in an `A`.
"""
struct CellField{A<:AbstractArray, G<:Grid} <: Field
    data::A
    grid::G
end

function CellField(::CPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = zeros(T, size(g))
    CellField{typeof(data), typeof(g)}(data, g)
end

function CellField(::GPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    CellField{typeof(data), typeof(g)}(data, g)
end

function CellField(T, ::CPU, g::RegularCartesianGrid)
    data = zeros(T, size(g))
    CellField{typeof(data), typeof(g)}(data, g)
end

function CellField(T, ::GPU, g::RegularCartesianGrid)
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    CellField{typeof(data), typeof(g)}(data, g)
end

"""
    FaceFieldX{A<:AbstractArray, G<:Grid} <: FaceField

An x-face-centered field defined on a grid `G` whose values are stored in an `A`.
"""
struct FaceFieldX{A<:AbstractArray, G<:Grid} <: FaceField
    data::A
    grid::G
end

function FaceFieldX(::CPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = zeros(T, size(g))
    FaceFieldX{typeof(data), typeof(g)}(data, g)
end

function FaceFieldX(::GPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    FaceFieldX{typeof(data), typeof(g)}(data, g)
end

"""
    FaceFieldY{A<:AbstractArray, G<:Grid} <: FaceField

A y-face-centered field defined on a grid `G` whose values are stored in an `A`.
"""
struct FaceFieldY{A<:AbstractArray, G<:Grid} <: FaceField
    data::A
    grid::G
end

function FaceFieldY(::CPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = zeros(T, size(g))
    FaceFieldY{typeof(data), typeof(g)}(data, g)
end

function FaceFieldY(::GPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    FaceFieldY{typeof(data), typeof(g)}(data, g)
end

"""
    FaceFieldZ{T, G<:Grid{T}} <: FaceField{G}

A z-face-centered field defined on a grid `G` whose values are stored in an `A`.
"""
struct FaceFieldZ{A<:AbstractArray, G<:Grid} <: FaceField
    data::A
    grid::G
end

function FaceFieldZ(::CPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = zeros(T, size(g))
    FaceFieldZ{typeof(data), typeof(g)}(data, g)
end

function FaceFieldZ(::GPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    FaceFieldZ{typeof(data), typeof(g)}(data, g)
end

"""
    EdgeField{T<:AbstractArray} <: Field

An edge-centered field defined on a grid `G` whose values are stored in an `A`.
"""
struct EdgeField{A<:AbstractArray, G<:Grid} <: Field
    data::A
    grid::G
end

function EdgeField(::CPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = zeros(T, size(g))
    EdgeField{typeof(data), typeof(g)}(data, g)
end

function EdgeField(::GPU, g::RegularCartesianGrid{T, <:AbstractRange}) where T <: AbstractFloat
    data = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    data .= 0.0
    EdgeField{typeof(data), typeof(g)}(data, g)
end

@inline size(f::Field) = size(f.grid)
@inline length(f::Field) = length(f.data)

@inline getindex(f::Field, inds...) = getindex(f.data, inds...)
# @inline getindex(f::Field, inds...) = f.data[inds...]

@inline lastindex(f::Field) = lastindex(f.data)
@inline lastindex(f::Field, dim) = lastindex(f.data, dim)

@inline setindex!(f::Field, v, inds...) = setindex!(f.data, v, inds...)
# @inline function setindex!(f::Field, v, inds...)
#     f.data[inds...] = v
# end

show(io::IO, f::Field) = show(io, f.data)

iterate(f::Field, state=1) = iterate(f.data, state)
# iterate(f::Field, state=1) = state > length(f) ? nothing : (f.data[state], state+1)

similar(f::CellField{T})  where {T} = CellField(f.metadata, f.grid, f.metadata.float_type)
similar(f::FaceFieldX{T}) where {T} = FaceFieldX(f.metadata, f.grid, f.metadata.float_type)
similar(f::FaceFieldY{T}) where {T} = FaceFieldY(f.metadata, f.grid, f.metadata.float_type)
similar(f::FaceFieldZ{T}) where {T} = FaceFieldZ(f.metadata, f.grid, f.metadata.float_type)
similar(f::EdgeField{T})  where {T} = EdgeField(f.metadata, f.grid, f.metadata.float_type)

set!(u::Field, v) = u.data .= convert(eltype(u.grid), v)
set!(u::Field, v::Field) = @. u.data = v.data

# set!(u::Field{G}, f::Function) where {G<:RegularCartesianGrid} = @. u.data = f(u.grid.xCA, u.grid.yCA, u.grid.zCA)

# Define +, -, and * on fields as element-wise calculations on their data. This
# is only true for fields of the same type, e.g. when adding a FaceFieldY to
# another FaceFieldY, otherwise some interpolation or averaging must be done so
# that the two fields are defined at the same point, so the operation which
# will not be commutative anymore.
for ft in (:CellField, :FaceFieldX, :FaceFieldY, :FaceFieldZ, :EdgeField)
    for op in (:+, :-, :*)
        @eval begin
            # +, -, * a Field by a Number on the left.
            function $op(num::Number, f::$ft)
                ff = similar(f)
                @. ff.data = $op(num, f.data)
                ff
            end

            # +, -, * a Field by a Number on the right.
            $op(f::$ft, num::Number) = $op(num, f)

            # Multiplying two fields together
            function $op(f1::$ft, f2::$ft)
                f3 = similar(f1)
                @. f3.data = $op(f1.data, f2.data)
                f3
            end
        end
    end
end


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

"""
    logvalof(r::NamedTuple{(...,:log,...)})::Real
    logvalof(r::LogDVal)::Real

Extract the log-density value from a result `r`.

Examples:

```julia
logvalof((..., log = log_value, ...)) == log_value
logvalof(LogDVal(log_value)) == log_value
```
"""
function logvalof end
export logvalof

function logvalof(d::Real)
    throw(ArgumentError("Can't get a logarithmic value from $d, unknown if it represents a lin or log value itself."))
end

@inline function logvalof(x::T) where {T<:NamedTuple}
    if hasfield(T, :logval) + hasfield(T, :logd) + hasfield(T, :log) > 1
        throw(ArgumentError("NamedTuples is ambiguous for logvalof contains fields $(join(map(string, filter(name -> name in (:logval, :logd, :log), fieldnames(T))), " and "))"))
    end
    if hasfield(T, :logval)
        x.logval
    elseif hasfield(T, :logd)
        x.logd
    elseif hasfield(T, :log)
        _logvalof_deprecated(x, Val(:log))
    else
        throw(ArgumentError("NamedTuple with fields $(fieldnames(T)) not supported by logvalof, doesn't have a field like :logval"))
    end
end

Base.@noinline function _logvalof_deprecated(x::NamedTuple, ::Val{name}) where name
    Base.depwarn("logvalof support for NamedTuple field $name is deprecated, use NamedTuples with field :logval instead", :logvalof)
    getfield(x, name)
end

"""
    LogDVal{T<:Real}

Represent the logarithm of the value of a statistical density at some point.
`LogDVal` provides means to unambiguously distinguish between linear and
logarithmic values of a density.

Constructor:

    LogDVal(logd::Real)

Use [`logvalof`](@ref) to extract the actual log-density value from
a `LogDVal`:

```julia
    logvalof(LogDVal(d)) == d
```
"""
struct LogDVal{T<:Real}
    logval::T
end

export LogDVal

logvalof(d::LogDVal) = d.logval


text:
mutable struct Entity
    genes::Array{Float64, 1}
    fitness::Float64
end

function generateEntity(genesLength)
   return Entity(rand(-1.0:0.0001:1.0, genesLength), 0)
end

function printEntity(entity1)
   for i in 1:length(entity1.genes)
      print("$(round(entity1.genes[i],digits=4)) ")
   end
   print("\t$(round(entity1.fitness,digits=4))\n")
end

function calculateFitness!(entity1)
   x = entity1.genes[1]
   y = entity1.genes[2]
   z = entity1.genes[3]
   w = entity1.genes[4]
   u = entity1.genes[5]
   v = entity1.genes[6]
   entity1.fitness = abs(4*x^2-6*x-3*y^3+0.5*y+3*z+8*w-6.1*u+2*v-10)
end

function crossover!(entity1, entity2, crossoverPoint)
    for i in 1:crossoverPoint
          x = entity1.genes[i]
          entity1.genes[i] = entity2.genes[i]
          entity2.genes[i] = x
    end
end

function mutate!(entity1, mutationPercentage)
    if rand(Float64)<mutationPercentage
	    mutationPoint = rand(1:length(entity1.genes))
            if entity1.genes[mutationPoint] > 0
   	    	entity1.genes[mutationPoint] = entity1.genes[mutationPoint]-0.1
            else
                entity1.genes[mutationPoint] = entity1.genes[mutationPoint]+0.1
            end
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 100
const MAXN = 500

function gaulegofpi(n)
    aₙ = BigFloat(1)
    bₙ = BigFloat(1) / sqrt(BigFloat(2))
    tₙ = parse(BigFloat, "0.25")
    pₙ = BigFloat(1)

    for i = 1:n
        aₙ₊₁ = (aₙ + bₙ) / BigFloat(2)
        bₙ = sqrt(aₙ * bₙ)
        tₙ = tₙ - pₙ * (aₙ - aₙ₊₁) * (aₙ - aₙ₊₁)
        pₙ = BigFloat(2) * pₙ
        aₙ = aₙ₊₁
    end

    return (aₙ + bₙ) * (aₙ + bₙ) / (BigFloat(4) * tₙ)
end

function main()
    n = GetDigit_n.getn(gaulegofpi, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(gaulegofpi, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(gaulegofpi, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, gaulegofpi(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
module TensorBoardLogger

using ProtoBuf
using CRC32c

#TODO: remove it. Only needed to compute histogram bins.
using StatsBase

# Protobuffer definitions for tensorboard
include("protojl/tensorflow.jl")
include("protojl/summary_pb.jl")
include("protojl/event_pb.jl")

# CRC Utils
include("utils.jl")

# Logging structures
mutable struct Logger
    logdir::String
    file::IOStream
    all_files::Dict{String, IOStream}
    global_step::Int
end

include("logging.jl")
include("Loggers/LogValue.jl")
include("Loggers/LogHistograms.jl")
include("Logger.jl")

#macro tb_log(name)
#    :(_tb_log($(esc(string(name))), $(esc(name))))
#end

export log_histogram, log_value, log_vector
export scalar_summary, histogram_summary, make_event
export Logger

export set_tb_logdir, reset_tb_logs, default_logging_session

end # module


text:
using JuMP
using Gurobi

# compute the gradient
function gradient(gamma1, gamma2, demandsVecCar, demandsVecCar0, demandsVecTruck, demandsVecTruck0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    gradi = zeros(2, numODpairs)
    for i = 1:numODpairs
        gradi[1, i] = 2 * gamma1 * (demandsVecCar[i] - demandsVecCar0[i]) + 2 * gamma2 * sum([(tapFlowVec[1, j] - observFlowVec[1, j]) * jacob[i, j, 1] for j = 1:numLinks])
        gradi[2, i] = 2 * gamma1 * (demandsVecTruck[i] - demandsVecTruck0[i]) + 2 * gamma2 * sum([(tapFlowVec[2, j] - observFlowVec[2, j]) * jacob[i, j, 2] for j = 1:numLinks])
    end
    return gradi
end

# compute a descent direction
function descDirec(gamma1, gamma2, demandsVecCar, demandsVecCar0, demandsVecTruck, demandsVecTruck0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    gradi = gradient(gamma1, gamma2, demandsVecCar, demandsVecCar0, demandsVecTruck, demandsVecTruck0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    h = similar(gradi)
    for k = 1:size(gradi)[1]
        for i = 1:size(gradi)[2]
                h[k, i] = -1 * gradi[k, i]
        end
    end
    return h
end

# compute a search direction
function searchDirec(demandsVecCar, demandsVecTruck, descDirect, epsilon_1)
    h = descDirect
    h_ = similar(h)
    for i = 1:size(h)[2]
        if (demandsVecCar[i] > epsilon_1) || (demandsVecCar[i] <= epsilon_1 && h[1, i] > 0)
            h_[1, i] = h[1, i]
        else
            h_[1, i] = 0
        end
        if (demandsVecTruck[i] > epsilon_1) || (demandsVecTruck[i] <= epsilon_1 && h[2, i] > 0)
            h_[2, i] = h[2, i]
        else
            h_[2, i] = 0
        end
    end
    return h_
end

# line search
function thetaMax(demandsVecCar, demandsVecTruck, searchDirect)
    h_ = searchDirect
    thetaList = Float64[]
    for i = 1:size(h_)[2]
        if h_[1, i] < 0
            push!(thetaList, - demandsVecCar[i]/h_[1, i])
        end
        if h_[2, i] < 0
            push!(thetaList, - demandsVecTruck[i]/h_[2, i])
        end
    end
    theta_max = minimum(thetaList)
    return theta_max
end

# Armijo line search and update
function objF(gamma1, gamma2, demandsVecCar, demandsVecCar0, demandsVecTruck, demandsVecTruck0, fcoeffs)
    demandsDicCar = demandsVecToDic(demandsVecCar)
    demandsDicTruck = demandsVecToDic(demandsVecTruck)
    tapFlowVec = tapMSA_Multi(demandsDicCar, demandsDicTruck, fcoeffs)[2]
    return gamma1 * sum([(demandsVecCar[i] - demandsVecCar0[i])^2 for i = 1:length(demandsVecCar)]) + gamma1 * sum([(demandsVecTruck[i] - demandsVecTruck0[i])^2 for i = 1:length(demandsVecTruck)]) + gamma2 * sum([(tapFlowVec[1, j] - tapFlowVecDict[0][1, j])^2 for j = 1:numLinks]) + gamma2 * sum([(tapFlowVec[2, j] - tapFlowVecDict[0][2, j])^2 for j = 1:numLinks])
end

function armijo(gamma1, gamma2, objFunOld, demandsVecCarOld, demandsVecTruckOld, demandsVecCar0, demandsVecTruck0, fcoeffs, searchDirec, thetaMax, Theta, N)
    demandsVecCarList = Array{Float64}[]
    demandsVecTruckList = Array{Float64}[]
    objFunList = Float64[]
    push!(demandsVecCarList, demandsVecCarOld)
    push!(demandsVecTruckList, demandsVecTruckOld)
    push!(objFunList, objFunOld)
    for n = 0:N
        demandsVecCarNew = similar(demandsVecCarOld)
        demandsVecTruckNew = similar(demandsVecTruckOld)
        for i = 1:length(demandsVecCarOld)
            demandsVecCarNew[i] = demandsVecCarOld[i] + (thetaMax/(Theta^n)) * searchDirec[1, i]
        end
        for i = 1:length(demandsVecTruckOld)
            demandsVecTruckNew[i] = demandsVecTruckOld[i] + (thetaMax/(Theta^n)) * searchDirec[2, i]
        end
    	push!(demandsVecCarList, demandsVecCarNew)
    	push!(demandsVecTruckList, demandsVecTruckNew)
    	push!(objFunList, objF(gamma1, gamma2, demandsVecCarNew, demandsVecCar0, demandsVecTruckNew, demandsVecTruck0, fcoeffs))
    end
    idx = indmin(objFunList)
    objFunNew = objFunList[idx]
    assert(objFunNew <= objFunOld)
    return demandsVecCarList[idx], demandsVecTruckList[idx], objFunNew
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
 _ALGO_HELP_MSG_ = """algorithm (ISRES | PRAXIS | SBPLX | SLSQP | LBFGS | AUGLAG | Newton | Anneal)

  Separate algorithms by commas to run them in succession; e.g. LBGFS,Newton

  Algorithms:
  ISRES   -   Improved stochastic ranking evolution strategy;
              nonlinear rained global optimization;
              note: requires arbitrary bounds on unknowns
  PRAXIS  -   PRincipal AXIS; gradient-free local optimization
  SBPLX   -   Subplex; variant of Nelder-Mead that is more efficient and robust;
              gradient-free local optimization
  SLSQP   -   Sequential quadratic programming for nonlineraly rained,
              gradient-based optimizationa
  LBFGS   -   low-storage BFGS; gradient-based unrained optimization
  AUGLAG  -   Augmented Lagrangian; rained optimization with LBFGS
  Newton  -   My implementation of the Newton-Raphson method
  Anneal  -   My implementation of simulated annealing

  Special syntax: try algorithm1,algorithm2 then backup_algorithms
  This first tries an algorithm that may be less robust but faster and
  if it fails then an algorithm that is more robust and slower.

  Note: in general, I've found SBPLX to be incredibly robust (for this problem).
  If you have no clue at an initial guess, my advice is to start with some
  simulated anneeling, hit with the SBPLX, and then finish it off with a little
  of Newton's method. I've tested this with intentionally pathological initial
  guesses and often converge to the solution anyway :)
""";

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
using Documenter

# include("../src/EquivariantOperators.jl")
# using .EquivariantOperators

include("../src/neural_networks.jl")

##
makedocs(
    sitename = "EquivariantOperators.jl",
    pages = ["index.md", "architecture.md", "publications.md", "tutorials.md"],
)

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
######### Jags line program example  ###########

using StatsPlots, Jags, Statistics

ProjDir = joinpath(dirname(@__FILE__))
cd(ProjDir) do

  line = "
  model {
    for (i in 1:n) {
          mu[i] <- alpha + beta*(x[i] - x.bar);
          y[i]   ~ dnorm(mu[i],tau);
    }
    x.bar   <- mean(x[]);
    alpha    ~ dnorm(0.0,1.0E-4);
    beta     ~ dnorm(0.0,1.0E-4);
    tau      ~ dgamma(1.0E-3,1.0E-3);
    sigma   <- 1.0/sqrt(tau);
  }
  "

  monitors = Dict(
    "alpha" => true,
    "beta" => true,
    "tau" => true,
    "sigma" => true
    )

  jagsmodel = Jagsmodel(
    name="line3",
    model=line,
    monitor=monitors,
    deviance=true, dic=true, popt=true,
    jagsthin=10, thin=1,
    pdir=ProjDir
    );

  println("\nJagsmodel that will be used:")
  jagsmodel |> display

  data = Dict(
    "x" => [1, 2, 3, 4, 5],
    "y" => [1, 3, 3, 3, 5],
    "n" => 5
  )

  inits = [
    Dict("alpha" => 0,"beta" => 0,"tau" => 1),
    Dict("alpha" => 1,"beta" => 2,"tau" => 1),
    Dict("alpha" => 3,"beta" => 3,"tau" => 2),
    Dict("alpha" => 5,"beta" => 2,"tau" => 5)
  ]

  println("Input observed data dictionary:")
  data |> display
  println("\nInput initial values dictionary:")
  inits |> display
  println()

  sim = jags(jagsmodel, data, inits, ProjDir)
  println()

  ## Plotting
  p = plot(sim)

end #cd


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

"""
    BAT.truncate_density(density::AbstractDensity, bounds::AbstractArray{<:Interval})::AbstractDensity

*Experimental feature, not part of stable public API.*

Truncate `density` to `bounds`, the resulting density will be effectively
zero outside of those bounds. In contrast `Distributions.truncated`,
`truncate_density` does *not* renormalize the density.

Requires `varshape(density) isa ArrayShape`.

Only supports densities that are essentially products of univariate
distributions, as well as posterior densities with such densities as priors.
"""
function truncate_density end
export truncate_density

function truncate_density(density::AbstractPosteriorDensity, bounds::AbstractArray{<:Interval})
    @argcheck varshape(density) isa ArrayShape
    PosteriorDensity(getlikelihood(density), truncate_density(getprior(density), bounds))
end

function truncate_density(density::ConstDensity, bounds::AbstractArray{<:Interval})
    old_bounds = var_bounds(density)
    b = HyperRectBounds(bounds)
    new_bounds = ismissing(old_bounds) ? b : var_bounds(density) ∩ b
    ConstDensity(density.value, new_bounds)
end

function truncate_density(density::DistributionDensity{<:MultivariateDistribution}, bounds::AbstractArray{<:Interval})
    r = truncate_dist_hard(density.dist, bounds)
    RenormalizedDensity(DistributionDensity(r.dist), r.logrenormf)
end

function truncate_density(density::RenormalizedDensity{<:DistributionDensity{<:MultivariateDistribution}}, bounds::AbstractArray{<:Interval})
    r = truncate_dist_hard(density.density.dist, bounds)
    RenormalizedDensity(DistributionDensity(r.dist), r.logrenormf + density.logrenormf)
end

"""
    BAT.truncate_dist_hard(dist::Distribution{Univariate}, bounds::Interval)::Distribution{Univariate}
    BAT.truncate_dist_hard(dist::Distribution{Multivariate}, bounds::AbstractArray{<:Interval})::Distribution{Multivariate}

*BAT-internal, not part of stable public API.*

Generalized variant of `Distributions.truncated` - also handles multivariate
distributions and operates on a best-effort basis: If distributions cannot
be truncated, may return the original distribution.

Returns a `NamedTuple`

```julia
    (dist = trunc_dist, logrenormf = logrenormf)
```

with the truncated distribution and the log-PDF amplitude difference to
the original (see [`BAT.trunc_logpdf_ratio`](@ref)).

Mainly used to implement [`BAT.truncate_density`](@ref).
"""
function truncate_dist_hard end

function truncate_dist_hard(dist::Distribution{Univariate}, bounds::AbstractArray{<:Interval})
    @argcheck length(eachindex(bounds)) == 1
    truncate_dist_hard(dist, first(bounds))
end

function truncate_dist_hard(dist::Distribution{Univariate}, bounds::Interval)
    min_lo = minimum(dist)
    max_hi = maximum(dist)
    lo = clamp(minimum(bounds), min_lo, max_hi)
    hi = clamp(max(lo, maximum(bounds)), min_lo, max_hi)

    trunc_dist = truncated(dist, lo, hi)
    logrenormf = trunc_logpdf_ratio(dist, trunc_dist)
    return (dist = trunc_dist, logrenormf = logrenormf)
end

function truncate_dist_hard(dist::Distributions.Truncated, bounds::Interval)
    # Note: Assumes that dist is result of trunctate, not of truncate_dist_hard:
    untrunc_dist = dist.untruncated

    min_lo = minimum(untrunc_dist)
    max_hi = maximum(untrunc_dist)
    lo = clamp(max(minimum(bounds), dist.lower), min_lo, max_hi)
    hi = clamp(max(lo, min(maximum(bounds), dist.upper)), min_lo, max_hi)
    trunc_dist = truncated(untrunc_dist, lo, hi)
    logrenormf = trunc_logpdf_ratio(dist, trunc_dist)
    return (dist = trunc_dist, logrenormf = logrenormf)
end

function truncate_dist_hard(d::Product, bounds::AbstractArray{<:Interval})
    @argcheck length(eachindex(bounds)) == length(d)
    r = truncate_dist_hard.(d.v, bounds)
    trunc_dists = map(x -> x.dist, r)
    logrenormf = sum(x.logrenormf for x in r)

    return (dist = Product(trunc_dists), logrenormf = logrenormf)
end

function truncate_dist_hard(d::StandardMvUniform, bounds::AbstractArray{<:Interval})
    @argcheck length(eachindex(bounds)) == length(d)
    n = length(eachindex(bounds))
    pd = product_distribution(Uniform.(fill(false, n), fill(true, n)))
    return truncate_dist_hard(pd, bounds)
end

function truncate_dist_hard(dist::ConstValueDist, bounds::AbstractVector{<:Interval})
    @argcheck length(eachindex(bounds)) == 0
    (dist = dist, logrenormf = 0)
end

function truncate_dist_hard(dist::NamedTupleDist{names}, bounds::AbstractArray{<:Interval}) where names
    @argcheck length(eachindex(bounds)) == totalndof(varshape(dist))
    distributions = values(dist)
    accessors = values(varshape(dist))

    r = map((dist, acc) -> truncate_dist_hard(dist, view(bounds, ValueShapes.view_idxs(eachindex(bounds), acc))), distributions, accessors)
    trunc_dist = NamedTupleDist(NamedTuple{names}(map(x -> x.dist, r)))
    logrenormf = sum(map(x -> x.logrenormf, r))
    (dist = trunc_dist, logrenormf = logrenormf)
end

function truncate_dist_hard(dist::ValueShapes.UnshapedNTD, bounds::AbstractArray{<:Interval}) where names
    @argcheck length(eachindex(bounds)) == length(dist)
    r = truncate_dist_hard(dist.shaped, bounds)
    (dist = unshaped(r.dist), logrenormf = r.logrenormf)
end

"""
    BAT.trunc_logpdf_ratio(orig_dist::Distribution{TP}, trunc_dist::Distribution{TP})::AbstractFloat

*BAT-internal, not part of stable public API.*

Computes the log-ratio between the amplitude of the PDF of a truncated
distribution and the original (untruncted) distribution, within the support
of the truncated one.

The PDF of both distributions must have the same shape within the support of
`trunc_dist` and may only differ in amplitude.

Mainly used to implement [`BAT.truncate_density`](@ref), in conjunction with
[`BAT.truncate_dist_hard`](@ref). The result contributes to the `logrenormf`
factor of a [`RenormalizedDensity`] that uses truncated distributions internally,
to ensure the density does not get renormalized.
"""
function trunc_logpdf_ratio end

function trunc_logpdf_ratio(orig_dist::Distribution, trunc_dist::Distribution)
    _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist)
end

function _trunc_logpdf_ratio_fallback(orig_dist::Distribution, trunc_dist::Distribution)
    x = rand(bat_determ_rng(), trunc_dist)
    +logpdf(orig_dist, x) - logpdf(trunc_dist, x)
end

function trunc_logpdf_ratio(orig_dist::Distributions.Truncated, trunc_dist::Distributions.Truncated)
    T = promote_type(typeof(trunc_dist.logtp), typeof(orig_dist.logtp))
    if orig_dist.untruncated == trunc_dist.untruncated
        convert(T, trunc_dist.logtp - orig_dist.logtp)
    else
        convert(T, _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist))
    end
end

function trunc_logpdf_ratio(orig_dist::Distribution, trunc_dist::Distributions.Truncated)
    T = typeof(trunc_dist.logtp)
    if orig_dist == trunc_dist.untruncated
        trunc_dist.logtp
    else
        convert(T, _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist))
    end
end


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ClimateMachine.Abstractions
import ClimateMachine.Abstractions: DiscontinuousSpectralElementGrid
using Impero, Printf, MPI, LinearAlgebra, Statistics
include(pwd() * "/sandbox/test_utils.jl")
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
Ω = Circle(0,1) × Circle(0,1)
n = 8
# functional 2D
grid = DiscontinuousSpectralElementGrid(Ω, elements = (2,3), polynomialorder = (n,n), array = ArrayType)
x = view(grid.vgeo, :, grid.x1id, :)   # x-direction
y = view(grid.vgeo, :, grid.x2id, :)   # y-direction
z = view(grid.vgeo, :, grid.x3id, :)   # z-direction

M = view(grid.vgeo, :, grid.Mid, :)    # mass matrix

ω1 = reshape(grid.ω[1], (length(grid.ω[1]), 1, 1))
ω2 = reshape(grid.ω[2], (1,length(grid.ω[2]), 1))
jacobian = M ./ reshape((ω1 .* ω2), (length(grid.ω[1]) * length(grid.ω[2]), 1))

##
using Plots
theme(:default)
gr(size=(300,300))
p1 = scatter(x[:], y[:], label = false, color = :black)
for i in 1:length(grid.vmap⁻)
    if div(i-1,n+1)%4 == 0
        color = :red
    elseif div(i-1,n+1)%4 == 1
        color = :orange
    elseif div(i-1,n+1)%4 == 2
        color = :purple
    else
        color = :green
    end
    scatter!([x[grid.vmap⁻[i]]], [y[grid.vmap⁻[i]]], color = color, label = false)
    display(p1)
    sleep(0.05)
end
##
for i in 1:size(x)[2]
    scatter(x[:], y[:], label = false, color = :black)
    p1 = scatter!(
            x[:,i],
            y[:,i],
            label = false,
            color = :red,
            title = "Element " * string(i))
    display(p1)
    sleep(0.5)
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
function initialize!(integrator, cache::ExplicitRKConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ExplicitRKConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  alg = unwrap_alg(integrator, nothing)
  @unpack A,c,α,αEEst,stages = cache
  @unpack kk = cache

  # Calc First
  kk[1] = integrator.fsalfirst

  # Calc Middle
  for i = 2:stages-1
    utilde = zero(kk[1])
    for j = 1:i-1
      utilde = utilde + A[j,i]*kk[j]
    end
    kk[i] = f(uprev + dt*utilde,p,t+c[i]*dt);
    integrator.destats.nf += 1
  end

  #Calc Last
  utilde = zero(kk[1])
  for j = 1:stages-1
    utilde = utilde + A[j,end]*kk[j]
  end
  kk[end] = f(uprev + dt*utilde,p,t+c[end]*dt)
  integrator.destats.nf += 1
  integrator.fsallast = kk[end] # Uses fsallast as temp even if not fsal

  # Accumulate Result
  utilde = α[1]*kk[1]
  for i = 2:stages
    utilde = utilde + α[i]*kk[i]
  end
  u = uprev + dt*utilde

  if integrator.opts.adaptive
    utilde = (α[1]-αEEst[1]).*kk[1]
    for i = 2:stages
      utilde = utilde + (α[i]-αEEst[i])*kk[i]
    end
    atmp = calculate_residuals(dt * utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if !isfsal(alg.tableau)
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
  end

  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::ExplicitRKCache)
  integrator.kshortsize = 2
  integrator.fsallast = cache.fsallast
  integrator.fsalfirst = cache.kk[1]
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ExplicitRKCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  alg = unwrap_alg(integrator, nothing)
  @unpack A,c,α,αEEst,stages = cache.tab
  @unpack kk,utilde,tmp,atmp = cache

  # Middle
  for i = 2:stages-1
    @.. utilde = zero(kk[1][1])
    for j = 1:i-1
      @.. utilde = utilde + A[j,i]*kk[j]
    end
    @.. tmp = uprev+dt*utilde
    f(kk[i],tmp,p,t+c[i]*dt)
    integrator.destats.nf += 1
  end

  #Last
  @.. utilde = zero(kk[1][1])
  for j = 1:stages-1
    @.. utilde = utilde + A[j,end]*kk[j]
  end
  @.. u = uprev + dt*utilde
  f(kk[end],u,p,t+c[end]*dt) #fsallast is tmp even if not fsal
  integrator.destats.nf += 1

  #Accumulate
  if !isfsal(alg.tableau)
    @.. utilde = α[1]*kk[1]
    for i = 2:stages
      @.. utilde = utilde + α[i]*kk[i]
    end
    @.. u = uprev + dt*utilde
  end

  if integrator.opts.adaptive
    @.. utilde = (α[1]-αEEst[1])*kk[1]
    for i = 2:stages
      @.. utilde = utilde + (α[i]-αEEst[i])*kk[i]
    end
    @.. tmp = dt*utilde
    calculate_residuals!(atmp, tmp, uprev, u,
                         integrator.opts.abstol, integrator.opts.reltol,
                         integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if !isfsal(alg.tableau)
    f(integrator.fsallast,u,p,t+dt)
    integrator.destats.nf += 1
  end
end


text:
@testset "Interpolation" begin
    @testset "T = $(T)" for (T,o,γ) in [(Float64,x -> one(Float64),1.0),(ComplexF64,x -> one(ComplexF64),exp(im*π/4))]
        L = 10.0
        f = x -> γ * sin(2π*x/L)

        N = 100

        rr = range(0, stop=L, length=1000)

        @testset "$(label)" for (label,R,tol) in [
            ("Finite-differences", FiniteDifferences(N, L/(N+1)), 0.01),
            ("Uniform, staggered finite-differences", StaggeredFiniteDifferences(N, L/(N-1/2)), 0.04),
            ("Non-uniform, staggered finite-differences", StaggeredFiniteDifferences(0.1, 0.2, 0.01, L), 0.04),
            ("Implicit finite-differences", ImplicitFiniteDifferences(N, L/(N+1)), 0.01),
            ("FE-DVR", FEDVR(range(0, stop=L, length=N), 4), 5e-6),
            ("B-splines", BSpline(LinearKnotSet(5, 0, L, N)), 5e-9)
        ]
            x = axes(R,1)
            fc = R \ f.(x)

            χ = R[rr,:]
            @test χ*fc ≈ f.(rr) atol=tol
            if !(R isa AbstractFiniteDifferences)
                # We don't bother testing the various
                # finite-differences for interpolations of constants,
                # since they do not yet support non-zero Dirichlet
                # boundary conditions.
                oc = R \ o.(x)
                @test norm(χ*oc .- 1) < 5e-14
            end
        end
    end
end


text:
module PyPlotTools

	import PyCall
	import PyPlot as plt

	const latex_column = 6.52437527778 # inches

	function __init__()
		global mpl_axes_grid = PyCall.pyimport("mpl_toolkits.axes_grid1")
		global rcParams = PyCall.PyDict(plt.matplotlib."rcParams")
		rcParams["image.interpolation"] = "none"
		rcParams["image.origin"] = "lower"
		rcParams["image.resample"] = false
		rcParams["text.latex.preamble"] = "\\usepackage{amsmath}\\usepackage{siunitx}"
		rcParams["savefig.bbox"] = "tight"
		return
	end

	function set_interactive()
		plt.ion()
		rcParams["lines.linewidth"] = 1
		rcParams["lines.markeredgewidth"] = 1
		rcParams["text.usetex"] = false
		rcParams["font.family"] = "sans-serif"
		rcParams["backend"] = "qt5agg"
		rcParams["figure.subplot.left"] = 0.1
		rcParams["figure.subplot.right"] = 0.9
		rcParams["figure.subplot.bottom"] = 0.1
		rcParams["figure.subplot.top"] = 0.9
		return
	end

	function set_eps()
		plt.ioff()
		rcParams["lines.linewidth"] = 2
		rcParams["lines.markeredgewidth"] = 1.25
		rcParams["text.usetex"] = true
		rcParams["font.family"] = "serif"
		rcParams["font.serif"] = "Computer Modern"
		plt.matplotlib.cm.viridis.set_bad("w", 1)
		rcParams["backend"] = "ps"
		rcParams["figure.subplot.left"] = 0
		rcParams["figure.subplot.right"] = 1
		rcParams["figure.subplot.bottom"] = 0
		rcParams["figure.subplot.top"] = 1
		return
	end

	function add_colourbar(fig, ax, image; size="5%", pad=0.05, kwargs...)
		divider = mpl_axes_grid.make_axes_locatable(ax)
		horizontal = false
		if haskey(kwargs, :orientation) && kwargs[:orientation] == "horizontal"
			cax = divider.new_vertical(;size, pad, pack_start=false) # Also: .append
			horizontal = true
		else
			cax = divider.new_horizontal(;size, pad, pack_start=false)
		end
		fig.add_axes(cax)
		colourbar = plt.colorbar(image, cax=cax; kwargs...)
		if horizontal
			# Set ticks to top since it makes more sense to look at the colourbar
			# first because looking at the image (What am I looking at?)
			cax.xaxis.set_tick_params(top=true, bottom=false, labeltop=true, labelbottom=false)
			cax.xaxis.set_label_position("top")
		end
		return cax, colourbar
	end

	function add_roi(ax, roi::NTuple{2, AbstractRange}; linewidth=1, edgecolor="r", facecolor="none")
		ax.add_patch(
			plt.matplotlib[:patches][:Rectangle](
				(roi[2].start - 1, roi[1].start - 1),
				roi[2].stop - roi[2].start,
				roi[1].stop - roi[1].start;
				transform=ax.transData,
				linewidth, edgecolor, facecolor
			)
		)
	end

end


text:
#=
    This file includes all the functions needed to check the power balance
    constraints
=#

function compute_power_balances!(h::Int64, y::Int64, s::Int64, mg::Microgrid)
    # Hydrogen
    checking!(h, y, s, mg, Hydrogen)
    # Heat
    checking!(h, y, s, mg, Heat)
    # Electricity
    checking!(h, y, s, mg, Electricity)
end

function power_balance(h::Union{Int64, UnitRange{Int64}}, y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, type::DataType)
    # Energy balance
    balance = 0.
    # Demands
    if !isempty(mg.demands)
        for a in mg.demands
            a.carrier isa type ? balance = balance .+ a.carrier.power[h,y,s] : nothing
        end
    end
    # Generations
    if !isempty(mg.generations)
        for a in mg.generations
            a.carrier isa type ? balance = balance .- a.carrier.power[h,y,s] : nothing
        end
    end
    # Storages
    if !isempty(mg.storages)
        for a in mg.storages
            a.carrier isa type ? balance = balance .- a.carrier.power[h,y,s] : nothing
        end
    end
    # Converters
    if !isempty(mg.converters)
        for a in mg.converters
            for c in a.carrier
                c isa type ? balance = balance .- c.power[h,y,s] : nothing
            end
        end
    end
    return balance
end

function checking!(h::Int64, y::Int64, s::Int64, mg::Microgrid, type::DataType)
    # Energy balance
    balance = power_balance(h, y, s, mg, type)
    # Grids
    if !isempty(mg.grids)
        for a in mg.grids
            if a.carrier isa type
                a.carrier.power[h,y,s] = min(a.powerMax, max(-a.powerMax, balance))
            end
        end
    else # If the balance equation is not fulfilled, systems are turned to zerp
        if balance > sum(a.carrier.power[h,y,s] for a in mg.demands if a.carrier isa type) + 1e-2 # 0.01 kW tolerance
            # Storage set to zero
            for a in mg.storages
                if a.carrier isa type
                    a.carrier.power[h,y,s]= 0.
                    a.soc[h+1,y,s] = max(0., a.soc[h,y,s] * (1. - a.η_self))
                end
            end
            # Converters set to zero
            for a in mg.converters
                for c in a.carrier
                    c.power[h,y,s] = 0.
                end
            end
        else
            println("Shedding energy demand!")
        end
    end
end


text:
using MechanicalSketch
import MechanicalSketch: O, label, background, sethue, setdash, fontface, fontsize, line
import MechanicalSketch: FS, EM, translate, color_with_lumin, midpoint, getrotation, rotate
import MechanicalSketch: empty_figure, finish, color_from_palette

let
    this_fig = empty_figure(filename = joinpath(@__DIR__, "test_4.png"))
    background(color_from_palette("seagreen1"))
    sethue(color_with_lumin(color_from_palette("darkblue"), 20))

    label("O", :NW, O)
    include("test_functions_4.jl")
    labelelledtri()
    rotate(π / 3)
    labelelledtri()
    MechanicalSketch.scale(0.5)
    rotate(2π / 3)
    labelelledtri()
    translate(O + (10EM, 0))
    labelelledtri()
    rotate(π)
    labelelledtri()
    label("-----------")
    finish()
end


text:
export
    canonise!,
    truncate!,
    compress!,
    compress

# This is for backwards compatibility
function compress(ϕ::AbstractMPS, Dcut::Int, tol::Number=1E-8, max_sweeps::Int=4, args...)
    ψ = copy(ϕ)
    compress!(ψ, Dcut, tol, max_sweeps, args...)
    ψ
end

function compress!(ϕ::AbstractMPS, Dcut::Int, tol::Number=1E-8, max_sweeps::Int=4, args...)
    # right canonise ϕ
    _left_sweep!(ϕ, args...)

    # Initial guess - truncated ϕ
    ψ = copy(ϕ)
    _right_sweep!(ϕ, Dcut, args...)

    # Create environment
    env = left_env(ϕ, ψ)

    # Variational compression
    overlap = Inf
    overlap_before = -Inf

    @info "Compressing state down to" Dcut

    for sweep ∈ 1:max_sweeps
        _left_sweep_var!!(ϕ, env, ψ, args...)
        overlap = _right_sweep_var!!(ϕ, env, ψ, args...)

        diff = abs(overlap_before - abs(overlap))
        @info "Convergence" diff

        if diff < tol
            @info "Finished in $sweep sweeps of $(max_sweeps)."
            return overlap
        else
            overlap_before = overlap
        end
    end
    overlap
end

function truncate!(ψ::AbstractMPS, s::Symbol, Dcut::Int=typemax(Int), args...)
    @assert s ∈ (:left, :right)
    if s == :right
        _right_sweep!(ψ, args...)
        _left_sweep!(ψ, Dcut, args...)
    else
        _left_sweep!(ψ, args...)
        _right_sweep!(ψ, Dcut, args...)
    end
end

canonise!(ψ::AbstractMPS, s::Symbol) = canonise!(ψ, Val(s))
canonise!(ψ::AbstractMPS, ::Val{:right}) = _left_sweep!(ψ, typemax(Int))
canonise!(ψ::AbstractMPS, ::Val{:left}) = _right_sweep!(ψ, typemax(Int))

function _right_sweep!(ψ::AbstractMPS, Dcut::Int=typemax(Int), args...)
    R = ones(eltype(ψ), 1, 1)
    for (i, A) ∈ enumerate(ψ)
        @matmul M̃[(x, σ), y] := sum(α) R[x, α] * A[α, σ, y]
        Q, R = qr_fact(M̃, Dcut, args...)
        R = R ./ maximum(abs.(R))
        @cast A[x, σ, y] := Q[(x, σ), y] (σ ∈ 1:size(A, 2))
        ψ[i] = A
    end
end

function _left_sweep!(ψ::AbstractMPS, Dcut::Int=typemax(Int), args...)
    R = ones(eltype(ψ), 1, 1)
    for i ∈ length(ψ):-1:1
        B = ψ[i]
        @matmul M̃[x, (σ, y)] := sum(α) B[x, σ, α] * R[α, y]
        R, Q = rq_fact(M̃, Dcut, args...)
        R = R ./ maximum(abs.(R))
        @cast B[x, σ, y] := Q[x, (σ, y)] (σ ∈ 1:size(B, 2))
        ψ[i] = B
    end
end

function _left_sweep_var!!(ϕ::AbstractMPS, env::Vector{<:AbstractMatrix}, ψ::AbstractMPS, args...)
    env[end] = ones(eltype(ϕ), 1, 1)

    for i ∈ length(ψ):-1:1
        L, R = env[i], env[i+1]

        # optimize site
        M = ψ[i]
        @tensor MM[x, σ, α] := L[x, β] * M[β, σ, α]
        @matmul MM[x, (σ, y)] := sum(α) MM[x, σ, α] * R[α, y]

        _, Q = rq_fact(MM, args...)
        @cast B[x, σ, y] := Q[x, (σ, y)] (σ ∈ 1:size(M, 2))

        # update ϕ and right environment
        ϕ[i] = B
        A = ψ[i]

        @tensor RR[x, y] := A[x, σ, α] * R[α, β] * conj(B)[y, σ, β] order = (β, α, σ)
        env[i] = RR
    end
    env[1][1]
end

function _right_sweep_var!!(ϕ::AbstractMPS, env::Vector{<:AbstractMatrix}, ψ::AbstractMPS, args...)
    env[1] = ones(eltype(ϕ), 1, 1)

    for (i, M) ∈ enumerate(ψ)
        L, R = env[i], env[i+1]

        # optimize site
        @tensor M̃[x, σ, α] := L[x, β] * M[β, σ, α]
        @matmul B[(x, σ), y] := sum(α) M̃[x, σ, α] * R[α, y]

        Q, _ = qr_fact(B, args...)
        @cast A[x, σ, y] := Q[(x, σ), y] (σ ∈ 1:size(M, 2))

        # update ϕ and left environment
        ϕ[i] = A
        B = ψ[i]

        @tensor LL[x, y] := conj(A[β, σ, x]) * L[β, α] * B[α, σ, y] order = (α, β, σ)
        env[i+1] = LL
    end
    env[end][1]
end

function _right_sweep(A::AbstractArray, Dcut::Int=typemax(Int), args...) where {T <: AbstractMPS}
    rank = ndims(A)
    ψ = MPS(eltype(A), rank)
    R = reshape(copy(A), (1, length(A)))

    for i ∈ 1:rank
        d = size(A, i)
        @cast M[(x, σ), y] := R[x, (σ, y)] (σ ∈ 1:d)
        Q, R = qr_fact(M, Dcut, args...)
        R = R ./ maximum(abs.(R))
        @cast B[x, σ, y] := Q[(x, σ), y] (σ ∈ 1:d)
        ψ[i] = B
    end
    ψ
end

function _left_sweep(A::AbstractArray, Dcut::Int=typemax(Int), args...) where {T <: AbstractMPS}
    rank = ndims(A)
    ψ = MPS(eltype(A), rank)
    R = reshape(copy(A), (length(A), 1))

    for i ∈ rank:-1:1
        d = size(A, i)
        @cast M[x, (σ, y)] := R[(x, σ), y] (σ ∈ 1:d)
        R, Q = rq_fact(M, Dcut, args...)
        R = R ./ maximum(abs.(R))
        @cast B[x, σ, y] := Q[x, (σ, y)] (σ ∈ 1:d)
        ψ[i] = B
    end
    ψ
end