text:
struct Config
    user::String                    # the OS username of the user running the server
    nodes::Dict{Type,Vector{Int}}   # the pids for the nodes on the cluster
    cpus::Dict{Int,Vector{Int}}     # the indices of the cpus per node
    auth::GitHub.Authorization      # the GitHub authorization used to post statuses/reports
    secret::String                  # the GitHub secret used to validate webhooks
    trackrepo::String               # the main Julia repo tracked by the server
    reportrepo::String              # the repo to which result reports are posted
    trigger::Regex                  # a regular expression to match comments against
    workdir::String                 # the server's work directory
    admin::String                   # GitHub handle of the server administrator
    bucket::Union{Nothing,String}   # AWS bucket to upload large files too
    testmode::Bool                  # if true, jobs will run as test jobs

    function Config(user, nodes, auth, secret;
                    workdir = pwd(),
                    cpus = Dict{Int,Vector{Int}}(),
                    trackrepo = "JuliaLang/julia",
                    reportrepo = "JuliaCI/NanosoldierReports",
                    trigger =  r"\@nanosoldier\s*`runbenchmarks\(.*?\)`",
                    admin = "",
                    bucket = nothing,
                    testmode = false)
        isempty(nodes) && throw(ArgumentError("need at least one node to work on"))
        return new(user, nodes, cpus, auth, secret, trackrepo,
                   reportrepo, trigger, workdir, admin, bucket, testmode)
    end
end

# the shared space in which child nodes can work
workdir(config::Config) = config.workdir

# the report repository
reportrepo(config::Config) = config.reportrepo

# the local directory of the report repository
reportdir(config::Config) = joinpath(workdir(config), split(reportrepo(config), "/")[2])

persistdir!(path) = (isdir(path) || mkdir(path); return path)

function persistdir!(config::Config)
    persistdir!(workdir(config))
    if isdir(reportdir(config))
        gitreset!(reportdir(config))
    else
        gitclone!(reportrepo(config), reportdir(config), config.auth)
    end
end

function nodelog(config::Config, node, message; error=nothing)
    time = now()
    if error !== nothing
        @error "[Node $node | $time]: Encountered error: $message" exception=error
    else
        @info "[Node $node | $time]: $message"
    end
    persistdir!(workdir(config))
    path = joinpath(workdir(config), "node$(node).log")
    open(path, "a") do file
        chmod(path, 0o660)
        println(file, time, " | ", node, " | ", message)
    end
end

# the list of CPUs for a given node
mycpus(config::Config, node=getpid()) = get(config.cpus, node, 1:(Sys.CPU_THREADS-1))


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Utilities used in the solver
#--------------------------------------------------------------------

"""
    Takes a set of tuples, and flattens them into
    a vector.
"""
function reshapeFromTuple(U::NTuple{N, Field}) where {N}
    return vcat(reshape.(U)...)
end

""" Takes a vector and reshapes them into tuples. Inverse of
    reshapeFromTuple.
"""
function reshapeToTuple(space::S, N::Int, x::Array{T,1})::NTuple{N, Field}  where {S, T}
    U = reshape(x, :, N)
    return Tuple(reshape(space, U[:, i]) for i in 1:N)
end

function enforcebc!(F::NTuple{N, Field{S}}, B::NTuple{N, Field{S}}) where {S,N}
    r = Field(first(F).space, (u,v)->v-u)
    for index_ in 1:N
        F[index_].value[end, :] = B[index_].value[end, :]
        F[index_].value[:, end] = B[index_].value[:, end]
    end
end

function enforceregularity!(F::NTuple{N, Field{S}}, A::NTuple{N, Field{S}}) where {S,N}
    r = Field(first(F).space, (u,v)->v-u)
    for index in CartesianIndices(r.value)
        if r.value[index] == 0.0
            for index_ in 1:N
                F[index_].value[index] = A[index_].value[index]
            end
        end
    end
end

ext: jl


text:
push!(LOAD_PATH, ".")
using WildBootTests
using StatFiles, StatsModels, DataFrames, DataFramesMeta, BenchmarkTools, Plots, CategoricalArrays, Random, StableRNGs

open("unittests.log", "w") do log  # use Github (Desktop) to detect changes in output

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\collapsed.dta"))
dropmissing!(df)
f = @formula(hasinsurance ~ 1 + selfemployed + post + post_self)
f = apply_schema(f, schema(f, df, Dict(:hasinsurance => ContinuousTerm)))
resp, predexog = modelcols(f, df)

println(log, "\nboottest post_self=.04, weight(webb)")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), auxwttype=WildBootTests.webb, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nboottest post_self=.04, weight(webb) reps(9999999) noci")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999999, auxwttype=WildBootTests.webb, getCI=false, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test))")

println(log, "\nregress hasinsurance selfemployed post post_self, cluster(year)")
println(log, "boottest (post_self=.05) (post=-.02), reps(9999) weight(webb)")
test = wildboottest([0 0 0 1; 0 0 1 0], [0.05; -0.02]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=WildBootTests.webb, rng=StableRNG(1231))
println(log, "F=$(teststat(test)) p=$(p(test))")

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))[:,[:wage; :tenure; :ttl_exp; :collgrad; :industry]]
dropmissing!(df)
desc = describe(df, :eltype)
for i in axes(desc, 1)  # needed only for lm()
  if desc[i,:eltype] == Float32
    sym = desc[i,:variable]
    @transform!(df, @byrow $sym=Float64($sym))
  end
end
f = @formula(wage ~ 1 + tenure + ttl_exp + collgrad)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

println(log, "\nconstraint 1 ttl_exp = .2")
println(log, "cnsreg wage tenure ttl_exp collgrad, constr(1) cluster(industry)")
println(log, "boottest tenure")
test = wildboottest([0 1 0 0], [0]; R1=[0 0 1 0], r1=[.2], resp, predexog, clustid=df.industry, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union), cluster(industry)")
println(log, "boottest tenure, ptype(equaltail)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage; :tenure; :ttl_exp; :collgrad; :industry; :union]]
dropmissing!(df)
f = @formula(wage ~ 1 + ttl_exp + collgrad)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ union)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=9999, ptype=WildBootTests.equaltail, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nboottest tenure, ptype(equaltail) reps(9999) weight(webb) stat(c)")
test = wildboottest([0 0 0 1], [.0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=9999, auxwttype=WildBootTests.webb, bootstrapc=true, ptype=WildBootTests.equaltail, rng=StableRNG(1231), gridmin=[-5], gridmax=[5], gridpoints=[100])
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
plot(plotpoints(test)...)

println(log, "\nboottest, ar")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, ARubin=true, reps=9999, rng=StableRNG(1231))
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
plot(plotpoints(test)...)

println(log, "\nboottest, ar nonull")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, ARubin=true, reps=9999, rng=StableRNG(1231), imposenull=false)
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nscoretest tenure")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=0)
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nwaldtest tenure")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=0, imposenull=false)
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nivregress liml wage (tenure = collgrad ttl_exp), cluster(industry)")
println(log, "boottest tenure")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage, :tenure, :ttl_exp, :collgrad, :industry]]
dropmissing!(df)
f = @formula(wage ~ 1)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ collgrad + ttl_exp)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 1], [0]; resp, predexog, predendog, inst, LIML=true, clustid=df.industry, small=false, reps=999, rng=StableRNG(1231))
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nivreg2 wage collgrad smsa race age (tenure = union married), cluster(industry) fuller(1)")
println(log, "boottest tenure, nograph weight(webb) reps(9999)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage, :tenure, :ttl_exp, :collgrad, :smsa, :race, :age, :union, :married, :industry]]
dropmissing!(df)
f = @formula(wage ~ 1 + collgrad + smsa + race + age)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ union + married)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 0 0 0 0 1], [0]; resp, predexog, predendog, inst, Fuller=1, clustid=df.industry, small=false, reps=9999, auxwttype=WildBootTests.webb, rng=StableRNG(1231))
println(log, "z=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nareg wage ttl_exp collgrad tenure [aw=hours] if occupation<., cluster(age) absorb(industry)")
println(log, "boottest tenure, cluster(age occupation) bootcluster(occupation)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage, :ttl_exp, :collgrad, :tenure, :age, :industry, :occupation, :hours]]
dropmissing!(df)
f = @formula(wage ~ ttl_exp + collgrad + tenure)  # constant unneeded in FE model
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
test = wildboottest([0 0 1], [0]; resp, predexog, clustid=Matrix(df[:, [:occupation, :age]]), nbootclustvar=1, nerrclustvar=2, obswt=df.hours, feid=df.industry, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\nglobal pix lnkm pixpetro pixdia pixwaterd pixcapdist pixmal pixsead pixsuit pixelev pixbdist")
println(log, "global geo lnwaterkm lnkm2split mean_elev mean_suit malariasuit petroleum diamondd")
println(log, "global poly capdistance1 seadist1 borderdist1")
println(log, "encode pixwbcode, gen(ccode)  // make numerical country identifier")
println(log, "qui areg lnl0708s centr_tribe lnpd0 \$pix \$geo \$poly, absorb(ccode)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(ccode)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(pixcluster)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(ccode pixcluster)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Work\Econometrics\Wild cluster\pixel-level-baseline-final.dta"))
pix  = [:lnkm, :pixpetro, :pixdia, :pixwaterd, :pixcapdist, :pixmal, :pixsead, :pixsuit, :pixelev, :pixbdist]
geo  = [:lnwaterkm, :lnkm2split, :mean_elev, :mean_suit, :malariasuit, :petroleum, :diamondd]
poly = [:capdistance1, :seadist1, :borderdist1]
df = df[:,[pix; geo; poly; :lnl0708s; :centr_tribe; :lnpd0; :pixwbcode; :pixcluster]]
dropmissing!(df)
df.ccode = levelcode.(categorical(df.pixwbcode, compress=true))
df.pixcode = levelcode.(categorical(df.pixcluster, compress=true))
f = Term(:lnl0708s) ~ sum(term.([:centr_tribe; :lnpd0; pix; geo; poly]))
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:ccode, :pixcode]]), nbootclustvar=1, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:pixcode, :ccode]]), nbootclustvar=1, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:pixcode, :ccode]]), nbootclustvar=2, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

println(log, "\ninfile coll merit male black asian year state chst using regm.raw, clear")
println(log, "qui regress coll merit male black asian i.year i.state if !inlist(state,34,57,59,61,64,71,72,85,88), cluster(state)	")
println(log, "generate individual = _n  // unique ID for each observation")
println(log, "boottest merit, nogr reps(9999) gridpoints(10)  // defaults to bootcluster(state)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) bootcluster(state year)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull bootcluster(state year)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) bootcluster(individual)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull bootcluster(individual)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Work\Econometrics\Wild cluster\regm.dta"))
df = DataFrame(coll=Bool.(df.coll), merit=Bool.(df.merit), male=Bool.(df.male), black=Bool.(df.black), asian=Bool.(df.asian), state=categorical(Int8.(df.state)), year=categorical(Int16.(df.year)))
dropmissing!(df)
df = df[df.state .∉ Ref([34,57,59,61,64,71,72,85,88]),:]
f = @formula(coll ~ 1 + merit + male + black + asian + year + state)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=levelcode.(df.state), gridpoints=[10], reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=levelcode.(df.state), reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=[levelcode.(df.year) levelcode.(df.state)], nbootclustvar=2, nerrclustvar=1, reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=[levelcode.(df.year) levelcode.(df.state)], nbootclustvar=2, nerrclustvar=1, reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid= [collect(1:nrow(df)) levelcode.(df.state)], nbootclustvar=1, nerrclustvar=1, reps=9999, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid= [collect(1:nrow(df)) levelcode.(df.state)], nbootclustvar=1, nerrclustvar=1, reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, "t=$(teststat(test)) p=$(p(test)) CI=$(CI(test))")

# d=download("https://raw.github.com/droodman/WildBootTests.jl/master/data/Levitt.dta", tempname() * ".dta")
# df = DataFrame(load(d))
# df = df[:, [:lViolentpop,: lPropertypop, ]]

# df = DataFrame(load(raw"D:\OneDrive\Documents\Macros\WildBootTests.jl"))
# df = df[:, [:proposition_vote, :treatment, :ideology1, :log_income, :Q1_immigration, :group_id1, :group_id2]]
# dropmissing!(df)
# f = @formula(proposition_vote ~ treatment + ideology1 + log_income + Q1_immigration + 1)
# f = apply_schema(f, schema(f, df, Dict(:Q1_immigration  => CategoricalTerm)))
# resp, predexog = modelcols(f, df)
# test = WildBootTests.wildboottest(([0 0 0 0 0 0 0 0 0 0 0 0 1], [0]); resp, predexog, clustid=Matrix(df[:, [:Q1_immigration, :group_id1, :group_id2]]), rng=StableRNG(1231))

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test
using Base: catch_stack

@testset "Basic exception stack handling" begin
    # Exiting the catch block normally pops the exception
    try
        error("A")
    catch
        @test length(catch_stack()) == 1
    end
    @test length(catch_stack()) == 0
    # Exiting via a finally block does not pop the exception
    try
        try
            error("A")
        finally
            @test length(catch_stack()) == 1
        end
    catch
        @test length(catch_stack()) == 1
    end
    # The combined try-catch-finally form obeys the same rules as above
    try
        error("A")
    catch
        @test length(catch_stack()) == 1
    finally
        @test length(catch_stack()) == 0
    end
    @test length(catch_stack()) == 0
    # Errors are pushed onto the stack according to catch block nesting
    try
        error("RootCause")
    catch
        @test length(catch_stack()) == 1
        try
            error("B")
        catch
            stack = catch_stack()
            @test length(stack) == 2
            @test stack[1][1].msg == "RootCause"
            @test stack[2][1].msg == "B"
        end
        # Stack pops correctly
        stack = catch_stack()
        @test length(stack) == 1
        @test stack[1][1].msg == "RootCause"
    end
end

@testset "Exception stack lowering special cases" begin
    # try block in value position
    val = try
        error("A")
    catch
        @test length(catch_stack()) == 1
        1
    end
    @test val == 1
    function test_exc_stack_tailpos()
        # try block in tail position
        try
            error("A")
        catch
            length(catch_stack())
        end
    end
    @test test_exc_stack_tailpos() == 1
    @test length(catch_stack()) == 0
end

@testset "Exception stacks - early exit from try or catch" begin
    # Exiting a catch block early with normal control flow — break, continue,
    # return, goto — will result in popping of the exception stack.
    function test_exc_stack_catch_return()
        try
            error("A")
        catch
            @test length(catch_stack()) == 1
            return
        end
    end
    test_exc_stack_catch_return()

    for i=1:1
        try
            error("A")
        catch
            @test length(catch_stack()) == 1
            break
        end
    end
    # Also test try-break-finally forms here. See #31766
    for i=1:1
        try
            error("A")
        catch
            @test length(catch_stack()) == 1
            break
        finally
            @test length(catch_stack()) == 0
        end
    end
    @test length(catch_stack()) == 0

    for i=1:1
        try
            error("A")
        catch
            @test length(catch_stack()) == 1
            continue
        end
    end
    for i=1:1
        try
            error("A")
        catch
            @test length(catch_stack()) == 1
            continue
        finally
            @test length(catch_stack()) == 0
        end
    end
    @test length(catch_stack()) == 0

    try
        error("A")
    catch
        @test length(catch_stack()) == 1
        @goto outofcatch
    end
    @label outofcatch
    try
        error("A")
    catch
        @test length(catch_stack()) == 1
        @goto outofcatch2
    finally
        @test length(catch_stack()) == 0
    end
    @label outofcatch2
    @test length(catch_stack()) == 0

    # Exiting from a try block in various ways should not affect the exception
    # stack state.
    try
        error("ExceptionInOuterTry")
    catch
        @test length(catch_stack()) == 1
        function test_exc_stack_try_return()
            try
                return
            catch
            end
        end
        test_exc_stack_try_return()
        for i=1:1
            try
                break
            catch
            end
        end
        for i=1:1
            try
                continue
            catch
            end
        end
        try
            @goto outoftry
        catch
        end
        @label outoftry
        @test length(catch_stack()) == 1
        @test catch_stack()[1][1] == ErrorException("ExceptionInOuterTry")
    end
end

@testset "Finally handling with exception stacks" begin
    # The lowering of finally is quite subtle when combined with break or
    # return because each finally block may be entered via multiple code paths
    # (eg, different occurrences of return), and these code paths must diverge
    # again once the finally block has completed. To complicate matters
    # further, the return code path must thread through every nested finally
    # block before actually returning, all the while preserving the information
    # about which variable to return.

    # Issue #34579
    (()-> begin
        try
            throw("err")
        catch
            # Explicit return => exception should be popped before finally block
            return
        finally
            @test length(Base.catch_stack()) == 0
        end
    end)()
    @test length(Base.catch_stack()) == 0

    while true
        try
            error("err1")
        catch
            try
                # Break target is outside catch block, but finally is inside =>
                # exception should not be popped inside finally block
                break
            finally
                @test length(Base.catch_stack()) == 1
            end
        end
    end
    @test length(Base.catch_stack()) == 0

    # Nested finally handling with `return`: each finally block should observe
    # only the active exceptions as according to its nesting depth.
    (() -> begin
        try
            try
                error("err1")
            catch
                try
                    try
                        error("err2")
                    catch
                        # This return needs to thread control flow through
                        # multiple finally blocks in the linearized IR.
                        return
                    end
                finally
                    # At this point err2 is dealt with
                    @test length(Base.catch_stack()) == 1
                    @test Base.catch_stack()[1][1] == ErrorException("err1")
                end
            end
        finally
            # At this point err1 is dealt with
            @test length(Base.catch_stack()) == 0
        end
    end)()
    @test length(Base.catch_stack()) == 0
end

@testset "Deep exception stacks" begin
    # Generate deep exception stack with recursive handlers.
    #
    # (Note that if you let this overflow the program stack (not the exception
    # stack) julia will crash. See #28577.)
    function test_exc_stack_deep(n)
        n != 1 || error("RootCause")
        try
            test_exc_stack_deep(n-1)
        catch
            error("n==$n")
        end
    end
    @test try
        test_exc_stack_deep(100)
    catch
        @test catch_stack()[1][1] == ErrorException("RootCause")
        length(catch_stack())
    end == 100
    @test length(catch_stack()) == 0
end

@testset "Exception stacks and Tasks" begin
    # Task switching should not affect exception state. See #12485.
    try
        error("A")
    catch
        t = @task try
            error("B")
        catch exc
            exc
        end
        yield(t)
        @test t.state == :done
        @test t.result == ErrorException("B")
        # Task exception state is preserved around task switches
        @test length(catch_stack()) == 1
        @test catch_stack()[1][1] == ErrorException("A")
    end
    @test length(catch_stack()) == 0
    # test rethrow() rethrows correct state
    bt = []
    try
        try
            error("A")
        catch
            bt = catch_backtrace()
            t = @task try
                error("B")
            catch exc
                exc
            end
            yield(t)
            @test t.state == :done
            @test t.result == ErrorException("B")
            @test bt == catch_backtrace()
            rethrow()
        end
    catch exc
        @test exc == ErrorException("A")
        @test bt == catch_backtrace()
    end
    @test length(catch_stack()) == 0
    # test rethrow with argument
    bt = []
    try
        try
            error("A")
        catch
            t = @task try
                error("B")
            catch exc
                exc
            end
            yield(t)
            @test t.state == :done
            @test t.result == ErrorException("B")
            bt = catch_backtrace()
            rethrow(ErrorException("C"))
        end
    catch exc
        @test exc == ErrorException("C")
        @test bt == catch_backtrace()
    end
    @test length(catch_stack()) == 0
    # Exception stacks on other tasks
    t = @task try
        error("A")
    catch
        error("B")
    end
    yield(t)
    @test t.state == :failed
    @test t.result == ErrorException("B")
    @test catch_stack(t, include_bt=false) == [ErrorException("A"), ErrorException("B")]
    # Exception stacks for tasks which never get the chance to start
    t = @task nothing
    @test (try
        @async Base.throwto(t, ErrorException("expected"))
        wait(t)
    catch e
        e
    end).task.exception == ErrorException("expected")
    @test length(catch_stack(t)) == 1
    @test length(catch_stack(t)[1][2]) > 0 # backtrace is nonempty
    # Exception stacks should not be accessed on concurrently running tasks
    t = @task ()->nothing
    @test_throws ErrorException("Inspecting the exception stack of a task which might "*
                                "be running concurrently isn't allowed.") catch_stack(t)
end

@testset "rethrow" begin
    @test try
        rethrow()
    catch exc
        exc
    end == ErrorException("rethrow() not allowed outside a catch block")
    @test try
        rethrow(ErrorException("A"))
    catch exc
        exc
    end == ErrorException("rethrow(exc) not allowed outside a catch block")
end

# issue #36527
function f36527()
    caught = false
    🏡 = Core.eval(Main, :(module asdf36527 end))
    try
        Core.eval(🏡, :(include_string($🏡, "@assert z36527 == 10")))
    catch ex
        GC.gc()
        catch_backtrace()
        caught = true
    end
    return caught
end

@test f36527()

# accessing an undefined var in tail position in a catch block
function undef_var_in_catch()
    try
        error("first error")
    catch
        __probably_n0t_defined__
    end
end
@test length(try
    undef_var_in_catch()
    []
catch
    catch_stack()
end) == 2

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Documenter
using Pkg

makedocs(
    modules = [Pkg],
    format = :html,
    sitename = "Pkg.jl",
    pages = Any[
        "index.md",
        "getting-started.md",
        "managing-packages.md",
        "environments.md",
        "creating-packages.md",
        "compatibility.md",
        "registries.md",
        "faq.md",
        "glossary.md",
        "api.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaLang/Pkg.jl",
    target = "build",
    julia = "nightly",
    deps = nothing,
    make = nothing,
)


text:
JSON.lower(vc::VirtualColumn) = non_nothing_dict(vc)

"""
    Expression(name::String, expression::String; outputType=nothing)

Create a virtual column using `expression`.

outputType should be String if provided.
"""
struct Expression <: VirtualColumn
    type::String
    name::String
    expression::String
    outputType
    function Expression(name, expression; outputType=nothing)
        outputType === nothing ||
            (outputType = uppercase(outputType)) ∈ ["LONG", "FLOAT", "DOUBLE", "STRING"] ||
            error("Invalid outputType")
        new("expression", name, expression, outputType)
    end
end


text:
using DifferentialEquations
using DataFrames
using Gadfly

a_=1.5
b_=1.0
c_=3.0
d_=1.0

f = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=a_ b=b_ c=c_ d=d_

function goLotka(u0=[1.0, 1.0], tspan=(0.0, 100.0), method=RK4)

  prob = ODEProblem(f, u0, tspan)
  sol = solve(prob, method(), dt=0.01)
  df=DataFrame(t=sol.t, x1=map(x->x[1],sol.u), y = map(x->x[2], sol.u), experiment="exp1")
  DataFrames.writetable("lotka.csv", df)

  plot(df,  x="x1", y="y")
end


text:
# wls object is not actively maintained

function calcWlsCoefficients(obj::WlsObject, inPt::Point)
  m = length(obj.b)
  datalen = size(obj.inputs, 1)
  firstTerm = zeros(m, m)
  secondTerm = zeros(m)

  for p in 1:datalen
    curPt = obj.inputs[p, :]
    w = obj.weightFunc(norm(inPt - curPt), obj.EPS)
    for i in 1:m
      for j in 1:m
        firstTerm[i, j] += w * obj.matrix[i, j](obj.vars => curPt)
      end
      secondTerm[i] += w * obj.b[i](obj.vars => curPt) * obj.outputs[p]
    end
  end

  coefficients = firstTerm \ secondTerm
  return coefficients
end

"""
    `(wo::WlsObject)(inputPoint::Point, value::Bool)`

Approximates the `WlsObject` at `inputPoint`.
"""
function (obj::WlsObject)(inPt::Point)
  c = calcWlsCoefficients(obj, inPt)
  poly = polynomial(c, obj.b)
  return poly(obj.vars => inPt)
end

function plotWls(obj::WlsObject)
  inputDim = size(obj.inputs, 2)
  if inputDim < 1 || inputDim > 2
    error("Cannot plot $inputDim dimensional graph")
  end
  x = [linspace(-2, 2, 21) for i in 1:inputDim]
  println(x)
  println(length(x))
  if inputDim == 2
    return plot(x[1], x[2], [obj([a, b], true) for b in x[1], a in x[2]], linetype=:surface)
  else
    return plot(x[1], [obj([a], true) for a in x[1]])
  end
end


text:
using FractionalDiffEq, Plots

# Analytical solution
analytical(x) = x.^1.8 .*mittleff(1.8, 2.8, -x.^1.8)
# Numerical solution
fun(x, y) = 1-y
prob = SingleTermFODEProblem(fun, 1.8, [0, 0], (0, 20))
sol = solve(prob, 0.01, PECE())
tspan = collect(0:0.01:20)
target = analytical.(tspan)

plot(sol, linewidth=5, label="Numerical", legend=:bottomright)
plot!(tspan, target, lw=3, ls=:dash, label="Analytical")


text:
using PackageCompiler

create_sysimage(:TOML, sysimage_path="./build/sys_toml.so", precompile_execution_file="./build/precompile_toml.jl")

max_stars_repo_name: periodicaidan/oryx-kbd-colors

max_issues_repo_name: periodicaidan/oryx-kbd-colors

max_forks_repo_name: periodicaidan/oryx-kbd-colors


text:
@testset "Mesh" begin
    tests = ["rectilinear_grid",
             "polytope_vertex_mesh",
             "./io"
            ]
    for t in tests
      include("$(t).jl")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using GaussianMixtureAlignment
using Test
using IntervalSets
using LinearAlgebra
using CoordinateTransformations

using GaussianMixtureAlignment: get_bounds, subranges, fullBlock

const GMA = GaussianMixtureAlignment
@testset "get bounds" begin
    μx = [3,0,0]
    μy = [-4,0,0]
    σ = ϕ = 1
    ndims = 3
    sqrt2 = √2

    x, y = IsotropicGaussian(μx, σ, ϕ), IsotropicGaussian(μy, σ, ϕ)

    # rotation distances, no translation
    # anti-aligned (no rotation) and aligned (180 degree rotation)
    lb, ub = get_bounds(x,y,2π,0,zeros(6))
    @test lb ≈ -GMA.overlap(1,2*σ^2,ϕ*ϕ, 1.) atol=1e-16
    @test ub ≈ -GMA.overlap(7^2,2*σ^2,ϕ*ϕ, 1.)
    lb, ub = get_bounds(x,y,2π,0,[0,0,π,0,0,0])
    @test lb ≈ ub ≈ -GMA.overlap(1,2*σ^2,ϕ*ϕ, 1.)
    # spheres with closest alignment at 90 degree rotation
    lb = get_bounds(x,y,π/√(3),0,zeros(6))[1]
    @test lb ≈ -GMA.overlap(5^2,2*σ^2,ϕ*ϕ, 1.)
    lb = get_bounds(x,y,π/(2*√(3)),0,[0,0,π/4,0,0,0])[1]
    @test lb ≈ -GMA.overlap(5^2,2*σ^2,ϕ*ϕ, 1.)

    # translation distance, no rotation
    # centered at origin
    lb, ub = get_bounds(x,y,0,2/√(3),zeros(6))
    @test lb ≈ -GMA.overlap(6^2,2*σ^2,ϕ*ϕ, 1.)
    @test ub ≈ -GMA.overlap(7^2,2*σ^2,ϕ*ϕ, 1.)
    # centered with translation of 1 in +x
    lb, ub = get_bounds(x,y,0,2/√(3),[0,0,0,1,0,0])
    @test lb ≈ -GMA.overlap(7^2,2*σ^2,ϕ*ϕ, 1.)
    @test ub ≈ -GMA.overlap(8^2,2*σ^2,ϕ*ϕ, 1.)
    # centered with translation of 3 in +y
    lb, ub = get_bounds(x,y,0,2/√(3),[0,0,0,0,3,0])
    @test lb ≈ -GMA.overlap((√(58)-1)^2,2*σ^2,ϕ*ϕ, 1.)
    @test ub ≈ -GMA.overlap(58,2*σ^2,ϕ*ϕ, 1.)

end

@testset "divide block" begin
    blk1 = NTuple{6,Tuple{Float64,Float64}}(((-π,π), (-π,π), (-π,π), (-1,1), (-1,1), (-1,1)))
    subblks1 = subranges(blk1, 2)
    @test length(subblks1) == 2^6
    for i=1:6
        intv = OpenInterval(0,0)
        for sblk in subblks1
            rng = sblk[i]
            intv = union(intv, ClosedInterval(rng[1], rng[2]))
        end
        @test intv == ClosedInterval(blk1[i][1], blk1[i][2])
    end
    # @show subblks1[length(subblks1)]

    blk2 = NTuple{6,Tuple{Float64,Float64}}(((0,π), (0,π), (0,π), (0,1), (0,1), (0,1)))
    subblks2 = subranges(blk2, 4)
    @test length(subblks2) == 4^6
    for i=1:6
        intv = OpenInterval(0,0)
        for sblk in subblks2
            rng = sblk[i]
            intv = union(intv, ClosedInterval(rng[1], rng[2]))
        end
        @test intv == ClosedInterval(blk2[i][1], blk2[i][2])
    end
end

@testset "GOGMA" begin
    # two sets of points, each forming a 3-4-5 triangle
    xpts = [[0.,0.,0.], [3.,0.,0.,], [0.,4.,0.]]
    ypts = [[1.,1.,1.], [1.,-2.,1.], [1.,1.,-3.]]
    σ = ϕ = 1.
    gmmx = IsotropicGMM([IsotropicGaussian(x, σ, ϕ) for x in xpts])
    gmmy = IsotropicGMM([IsotropicGaussian(y, σ, ϕ) for y in ypts])

    # aligning a GMM to itself
    bigblock = fullBlock(gmmx, gmmx)
    # @show bigblock.lowerbound, bigblock.upperbound
    @test bigblock.lowerbound ≈ -length(gmmx.gaussians)^2 # / √(4π)^3

    blk = fullBlock(gmmx, gmmx, NTuple{6,Tuple{Float64,Float64}}(((0,π), (0,π), (0,π), (0,2), (0,2), (0,2))))
    lb = blk.lowerbound
    ub = blk.upperbound
    for i = 1:20
        blk = fullBlock(gmmx, gmmx, subranges(blk.ranges, 2)[1])
        @test blk.lowerbound >= lb
        @test blk.upperbound <= ub
        lb = blk.lowerbound
        ub = blk.upperbound
    end

    # make sure this runs without an error
    res = gogma_align(gmmx, gmmy, maxblocks=1E5)
    res = tiv_gogma_align(gmmx, gmmy)

    mgmmx = IsotropicMultiGMM(Dict(:x => gmmx, :y => gmmy))
    mgmmy = IsotropicMultiGMM(Dict(:y => gmmx, :x => gmmy))
    res = gogma_align(mgmmx, mgmmy, maxblocks=1E5)
    res = tiv_gogma_align(mgmmx, mgmmy)

    # ROCS alignment should work perfectly
    @test isapprox(rocs_align(gmmx, gmmy).minimum, -overlap(gmmx,gmmx); atol=1E-12)
end

@testset "directional GOGMA" begin
    xpt = [1.,0.,0.]
    ypt = [1.,0.,0.]
    xdir = [0.,1.,0.]
    ydir = [1.,0.,0.]
    σ = ϕ = 1.
    x = IsotropicGaussian(xpt, σ, ϕ, [xdir])
    y = IsotropicGaussian(ypt, σ, ϕ, [ydir])
    @test get_bounds(x,y,0.,0.,zeros(6)) == (-0.5,-0.5)
    @test get_bounds(x,y,π/(3*GMA.sqrt3),0.,zeros(6)) == (-0.75,-0.5)
    @test get_bounds(x,y,π/GMA.sqrt3,0.,zeros(6)) == (-1.0,-0.5)

    xpt = [1.,0.,0.]
    ypt = [1.,0.,0.]
    xdir = [-1.,0.,0.]
    ydir = [1.,0.,0.]
    σ = ϕ = 1.
    x = IsotropicGaussian(xpt, σ, ϕ, [xdir])
    y = IsotropicGaussian(ypt, σ, ϕ, [ydir])
    @test get_bounds(x,y,0.,0.,zeros(6)) == (0.,0.)

    xpts = [[0.,0.,0.], [3.,0.,0.,], [0.,4.,0.]]
    ypts = [[1.,1.,1.], [1.,-2.,1.], [1.,1.,-3.]]
    σ = ϕ = 1.
    xdirs = [[0.,-1.,0], [1.,0.,0.], [0.,1.,0.]]
    ydirs = [[0.,0.,1.], [0.,-1.,0.], [0.,0.,-1]]
    dgmmx = IsotropicGMM([IsotropicGaussian(x, σ, ϕ, [xdirs[i]]) for (i,x) in enumerate(xpts)])
    dgmmy = IsotropicGMM([IsotropicGaussian(y, σ, ϕ, [ydirs[i]]) for (i,y) in enumerate(ypts)])
    objminxy = tiv_gogma_align(dgmmx, dgmmy).upperbound
    objminxx = tiv_gogma_align(dgmmx, dgmmx).upperbound
    objminyy = tiv_gogma_align(dgmmy, dgmmy).upperbound
    @test objminxy ≈ objminxx ≈ objminyy

    randdgmmy = IsotropicGMM([IsotropicGaussian(y, σ, ϕ, [2*rand(3).-1]) for (i,y) in enumerate(ypts)])
    randobjminxy = tiv_gogma_align(dgmmx, randdgmmy).upperbound
    @test randobjminxy > objminxy

    # when σ is very small, only very well aligned Gaussians will contribute to the overlap
    σ = 0.001
    ddgmmx = IsotropicGMM([IsotropicGaussian(x, σ, ϕ, [xdirs[i], normalize(2*rand(3).-1)]) for (i,x) in enumerate(xpts)])
    ddgmmy = IsotropicGMM([IsotropicGaussian(y, σ, ϕ, [ydirs[i]]) for (i,y) in enumerate(ypts)])
    ddobjmin = tiv_gogma_align(ddgmmx, ddgmmy).upperbound
    @test ddobjmin ≈ -3.0

    mgmmx = IsotropicMultiGMM(Dict(:one => IsotropicGMM([IsotropicGaussian(xpts[1], σ, ϕ, [xdirs[1]])]),
                                   :two => IsotropicGMM([IsotropicGaussian(xpts[2], σ, ϕ, [xdirs[2]])]),
                                   :three => IsotropicGMM([IsotropicGaussian(xpts[3], σ, ϕ, [xdirs[3]])])))
    mgmmy = IsotropicMultiGMM(Dict(:one => IsotropicGMM([IsotropicGaussian(ypts[1], σ, ϕ, [ydirs[1]])]),
                                   :two => IsotropicGMM([IsotropicGaussian(ypts[2], σ, ϕ, [ydirs[2]])]),
                                   :three => IsotropicGMM([IsotropicGaussian(ypts[3], σ, ϕ, [ydirs[3]])])))
    mobjminxy = tiv_gogma_align(mgmmx, mgmmy).upperbound
    @test mobjminxy ≈ -3.0
end

@testset "TIV-GOGMA (perfect alignment)" begin
    for i=1:10
        randpts = 10*rand(3,50)
        randtform = AffineMap(10*rand(6)...)
        gmmx = IsotropicGMM([IsotropicGaussian(randpts[:,i],1,1) for i=1:size(randpts,2)])
        gmmy = randtform(gmmx)
        min_overlap_score = -overlap(gmmx,gmmx)
        res = tiv_gogma_align(gmmx,gmmy,0.5,0.5; maxstagnant=1E3)
        @test isapprox(res.upperbound, min_overlap_score; rtol=0.01)
        @test isapprox(overlap(res.tform(gmmx), gmmy), -min_overlap_score; rtol=0.01)
    end

    for i=1:10
        randpts = 10*rand(3,50)
        randtform = AffineMap(10*rand(6)...)
        mgmmx = IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:Int(size(randpts,2)/5)]) for j=1:5]))
        mgmmy = randtform(mgmmx)
        min_overlap_score = -overlap(mgmmx,mgmmx)
        res = tiv_gogma_align(mgmmx,mgmmy,0.5,0.5; maxstagnant=1E3)
        @test isapprox(res.upperbound, min_overlap_score; rtol=0.01)
        @test isapprox(overlap(res.tform(mgmmx), mgmmy), -min_overlap_score; rtol=0.01)
    end
end

@testset "ROCS (perfect alignment)" begin
    for i=1:10
        randpts = 10*rand(3,50)
        randtform = AffineMap(10*rand(6)...)
        gmmx = IsotropicGMM([IsotropicGaussian(randpts[:,i],1,1) for i=1:size(randpts,2)])
        gmmy = randtform(gmmx)
        min_overlap_score = -overlap(gmmx,gmmx)
        rocs_res = rocs_align(gmmx,gmmy)
        ovlp, tform  = rocs_res.minimum, rocs_res.tform
        @test isapprox(ovlp, min_overlap_score; atol=1E-12)
        @test isapprox(overlap(tform(gmmx), gmmy), -min_overlap_score; atol=1E-12)
    end

    for i=1:10
        randpts = 10*rand(3,50)
        randtform = AffineMap(10*rand(6)...)
        mgmmx = IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:Int(size(randpts,2)/5)]) for j=1:5]))
        mgmmy = randtform(mgmmx)
        min_overlap_score = -overlap(mgmmx,mgmmx)
        rocs_res = rocs_align(mgmmx,mgmmy)
        ovlp, tform  = rocs_res.minimum, rocs_res.tform
        @test isapprox(ovlp, min_overlap_score; atol=1E-12)
        @test isapprox(overlap(tform(mgmmx), mgmmy), -min_overlap_score; atol=1E-12)
    end
end

# @testset "TIV-GOGMA (missing data alignment)" begin
#     for i=1:10
#         randpts = 10*rand(3,50)
#         randtform = AffineMap(10*rand(6)...)
#         gmmx = IsotropicGMM([IsotropicGaussian(randpts[:,i],1,1) for i=1:size(randpts,2)])
#         gmmy = randtform(IsotropicGMM(gmmx.gaussians[1:40]))
#         min_overlap_score = -overlap(gmmx,inv(randtform)(gmmy))
#         res = tiv_gogma_align(gmmx,gmmy,0.5,0.5; maxstagnant=1E4)
#         @test res.upperbound <= min_overlap_score*0.99
#         @test -overlap(res.tform(gmmx), gmmy) <= min_overlap_score*0.99
#     end

#     for i=1:10
#         randpts = 10*rand(3,50)
#         randtform = AffineMap(10*rand(6)...)
#         mgmmx = IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:10]) for j=1:5]))
#         mgmmy = randtform(IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:8]) for j=1:5])))
#         min_overlap_score = -overlap(mgmmx,inv(randtform)(mgmmy))
#         res = tiv_gogma_align(mgmmx,mgmmy,0.5,0.5; maxstagnant=1E4)
#         @test res.upperbound <= min_overlap_score*0.99
#         @test -overlap(res.tform(mgmmx), mgmmy) <= min_overlap_score*0.99
#     end
# end

# @testset "ROCS-GOGMA (missing data alignment)" begin
#     for i=1:10
#         randpts = 10*rand(3,50)
#         randtform = AffineMap(10*rand(6)...)
#         gmmx = IsotropicGMM([IsotropicGaussian(randpts[:,i],1,1) for i=1:size(randpts,2)])
#         gmmy = randtform(IsotropicGMM(gmmx.gaussians[1:40]))
#         min_overlap_score = -overlap(gmmx,inv(randtform)(gmmy))
#         res = rocs_gogma_align(gmmx,gmmy; maxstagnant=1E4)
#         @test res.upperbound <= min_overlap_score*0.99
#         @test -overlap(res.tform(gmmx), gmmy) <= min_overlap_score*0.99
#     end

#     for i=1:10
#         randpts = 10*rand(3,50)
#         randtform = AffineMap(10*rand(6)...)
#         mgmmx = IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:10]) for j=1:5]))
#         mgmmy = randtform(IsotropicMultiGMM(Dict([Symbol(j)=>IsotropicGMM([IsotropicGaussian(randpts[:,i+10*(j-1)],1,1) for i=1:8]) for j=1:5])))
#         min_overlap_score = -overlap(mgmmx,inv(randtform)(mgmmy))
#         res = rocs_gogma_align(mgmmx, mgmmy; maxstagnant=1E4)
#         @test res.upperbound <= min_overlap_score*0.99
#         @test -overlap(res.tform(mgmmx), mgmmy) <= min_overlap_score*0.99
#     end
# end


text:
using MinimaxCalibratedEBayes
using Test
using Plots

mhist =  MCEBHistogram(-6:0.02:6)
marginal_target = MarginalDensityTarget(StandardNormalSample(0.0))

marginal_target_cb = DiscretizedAffineEstimator(mhist, ButuceaComte(marginal_target, 10_000))

#test for symmetry
@test marginal_target_cb(-5.0001) == marginal_target_cb(+5.0001)

pgfplotsx()
plot(cb_disc)

prior_target = PriorDensityTarget(0.0)

prior_target_cb = DiscretizedAffineEstimator(mhist, ButuceaComte(prior_target, 10_000))

@test prior_target_cb(-5.0001) == prior_target_cb(+5.0001)

plot(prior_target_cb)

gmix_tmp = GaussianMixturePriorClass(0.1, [0.0], Gurobi.Optimizer)

Z = DiscretizedStandardNormalSample(1.0, mhist)

res = worst_case_bias(marginal_target_cb,
                 Z,
                 gmix_tmp,
                 marginal_target)

tmp = JuMP.value.(res[1][:πs])

findmax(tmp)

res2 = worst_case_bias(prior_target_cb,
                 Z,
                 gmix_tmp,
                 prior_target) #point is, I can check bias...

sum(prior_target_cb.Q .* prob_mass) - pdf(Normal(0.0, 0.1), 0.0)
prob_mass = pdf(marginalize(Normal(0.0, 0.1),Z))

sum(prob_mass)

integrate()

# check bias now.


text:
#=
Should not be constructed by the user - acts like a vector (used in eigsolve)
I think it makes sense to see these things as an actual state instead of return an array of B tensors (what we used to do)
This will allow us to plot energy density (finite qp) and measure observeables.
=#

struct FiniteQP{S<:FiniteMPS,T1,T2}
    # !(left_gs === right_gs) => domain wall excitation
    left_gs::S
    right_gs::S

    VLs::Vector{T1} # AL' VL = 0 (and VL*X = B)
    Xs::Vector{T2} # contains variational parameters

end

function rand_quasiparticle(left_gs::FiniteMPS,right_gs=left_gs;excitation_space=oneunit(virtualspace(left_gs,1)))
    #find the left null spaces for the TNS
    VLs = [adjoint(rightnull(adjoint(v))) for v in left_gs.AL]
    Xs = [TensorMap(rand,eltype(left_gs.AL[1]),space(VLs[loc],3)',excitation_space'*space(right_gs.AR[ loc],3)') for loc in 1:length(left_gs)]

    FiniteQP(left_gs,right_gs,VLs,Xs)
end

struct InfiniteQP{S<:InfiniteMPS,T1,T2}
    # !(left_gs === right_gs) => domain wall excitation
    left_gs::S
    right_gs::S

    VLs::Vector{T1} # AL' VL = 0 (and VL*X = B)
    Xs::Vector{T2} # contains variational parameters

    momentum::Float64
end

function rand_quasiparticle(left_gs::InfiniteMPS,right_gs=left_gs;excitation_space = oneunit(space(left_gs.AL[1],1)),momentum=0.0)
    #find the left null spaces for the TNS
    VLs = [adjoint(rightnull(adjoint(v))) for v in left_gs.AL]
    Xs = [TensorMap(rand,eltype(left_gs.AL[1]),space(VLs[loc],3)',excitation_space'*space(right_gs.AR[ loc+1],1)) for loc in 1:length(left_gs)]

    InfiniteQP(left_gs,right_gs,VLs,Xs,momentum)
end

const QP = Union{InfiniteQP,FiniteQP};

utilleg(v::QP) = space(v.Xs[1],2)
Base.copy(a::QP) = copyto!(similar(a),a)
function Base.copyto!(a::QP,b::QP)
    for (i,j) in zip(a.Xs,b.Xs)
        copyto!(i,j)
    end
    a
end
function Base.getproperty(v::QP,s::Symbol)
    if s == :trivial
        return v.left_gs === v.right_gs
    else
        return getfield(v,s)
    end
end

function Base.:-(v::QP,w::QP)
    t = similar(v)
    t.Xs[:] = (v.Xs-w.Xs)[:]
end
LinearAlgebra.dot(v::QP, w::QP) = sum(dot.(v.Xs, w.Xs))
LinearAlgebra.norm(v::QP) = norm(norm.(v.Xs))
LinearAlgebra.normalize!(w::QP) = rmul!(w,1/norm(w));
Base.length(v::QP) = length(v.Xs)
Base.eltype(v::QP) = eltype(eltype(v.Xs)) # - again debateable, need scaltype
Base.similar(v::InfiniteQP,t=eltype(v)) = InfiniteQP(v.left_gs,v.right_gs,v.VLs,map(e->similar(e,t),v.Xs),v.momentum)
Base.similar(v::FiniteQP,t=eltype(v)) = FiniteQP(v.left_gs,v.right_gs,v.VLs,map(e->similar(e,t),v.Xs))
function LinearAlgebra.mul!(w::QP, a, v::QP)
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.mul!(i, a, j)
    end
    return w
end

function LinearAlgebra.mul!(w::QP, v::QP, a)
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.mul!(i, j, a)
    end
    return w
end
function LinearAlgebra.rmul!(v::QP, a)
    for x in v.Xs
        LinearAlgebra.rmul!(x, a)
    end
    return v
end

function LinearAlgebra.axpy!(a, v::QP, w::QP)
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.axpy!(a, j, i)
    end
    return w
end
function LinearAlgebra.axpby!(a, v::QP, b, w::QP)
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.axpby!(a, j, b, i)
    end
    return w
end

Base.:*(v::QP, a) = mul!(similar(v),a,v)
Base.:*(a, v::QP) = mul!(similar(v),a,v)

Base.zero(v::QP) = v*0;
Base.getindex(v::QP,i::Int) = v.VLs[i]*v.Xs[i];
function Base.setindex!(v::QP,B,i::Int)
    v.Xs[i] = v.VLs[i]'*B
    v
end

function Base.convert(::Type{<:FiniteMPS},v::FiniteQP)
    #very slow and clunky, but shouldn't be performance critical anyway

    elt = eltype(v)
    function simplefuse(temp)
        frontmap = isomorphism(Matrix{elt},fuse(space(temp,1)*space(temp,2)),space(temp,1)*space(temp,2))
        backmap = isomorphism(Matrix{elt},space(temp,5)'*space(temp,4)',fuse(space(temp,5)'*space(temp,4)'))

        @tensor tempp[-1 -2;-3] := frontmap[-1,1,2]*temp[1,2,-2,3,4]*backmap[4,3,-3]
    end

    utl = utilleg(v); ou = oneunit(utl); utsp = ou ⊕ ou;
    upper = isometry(Matrix{elt},utsp,ou); lower = leftnull(upper);
    upper_I = upper*upper'; lower_I = lower*lower'; uplow_I = upper*lower';

    Ls = v.left_gs.AL[1:end];
    Rs = v.right_gs.AR[1:end];

    #step 0 : fuse the utility leg of B with the first leg of B
    Bs = map(1:length(v)) do i
        t = v[i]
        frontmap = isomorphism(Matrix{elt},fuse(utl*space(t,1)),utl*space(t,1));
        @tensor tt[-1 -2;-3]:=t[1,-2,2,-3]*frontmap[-1,2,1]
    end

    #step 1 : pass utl through Ls
    passer = isomorphism(Matrix{elt},utl,utl);
    Ls = map(Ls) do L
        @tensor temp[-1 -2 -3 -4;-5]:=L[-2,-3,-4]*passer[-1,-5]
        simplefuse(temp)
    end

    #step 2 : embed all Ls/Bs/Rs in the same space
    superspaces = map(zip(Ls,Rs)) do (L,R)
        supremum(space(L,1),space(R,1))
    end
    push!(superspaces,supremum(_lastspace(Ls[end])',_lastspace(Rs[end])'))
    for i in 1:(length(v)+1)
        Lf = isometry(Matrix{elt},superspaces[i],i <= length(v) ? _firstspace(Ls[i]) : _lastspace(Ls[i-1])')
        Rf = isometry(Matrix{elt},superspaces[i],i <= length(v) ? _firstspace(Rs[i]) : _lastspace(Rs[i-1])')

        if i <= length(v)
            @tensor Ls[i][-1 -2;-3] := Lf[-1,1]*Ls[i][1,-2,-3]
            @tensor Rs[i][-1 -2;-3] := Rf[-1,1]*Rs[i][1,-2,-3]
            @tensor Bs[i][-1 -2;-3] := Lf[-1,1]*Bs[i][1,-2,-3]
        end

        if i>1
            @tensor Ls[i-1][-1 -2;-3] := Ls[i-1][-1 -2;1]*conj(Lf[-3,1])
            @tensor Rs[i-1][-1 -2;-3] := Rs[i-1][-1 -2;1]*conj(Rf[-3,1])
            @tensor Bs[i-1][-1 -2;-3] := Bs[i-1][-1 -2;1]*conj(Rf[-3,1])
        end
    end

    #step 3 : fuse the correct *_I with the correct tensor (and enforce boundary conditions)
    for i in 1:length(v)
        @tensor temp[-1 -2 -3 -4; -5] := Ls[i][-2,-3,-4]*upper_I[-1,-5]
        Ls[i] = simplefuse(temp) * (i<length(v));

        @tensor temp[-1 -2 -3 -4; -5] := Rs[i][-2,-3,-4]*lower_I[-1,-5]
        Rs[i] = simplefuse(temp) * (i>1);

        @tensor temp[-1 -2 -3 -4; -5] := Bs[i][-2,-3,-4]*uplow_I[-1,-5]
        Bs[i] = simplefuse(temp);
    end

    return FiniteMPS(Ls+Rs+Bs)
end


text:
const OPTIONKEY = :ShowGraphviz_option

ShowGraphviz.getoption(io::IO) = get(io, OPTIONKEY, ``)::Cmd
ShowGraphviz.setoption(io::IO, option) = IOContext(io, OPTIONKEY, `$option`)
function ShowGraphviz.addoption(io::IO, option)
    old = ShowGraphviz.getoption(io)
    new = `$old $option`
    return ShowGraphviz.setoption(io, new)
end

function run_dot(output::IO, input::IO, option = ``)
    cfg = CONFIG.dot_option
    ctx = ShowGraphviz.getoption(output)
    dot() do cmd
        cmd = `$cmd $cfg $ctx $option`
        @debug "Run: $cmd"
        communicate(cmd; stdin = input, stdout = output)
    end
    return output
end

function dot_to_iobuffer(dot)
    io = IOBuffer()
    show(io, MIME"text/vnd.graphviz"(), dot)
    seekstart(io)
    return io
end

Base.show(io::IO, ::MIME"text/vnd.graphviz", dot::ShowGraphviz.DOT) = print(io, dot.source)
dot_to_iobuffer(dot::ShowGraphviz.DOT) = IOBuffer(dot.source)

"""
    ShowGraphviz.show(io::IO, mime, x; option = ``)

Given a Julia value `x` that defines `show` on `text/vnd.graphviz` MIME type,
convert it into mime type using `dot` program.  Command line options for `dot`
can be set by `option` keyword argument.

It can be used even when `show` of `x` is not defined with `@derive` or
`@deriveall`.
"""
ShowGraphviz.show(io::IO, mime::AbstractString, dot; option = ``) =
    ShowGraphviz.show(ShowGraphviz.addoption(io, option), MIME(mime), dot)

# Run `dot -T?`
# See also https://graphviz.org/docs/outputs/
const MIME_TO_DOT_OUTPUT = [
    "image/png" => "png",
    "image/gif" => "gif",
    "image/svg+xml" => "svg",
    "application/pdf" => "pdf",
]

for (mime, output) in MIME_TO_DOT_OUTPUT
    @eval ShowGraphviz.show(io::IO, ::$(typeof(MIME(mime))), dot) =
        run_dot(io, dot_to_iobuffer(dot), $(`-T$output`))
end

function derive_impl(T, mimes::Vector{String})
    exprs = map(mimes) do mime
        quote
            $Base.show(io::$IO, m::$(typeof(MIME(mime))), x::$T) =
                $ShowGraphviz.show(io, m, x)
        end
    end
    return Expr(:block, exprs...)
end

"""
    ShowGraphviz.@derive Type "mime₁" "mime₂" … "mimeₙ"
"""
macro derive(T, mimes...)
    mimes = collect(String, mimes)
    return esc(derive_impl(T, mimes))
end

"""
    ShowGraphviz.@deriveall Type
"""
macro deriveall(T)
    mimes = map(first, MIME_TO_DOT_OUTPUT)
    return esc(derive_impl(T, mimes))
end

macro dot_str(source)
    esc(:($ShowGraphviz.DOT($(QuoteNode(source)))))
end

function ShowGraphviz.DOT(source::IO)
    io = IOBuffer()
    write(io, source)
    return ShowGraphviz.DOT(String(take!(io)))
end

function ShowGraphviz.DOT(dot)
    io = dot_to_iobuffer(dot)
    return ShowGraphviz.DOT(String(take!(io)))
end

@deriveall ShowGraphviz.DOT

function Base.show(io::IO, ::MIME"text/plain", dot::ShowGraphviz.DOT)
    println(io, ShowGraphviz, ".dot\"\"\"")
    highlight_dot(io, dot.source)
    print(io, "\"\"\"")
end


text:
module DataAnalysis

include("FTUtils.jl")
include("HDF5Utils.jl")
include("TraceUtils.jl")
include("FilterUtils.jl")

end


text:
function simulate(pomcp::POMCPOWPlanner, h_node::POWTreeObsNode{B,A,O}, s::S, d) where {B,S,A,O}

    tree = h_node.tree
    h = h_node.node

    sol = pomcp.solver

    if POMDPs.isterminal(pomcp.problem, s) || d <= 0
        return (0.0,1)
    end

    if sol.enable_action_pw
        total_n = tree.total_n[h]
        if length(tree.tried[h]) <= sol.k_action*total_n^sol.alpha_action
            if h == 1
                a = next_action(pomcp.next_action, pomcp.problem, tree.root_belief, POWTreeObsNode(tree, h))
            else
                a = next_action(pomcp.next_action, pomcp.problem, StateBelief(tree.sr_beliefs[h]), POWTreeObsNode(tree, h))
            end
            if !sol.check_repeat_act || !haskey(tree.o_child_lookup, (h,a))
                push_anode!(tree, h, a,
                            init_N(pomcp.init_N, pomcp.problem, POWTreeObsNode(tree, h), a),
                            init_V(pomcp.init_V, pomcp.problem, POWTreeObsNode(tree, h), a),
                            sol.check_repeat_act)
            end
        end
    else # run through all the actions
        if isempty(tree.tried[h])
            if h == 1
                action_space_iter = POMDPs.actions(pomcp.problem, tree.root_belief)
            else
                action_space_iter = POMDPs.actions(pomcp.problem, StateBelief(tree.sr_beliefs[h]))
            end
            anode = length(tree.n)
            for a in action_space_iter
                push_anode!(tree, h, a,
                            init_N(pomcp.init_N, pomcp.problem, POWTreeObsNode(tree, h), a),
                            init_V(pomcp.init_V, pomcp.problem, POWTreeObsNode(tree, h), a),
                            false)
            end
        end
    end
    total_n = tree.total_n[h]

    best_node = select_best(pomcp.criterion, h_node, pomcp.solver.rng)
    a = tree.a_labels[best_node]

    new_node = false
    total_counts=0
    total_reward=0.0
    flag_csyes=false
    if tree.n_a_children[best_node] <= sol.k_observation*(tree.n[best_node]^sol.alpha_observation)

        sp, o, r = @gen(:sp, :o, :r)(pomcp.problem, s, a, sol.rng)

        if sol.check_repeat_obs && haskey(tree.a_child_lookup, (best_node,o))
            hao = tree.a_child_lookup[(best_node, o)]
        else
            new_node = true
            hao = length(tree.sr_beliefs) + 1
            push!(tree.sr_beliefs,
                  init_node_sr_belief(pomcp.node_sr_belief_updater,
                                      pomcp.problem, s, a, sp, o, r))
            """
            TODO: init_node是当o为新节点时，在创建新节点的同时把s'和w(s')加入到o中。
            因此在这我们需要进行补偿性采样，把s'和weight加入到其他的所有o中，然后在其他o中也进行补偿性采样。
            """
            """
            Step 1:把s'和weight加入到其他的所有o中.
            """
            """
            Step 2:从任意一个其他o中，取出所有的(s',r)对中的s'，计算w，加入到现在的o的belief中
            """
            if(length(tree.generated_ltc[best_node])!=0)#如果存在其他o，执行step1.2d<85 &&
                first_pair=first(tree.generated_ltc[best_node])
                first_pair_hao=first_pair.second
                if(length(tree.sr_beliefs[first_pair_hao].dist.items)<100) #step1
                    for pair in tree.generated_ltc[best_node]
                        o_temp=pair.first
                        hao_temp=pair.second
                        push_weighted!(tree.sr_beliefs[hao_temp], pomcp.node_sr_belief_updater, s, sp, r)
                    end
                end
            #    step2
                for i in 1:length(tree.sr_beliefs[first_pair_hao].dist.items)#step2
                    if(i>100)
                        break
                    end
                    sp_temp=tree.sr_beliefs[first_pair_hao].dist.items[i][1]
                    r_temp=tree.sr_beliefs[first_pair_hao].dist.items[i][2]
                    push_weighted!(tree.sr_beliefs[hao], pomcp.node_sr_belief_updater, s, sp_temp, r_temp)
                end
            end

            push!(tree.total_n, 0)
            push!(tree.tried, Int[])
            push!(tree.o_labels, o)

            if sol.check_repeat_obs
                tree.a_child_lookup[(best_node, o)] = hao
            end
            tree.n_a_children[best_node] += 1
            push!(tree.generated_ltc[best_node], o=>hao)
        end
        push!(tree.generated[best_node], o=>hao)
            """
            Step3:
            对其他每个o进行补偿性采样，
            """
            # if d<90
            #     if new_node
            #         for j in 1:(length(tree.generated_ltc[best_node])-1)
            #             pair=tree.generated_ltc[best_node][j]
            #             o_temp=pair.first
            #             hao_temp=pair.second
            #             before_weight=tree.sr_beliefs[hao_temp].dist.cdf[end-1]
            #             before_sample_times=tree.total_n[hao_temp]
            #             new_added_weight=tree.sr_beliefs[hao_temp].dist.cdf[end]-before_weight
            #             cs_times=div(before_sample_times*new_added_weight,before_weight)#下取整
            #             if (cs_times!=0)
            #                 for i in 1:cs_times
            #                     R_temp,counts_temp=simulate_ow(pomcp, POWTreeObsNode(tree, hao_temp), sp, d-1)
            #                     R = r + POMDPs.discount(pomcp.problem)*R_temp
            #                     total_counts+=counts_temp
            #                     total_reward+=R*counts_temp
            #                 end
            #             end
            #         end
            #     end
            # end
    else

        sp, r = @gen(:sp, :r)(pomcp.problem, s, a, sol.rng)

    end

    if r == Inf
        @warn("POMCPOW: +Inf reward. This is not recommended and may cause future errors.")
    end

    if new_node
        R = r + POMDPs.discount(pomcp.problem)*estimate_value(pomcp.solved_estimate, pomcp.problem, sp, POWTreeObsNode(tree, hao), d-1)
        total_counts+=1
        total_reward+=R
        """
        total_counts++
        total_reward+=R
        """
    else
        pair = rand(sol.rng, tree.generated[best_node])
        o = pair.first
        hao = pair.second
        """
        To-do: 当o不是新节点时，需要把这轮从G中得到的s'，加入到当前a下的所有o分支的B(hao)中，同时根据每个分支的o计算相应的权重
        然后在每个O分支下，对新粒子进行补偿性采样。
        """
        """
        Step 1:把s'和weight加入到所有o中,包括自己.
        """
        if(length(tree.sr_beliefs[hao].dist.items)<100) #限制粒子数为1000
            for pair in tree.generated_ltc[best_node]
                o_temp=pair.first
                hao_temp=pair.second
                push_weighted!(tree.sr_beliefs[hao_temp], pomcp.node_sr_belief_updater, s, sp, r)
            end
            """
            Step2:补偿性采样
            """
            #  for pair in tree.generated_ltc[best_node]
            #      o_temp=pair.first
            #      hao_temp=pair.second
            #      before_weight=tree.sr_beliefs[hao_temp].dist.cdf[end-1]
            #      before_sample_times=tree.total_n[hao_temp]
            #      new_added_weight=tree.sr_beliefs[hao_temp].dist.cdf[end]-tree.sr_beliefs[hao_temp].dist.cdf[end-1]
            #      cs_times=div(before_sample_times*new_added_weight,before_weight)#下取整
            #      if (cs_times!=0)
            #          for i in 1:cs_times
            #              R_temp,counts_temp=simulate(pomcp, POWTreeObsNode(tree, hao_temp), sp, d-1)
            #              R = r + POMDPs.discount(pomcp.problem)*R_temp
            #              total_reward+=R*counts_temp
            #              total_counts+=counts_temp
            #          end
            #      end
            #  end
        else
            push_weighted!(tree.sr_beliefs[hao], pomcp.node_sr_belief_updater, s, sp, r)#默认把一个粒子加入到当前o
        end
      #  push_weighted!(tree.sr_beliefs[hao], pomcp.node_sr_belief_updater, s, sp, r)#默认把一个粒子加入到当前o
        temp_lenght=length(tree.sr_beliefs[hao].dist.items)
        nums_v=length(tree.total_n[hao])
        if(d<90&&nums_v>1&&nums_v<10)#temp_lenght<10) #d<30
            before_weight=tree.sr_beliefs[hao].dist.cdf[temp_lenght-1]
            before_sample_times=tree.total_n[hao]
            new_added_weight=tree.sr_beliefs[hao].dist.cdf[temp_lenght]-before_weight
            cs_times=div(before_sample_times*new_added_weight,before_weight)#下取整
            if (cs_times!=0)
                flag_csyes=true
                if(cs_times>3)
                    cs_times=2
                end
                for i in 1:cs_times
                  #  sp, r = rand(sol.rng, tree.sr_beliefs[hao])
                    R_temp,counts_temp=simulate(pomcp, POWTreeObsNode(tree, hao), sp, d-1)
                    R = r + POMDPs.discount(pomcp.problem)*R_temp
                    total_counts+=counts_temp
                    total_reward+=R*counts_temp
                end
            end
        end

       """
       Step3:随机采样

       """
       sp, r = rand(sol.rng, tree.sr_beliefs[hao])
       if(flag_csyes)
       R_temp,counts_temp=simulate(pomcp, POWTreeObsNode(tree, hao), sp, d-1)
       else
       R_temp,counts_temp=simulate(pomcp, POWTreeObsNode(tree, hao), sp, d-1)
       end

       R = r + POMDPs.discount(pomcp.problem)*R_temp
       total_counts+=counts_temp
       total_reward+=R*counts_temp
        # sp, r = rand(sol.rng, tree.sr_beliefs[hao])
        # R = r + POMDPs.discount(pomcp.problem)*simulate(pomcp, POWTreeObsNode(tree, hao), sp, d-1)
    end

    """
    返回Reward和次数
    """
    tree.n[best_node] += total_counts
    tree.total_n[h] += total_counts
    if tree.v[best_node] != -Inf
        tree.v[best_node] = ((tree.n[best_node]-total_counts)*tree.v[best_node]+total_reward)/tree.n[best_node]
    end

    return (total_reward/total_counts,total_counts)

    # tree.n[best_node] += 1
    # tree.total_n[h] += 1
    # if tree.v[best_node] != -Inf
    #     tree.v[best_node] += (R-tree.v[best_node])/tree.n[best_node]
    # end
    # return R
end

function simulate_ow(pomcp::POMCPOWPlanner, h_node::POWTreeObsNode{B,A,O}, s::S, d) where {B,S,A,O}

    tree = h_node.tree
    h = h_node.node

    sol = pomcp.solver

    if POMDPs.isterminal(pomcp.problem, s) || d <= 0
        return (0.0,1)
    end

    if sol.enable_action_pw
        total_n = tree.total_n[h]
        if length(tree.tried[h]) <= sol.k_action*total_n^sol.alpha_action
            if h == 1
                a = next_action(pomcp.next_action, pomcp.problem, tree.root_belief, POWTreeObsNode(tree, h))
            else
                a = next_action(pomcp.next_action, pomcp.problem, StateBelief(tree.sr_beliefs[h]), POWTreeObsNode(tree, h))
            end
            if !sol.check_repeat_act || !haskey(tree.o_child_lookup, (h,a))
                push_anode!(tree, h, a,
                            init_N(pomcp.init_N, pomcp.problem, POWTreeObsNode(tree, h), a),
                            init_V(pomcp.init_V, pomcp.problem, POWTreeObsNode(tree, h), a),
                            sol.check_repeat_act)
            end
        end
    else # run through all the actions
        if isempty(tree.tried[h])
            if h == 1
                action_space_iter = POMDPs.actions(pomcp.problem, tree.root_belief)
            else
                action_space_iter = POMDPs.actions(pomcp.problem, StateBelief(tree.sr_beliefs[h]))
            end
            anode = length(tree.n)
            for a in action_space_iter
                push_anode!(tree, h, a,
                            init_N(pomcp.init_N, pomcp.problem, POWTreeObsNode(tree, h), a),
                            init_V(pomcp.init_V, pomcp.problem, POWTreeObsNode(tree, h), a),
                            false)
            end
        end
    end
    total_n = tree.total_n[h]

    best_node = select_best(pomcp.criterion, h_node, pomcp.solver.rng)
    a = tree.a_labels[best_node]

    new_node = false
    total_counts=0
    total_reward=0.0
    if tree.n_a_children[best_node] <= sol.k_observation*(tree.n[best_node]^sol.alpha_observation)

        sp, o, r = @gen(:sp, :o, :r)(pomcp.problem, s, a, sol.rng)

        if sol.check_repeat_obs && haskey(tree.a_child_lookup, (best_node,o))
            hao = tree.a_child_lookup[(best_node, o)]
        else
            new_node = true
            hao = length(tree.sr_beliefs) + 1
            push!(tree.sr_beliefs,
                  init_node_sr_belief(pomcp.node_sr_belief_updater,
                                      pomcp.problem, s, a, sp, o, r))

            push!(tree.total_n, 0)
            push!(tree.tried, Int[])
            push!(tree.o_labels, o)

            if sol.check_repeat_obs
                tree.a_child_lookup[(best_node, o)] = hao
            end
            tree.n_a_children[best_node] += 1

        end
        push!(tree.generated[best_node], o=>hao)

    else

        sp, r = @gen(:sp, :r)(pomcp.problem, s, a, sol.rng)

    end

    if r == Inf
        @warn("POMCPOW: +Inf reward. This is not recommended and may cause future errors.")
    end

    if new_node
        R = r + POMDPs.discount(pomcp.problem)*estimate_value(pomcp.solved_estimate, pomcp.problem, sp, POWTreeObsNode(tree, hao), d-1)
        total_counts+=1
        total_reward+=R
        """
        total_counts++
        total_reward+=R
        """
    else
        pair = rand(sol.rng, tree.generated[best_node])
        o = pair.first
        hao = pair.second

        push_weighted!(tree.sr_beliefs[hao], pomcp.node_sr_belief_updater, s, sp, r)#默认把一个粒子加入到当前o

       """
       Step3:随机采样

       """
       sp, r = rand(sol.rng, tree.sr_beliefs[hao])

       R_temp,counts_temp=simulate_ow(pomcp, POWTreeObsNode(tree, hao), sp, d-1)

       R = r + POMDPs.discount(pomcp.problem)*R_temp
       total_counts+=counts_temp
       total_reward+=R*counts_temp
        # sp, r = rand(sol.rng, tree.sr_beliefs[hao])
        # R = r + POMDPs.discount(pomcp.problem)*simulate(pomcp, POWTreeObsNode(tree, hao), sp, d-1)
    end

    """
    返回Reward和次数
    """
    tree.n[best_node] += total_counts
    tree.total_n[h] += total_counts
    if tree.v[best_node] != -Inf
        tree.v[best_node] = ((tree.n[best_node]-total_counts)*tree.v[best_node]+total_reward)/tree.n[best_node]
    end

    return (total_reward/total_counts,total_counts)

    # tree.n[best_node] += 1
    # tree.total_n[h] += 1
    # if tree.v[best_node] != -Inf
    #     tree.v[best_node] += (R-tree.v[best_node])/tree.n[best_node]
    # end
    # return R
end

max_stars_repo_name: Lttcc/CSPOMCPOW

max_issues_repo_name: Lttcc/CSPOMCPOW

max_forks_repo_name: Lttcc/CSPOMCPOW


text:
using PackageCompiler
using Glob

using TurbulenceConvection
using CalibrateEDMF
using EnsembleKalmanProcesses
using OrdinaryDiffEq
using NCDatasets
using StochasticDiffEq
using ClimaCore

cedmf = pkgdir(CalibrateEDMF)

if isempty(Glob.glob("CEDMF.so"))
    create_sysimage(
        [
            "TurbulenceConvection",
            "EnsembleKalmanProcesses",
            "OrdinaryDiffEq",
            "NCDatasets",
            "StochasticDiffEq",
            "ClimaCore",
        ];
        sysimage_path = "CEDMF.so",
        precompile_execution_file = joinpath(cedmf, "test", "runtests.jl"),
    )
end


text:
"""Utilities for testing search structures."""
module Testing

using Test
using Neighborhood
using Neighborhood: bruteforcesearch

export cmp_search_results, cmp_bruteforce, search_allfuncs, check_search_results,
       test_bulksearch

"""
Get arguments tuple to `search`, using the 3-argument version if `skip=nothing`.
"""
get_search_args(ss, query, t, skip) = isnothing(skip) ? (ss, query, t) : (ss, query, t, skip)

"""
    cmp_search_results(results...)::Bool

Compare two or more sets of search results (`(idxs, ds)` tuples) and check that
they are identical up to ordering.
"""
function cmp_search_results(results::Tuple{Vector, Vector}...)
    length(results) < 2 && error("Expected at least two sets of results")

    idxs1, ds1 = results[1]
    rest = results[2:end]

    idxset = Set(idxs1)
    dist_map = Dict(i => d for (i, d) in zip(idxs1, ds1))

    for (idxs_i, ds_i) in rest
        Set(idxs_i) == idxset || return false
        all(dist_map[i] == d for (i, d) in zip(idxs_i, ds_i)) || return false
    end

    return true
end

"""
    cmp_bruteforce(results, data, metric, query, t[, skip])::Bool

Check whether `results` returned from [`search`](@ref) match those computed
with [`Neighborhood.bruteforcesearch`](@ref)`(data, metric, query, t[, skip])`
(up to order). `skip` may be `nothing`, which calls the 4-argument method.

**Caution:** results of a `NeighborNumber` search are only expected to match if
the distances from `query` to each point in `data` are all distinct, otherwise
there may be some ambiguity in which data points are included.
"""
function cmp_bruteforce(results, data, metric, query, t, skip=nothing)
    bf = bruteforcesearch(data, metric, query, t, skip)
    return cmp_search_results(results, bf)
end

"""
    search_allfuncs(ss, query, t[, skip])

Call [`search`](@ref)`(ss, query, t[, skip])` and check that the result matches
those for [`isearch`](@ref) and [`knn`](@ref)/[`inrange`](@ref) (depending on
search type) for the equivalent arguments.

`skip` may be `nothing`, in which case the 3-argument methods of all functions
will be called. Uses `Test.@test` internally.
"""
function search_allfuncs(ss, query, t, skip=nothing)
    args = get_search_args(ss, query, t, skip)
    idxs, ds = result = search(args...)

    @test Set(isearch(args...)) == Set(idxs)
    cmp_search_results(result, _alt_search_func(args...))

    return result
end

# Call inrange() or knn() given arguments to search()
function _alt_search_func(ss, query, t::WithinRange, args...)
    inrange(ss, query, t.r, args...)
end
function _alt_search_func(ss, query, t::NeighborNumber, args...)
    knn(ss, query, t.k, args...)
end

"""
    check_search_results(data, metric, results, query, t[, skip])

Check that `results = search(ss, query, t[, skip])` make sense for a search
structure `ss` with data `data` and metric `metric`.

Note that this does not calculate the known correct value to compare to (which
may be expensive for large data sets), just that the results have the
expected properties. Use [`cmp_bruteforce`] for a more exact test.
`skip` may be `nothing`, in which case the 3-argument methods of all functions
will be called. Uses `Test.@test` internally.

Checks the following:
* `results` is a 2-tuple of `(idxs, ds)`.
* `ds[i] == metric(query, data[i])`.
* `skip(i)` is false for all `i` in `idxs`.
* For `t::NeighborNumber`:
  * `length(idxs) <= t.k`.
* For `t::WithinRange`:
  * `d <= t.r` for all `d` in `ds`.
"""
function check_search_results(data, metric, results, query, t, skip=nothing)
    idxs, ds = results
    @test ds == [metric(query, data[i]) for i in idxs]
    !isnothing(skip) && @test !any(map(skip, idxs))
    _check_search_results(data, metric, results, query, t, skip)
end

function _check_search_results(data, metric, (idxs, ds), query, t::NeighborNumber, skip)
    @test length(idxs) <= t.k
end

function _check_search_results(data, metric, (idxs, ds), query, t::WithinRange, skip)
    @test all(<=(t.r), ds)
end

"""
    test_bulksearch(ss, queries, t[, skip=nothing])

Test that [`bulksearch`](@ref) gives the same results as individual applications
of [`search`](@ref).

`skip` may be `nothing`, in which case the 3-argument methods of both functions
will be called. Uses `Test.@test` internally.
"""
function test_bulksearch(ss, queries, t, skip=nothing)
    args = get_search_args(ss, queries, t, skip)
    bidxs, bds = bulksearch(args...)

    @test bulkisearch(args...) == bidxs

    for (i, query) in enumerate(queries)
        result = if isnothing(skip)
            search(ss, query, t)
        else
            iskip = j -> skip(i, j)
            search(ss, query, t, iskip)
        end
        @test result == (bidxs[i], bds[i])
    end
end

end  # module


text:
module StdLib

using ..CxxWrapCore

abstract type CppBasicString <: AbstractString end

# These are defined in C++, but the functions need to exist to add methods
function append end
function cppsize end
function cxxgetindex end
function cxxsetindex! end
function push_back end
function resize end

import Libdl
@wrapmodule(CxxWrapCore.libcxxwrap_julia_stl, :define_cxxwrap_stl_module, Libdl.RTLD_GLOBAL)

function __init__()
  @initcxx
end

# Pass-through for fundamental types
_append_dispatch(v::StdVector,a::Vector,::Type{CxxWrapCore.IsNormalType}) = append(v,a)
# For C++ types, convert the array to an array of references, so the pointers can be read directly from a contiguous array on the C++ side
_append_dispatch(v::StdVector{T}, a::Vector{<:T},::Type{CxxWrapCore.IsCxxType}) where {T} = append(v,CxxWrapCore.CxxRef.(a))
# Choose the correct append method depending on the type trait
append(v::StdVector{T}, a::Vector{<:T}) where {T} = _append_dispatch(v,a,CxxWrapCore.cpp_trait_type(T))

Base.ncodeunits(s::CppBasicString)::Int = cppsize(s)
Base.codeunit(s::StdString) = UInt8
Base.codeunit(s::StdWString) = Cwchar_t == Int32 ? UInt32 : UInt16
Base.codeunit(s::CppBasicString, i::Integer) = reinterpret(codeunit(s), s[i])
Base.isvalid(s::CppBasicString, i::Integer) = (0 < i <= ncodeunits(s))
function Base.iterate(s::CppBasicString, i::Integer=1)
  if !isvalid(s,i)
    return nothing
  end
  return(convert(Char,codeunit(s,i)),i+1)
end

function StdWString(s::String)
  char_arr = transcode(Cwchar_t, s)
  StdWString(char_arr, length(char_arr))
end

function StdVector(v::Vector{T}) where {T}
  if (CxxWrapCore.cpp_trait_type(T) == CxxWrapCore.IsCxxType)
    return StdVector(CxxRef.(v))
  end
  result = StdVector{T}()
  append(result, v)
  return result
end

function StdVector(v::Vector{CxxRef{T}}) where {T}
  result = isconcretetype(T) ? StdVector{supertype(T)}() : StdVector{T}()
  append(result, v)
  return result
end

function StdVector(v::Vector{Bool})
  result = StdVector{CxxBool}()
  append(result, convert(Vector{CxxBool}, v))
  return result
end

Base.IndexStyle(::Type{<:StdVector}) = IndexLinear()
Base.size(v::StdVector) = (Int(cppsize(v)),)
Base.getindex(v::StdVector, i::Int) = cxxgetindex(v,i)[]
Base.setindex!(v::StdVector{T}, val, i::Int) where {T} = cxxsetindex!(v, convert(T,val), i)

function Base.push!(v::StdVector, x)
  push_back(v, x)
  return v
end

function Base.resize!(v::StdVector, n::Integer)
  resize(v, n)
  return v
end

Base.empty!(v::StdVector) = Base.resize!(v, 0)

function Base.append!(v::StdVector, a::Vector)
  append(v, a)
  return v
end

function Base.append!(v::StdVector{CxxBool}, a::Vector{Bool})
  append(v, convert(Vector{CxxBool}, a))
  return v
end

Base.String(s::StdString) = unsafe_string(reinterpret(Ptr{Cchar},c_str(s).cpp_object))
function Base.String(s::StdWString)
  chars = unsafe_wrap(Vector{Cwchar_t}, reinterpret(Ptr{Cwchar_t},c_str(s).cpp_object), (cppsize(s),))
  return transcode(String, chars)
end
Base.show(io::IO, s::CppBasicString) = show(io, String(s))
Base.cmp(a::CppBasicString, b::String) = cmp(String(a),b)
Base.cmp(a::String, b::CppBasicString) = cmp(a,String(b))

# Make sure functions taking a C++ string as argument can also take a Julia string
CxxWrapCore.map_julia_arg_type(x::Type{<:StdString}) = AbstractString
Base.convert(::Type{T}, x::String) where {T<:StdString} = StdString(x)
Base.cconvert(::Type{CxxWrapCore.ConstCxxRef{StdString}}, x::String) = StdString(x)
Base.unsafe_convert(::Type{CxxWrapCore.ConstCxxRef{StdString}}, x::StdString) = ConstCxxRef(x)

function StdValArray(v::Vector{T}) where {T}
  return StdValArray{T}(v, length(v))
end

Base.IndexStyle(::Type{<:StdValArray}) = IndexLinear()
Base.size(v::StdValArray) = (Int(cppsize(v)),)
Base.getindex(v::StdValArray, i::Int) = cxxgetindex(v,i)[]
Base.setindex!(v::StdValArray{T}, val, i::Int) where {T} = cxxsetindex!(v, convert(T,val), i)

end


text:
function NewtonR()
    println("Using Newton's method!!!!!")
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.QCircuits.Gates: appendparams!, ParamT, appendRandParams!,
                                       U3

function getparam(g)
    param = ParamT[]
    appendparams!(param, g)

    return param
end
function getRandParam(g)
    param = ParamT[]
    appendRandParams!(param, g)

    return param
end

qr = QuantumRegister(1)
g = U3(qr[0], 1.0, 2.0, 3.0)
@test getparam(g) == ParamT[]
@test getRandParam(g) == ParamT[]
@test length(g) == 0
@test_throws AssertionError setparameters!(g, [2.0, 1.0, 4.0])
@test_throws AssertionError setparameters!(g, [2.0])
setparameters!(g, [])

g = U3(qr[0], ParameterT(1.0), ParameterT(2.0), ParameterT(3.0))
@test getparam(g) == ParamT[1.0, 2.0, 3.0]
@test length(getRandParam(g)) == 3
@test getparam(g) != getRandParam(g)
@test length(g) == 3
setparameters!(g, [2.0, 1.0, 4.0])
@test getparam(g) == ParamT[2.0, 1.0, 4.0]
@test_throws AssertionError setparameters!(g, [2.0, 1.0, 4.0, 4.0])
@test_throws AssertionError setparameters!(g, [2.0, 1.0])
setparameters!(g, [-2.0, -1.0, -4.0])
@test getparam(g) == ParamT[-2.0, -1.0, -4.0]

g = U3(qr[0])
@test length(getparam(g)) == 3
@test length(getRandParam(g)) == 3
@test getparam(g) != getRandParam(g)
@test length(g) == 3
setparameters!(g, [2.0, 1.0, 4.0])
@test getparam(g) == ParamT[2.0, 1.0, 4.0]

qr = QuantumRegister(1)
g = U3(qr[0], ParameterT(1.0), ParameterT(2.0), 3.0)
@test getparam(g) == ParamT[1.0, 2.0]
@test length(getRandParam(g)) == 2
@test getparam(g) != getRandParam(g)
@test length(g) == 2
setparameters!(g, [2.0, 1.0])
@test getparam(g) == ParamT[2.0, 1.0]
@test_throws AssertionError setparameters!(g, [2.0, 1.0, 4.0])

################################################################################
#                Complex Params                                                #
################################################################################
using QuantumCircuits.QCircuits.ComplexGates
using QuantumCircuits.QCircuits.ComplexGates: U4_params

qr = QuantumRegister(2)
g = U4(qr[0], qr[1])
p = getparam(g)
@test length(p) == U4_params
rp = getRandParam(g)
@test length(rp) == U4_params
@test p != rp
@test length(g) == U4_params
np = [i for i in 1:U4_params]
@test p != np
setparameters!(g, np)
@test getparam(g) == np

qc = QCircuit(2)
qc.u4(0, 1)
p = getparameters(qc)
@test length(p) == U4_params
rp = getRandParameters(qc)
@test length(rp) == U4_params
@test p != rp
np = [i for i in 1:U4_params]
@test p != np
setparameters!(qc, np)
@test getparameters(qc) == np

dqc = decompose(qc)
@test length(getparameters(dqc)) == U4_params
@test getparameters(dqc) == np
np = [-i for i in 1:U4_params]
setparameters!(dqc, np)
@test getparameters(dqc) == np
bindparameters!(dqc)
@test getparameters(dqc) == Number[]

# @test length(getparameters(qc)) == U4_params
# bindparameters!(qc)
# @test getparameters(qc) == Number[]

################################################################################

qr = QuantumRegister(3, "q")
qc = QCircuit(qr)
qc.u3(qr)
n = 3*3
p = getparameters(qc)
@test length(p) == n
rp = getRandParameters(qc)
@test length(rp) == n
@test p != rp
setparameters!(qc, rp)
@test getparameters(qc) == rp

np = [i for i in 1:n]
@test rp != np
setparameters!(qc, np)
@test getparameters(qc) == np

bindparameters!(qc)
@test getparameters(qc) == Number[]


text:
@testset "ArviZ.Dataset" begin
    data = load_arviz_data("centered_eight")
    dataset = data.posterior
    @test dataset isa ArviZ.Dataset

    @testset "construction" begin
        pydataset = PyObject(dataset)
        @test ArviZ.Dataset(pydataset) isa ArviZ.Dataset
        @test PyObject(ArviZ.Dataset(pydataset)) === pydataset
        @test ArviZ.Dataset(dataset) === dataset
        @test_throws ArgumentError ArviZ.Dataset(py"PyNullObject()")
        @test hash(dataset) == hash(pydataset)

        vars = Dict("x" => ("dimx", randn(3)), ("y" => (("dimy_1", "dimy_2"), randn(3, 2))))
        coords = Dict(
            "dimx" => [1, 2, 3], "dimy_1" => ["a", "b", "c"], "dimy_2" => ["d", "e"]
        )
        attrs = Dict("prop1" => 1, "prop2" => "propval")
        @inferred ArviZ.Dataset(; data_vars=vars, coords, attrs)
        ds = ArviZ.Dataset(; data_vars=vars, coords, attrs)
        @test ds isa ArviZ.Dataset
        vars2, kwargs = ArviZ.dataset_to_dict(ds)
        for (k, v) in vars
            @test k ∈ keys(vars2)
            @test vars2[k] ≈ v[2]
        end
        @test kwargs.coords == coords
        for (k, v) in attrs
            @test k ∈ keys(kwargs.attrs)
            @test kwargs.attrs[k] == v
        end
    end

    @testset "properties" begin
        @test length(propertynames(dataset)) > 1
        @test dataset["mu"] isa PyObject
        @test dataset["mu"].values == py"$(dataset)['mu'].values"
    end

    @testset "conversion" begin
        @test pyisinstance(PyObject(dataset), ArviZ.xarray.Dataset)
        dataset4 = convert(ArviZ.Dataset, PyObject(dataset))
        @test dataset4 isa ArviZ.Dataset
        @test PyObject(dataset4) === PyObject(dataset)

        # TODO: improve this test
        @test convert(ArviZ.Dataset, [1.0, 2.0, 3.0, 4.0]) isa ArviZ.Dataset
    end

    @testset "show(::ArviZ.Dataset)" begin
        @testset "$mimetype" for mimetype in ("plain", "html")
            text = repr(MIME("text/$(mimetype)"), dataset)
            @test text isa String
            @test !occursin("<xarray.Dataset>", text)
            @test !occursin("<xarray.Dataset>", text)
            @test occursin("Dataset (xarray.Dataset)", text)
        end
    end
end

@testset "ArviZ.convert_to_dataset" begin
    rng = MersenneTwister(42)

    @testset "ArviZ.convert_to_dataset(::ArviZ.Dataset; kwargs...)" begin
        data = load_arviz_data("centered_eight")
        dataset = data.posterior
        @test ArviZ.convert_to_dataset(dataset) isa ArviZ.Dataset
        @test ArviZ.convert_to_dataset(dataset) === dataset
    end

    @testset "ArviZ.convert_to_dataset(::InferenceData; kwargs...)" begin
        A = Dict("A" => randn(rng, 2, 10, 2))
        B = Dict("B" => randn(rng, 2, 10, 2))
        dataA = ArviZ.convert_to_dataset(A)
        dataB = ArviZ.convert_to_dataset(B)
        idata = InferenceData(; posterior=dataA, prior=dataB)

        ds1 = ArviZ.convert_to_dataset(idata)
        @test ds1 isa ArviZ.Dataset
        @test "A" ∈ [ds1.keys()...]

        ds2 = ArviZ.convert_to_dataset(idata; group=:prior)
        @test ds2 isa ArviZ.Dataset
        @test "B" ∈ [ds2.keys()...]
    end

    @testset "ArviZ.convert_to_dataset(::Dict; kwargs...)" begin
        data = Dict("x" => randn(rng, 4, 100), "y" => randn(rng, 4, 100))
        dataset = ArviZ.dict_to_dataset(data)
        @test dataset isa ArviZ.Dataset
        @test "x" ∈ [dataset.keys()...]
        @test "y" ∈ [dataset.keys()...]
    end
end

@testset "ArviZ.convert_to_constant_dataset" begin
    @testset "ArviZ.convert_to_constant_dataset(::Dict)" begin
        data = Dict("x" => randn(4, 5), "y" => ["a", "b", "c"])
        dataset = ArviZ.convert_to_constant_dataset(data)
        @test dataset isa ArviZ.Dataset
        @test "x" ∈ dataset.keys()
        @test "y" ∈ dataset.keys()
        @test Set(dataset.coords) == Set(["x_dim_0", "x_dim_1", "y_dim_0"])
        @test collect(dataset._variables["x"].values) == data["x"]
        @test collect(dataset._variables["y"].values) == data["y"]
    end

    @testset "ArviZ.convert_to_constant_dataset(::Dict; kwargs...)" begin
        data = Dict("x" => randn(4, 5), "y" => ["a", "b", "c"])
        coords = Dict("xdim1" => 1:4, "xdim2" => 5:9, "ydim1" => ["d", "e", "f"])
        dims = Dict("x" => ["xdim1", "xdim2"], "y" => ["ydim1"])
        library = "MyLib"
        dataset = ArviZ.convert_to_constant_dataset(data)
        attrs = Dict("prop" => "propval")

        dataset = ArviZ.convert_to_constant_dataset(data; coords, dims, library, attrs)
        @test dataset isa ArviZ.Dataset
        @test "x" ∈ dataset.keys()
        @test "y" ∈ dataset.keys()
        @test Set(dataset.coords) == Set(["xdim1", "xdim2", "ydim1"])
        @test collect(dataset._variables["xdim1"].values) == coords["xdim1"]
        @test collect(dataset._variables["xdim2"].values) == coords["xdim2"]
        @test collect(dataset._variables["ydim1"].values) == coords["ydim1"]
        @test collect(dataset["x"].coords) == ["xdim1", "xdim2"]
        @test collect(dataset["y"].coords) == ["ydim1"]
        @test collect(dataset._variables["x"].values) == data["x"]
        @test collect(dataset._variables["y"].values) == data["y"]
        @test dataset.attrs["prop"] == attrs["prop"]
        @test dataset.attrs["inference_library"] == library
    end

    @testset "ArviZ.convert_to_constant_dataset(::NamedTuple; kwargs...)" begin
        data = (x=randn(4, 5), y=["a", "b", "c"])
        coords = (xdim1=1:4, xdim2=5:9, ydim1=["d", "e", "f"])
        dims = (x=["xdim1", "xdim2"], y=["ydim1"])
        library = "MyLib"
        dataset = ArviZ.convert_to_constant_dataset(data)
        attrs = (prop="propval",)

        dataset = ArviZ.convert_to_constant_dataset(data; coords, dims, library, attrs)
        @test dataset isa ArviZ.Dataset
        @test "x" ∈ dataset.keys()
        @test "y" ∈ dataset.keys()
        @test Set(dataset.coords) == Set(["xdim1", "xdim2", "ydim1"])
        @test collect(dataset._variables["xdim1"].values) == coords.xdim1
        @test collect(dataset._variables["xdim2"].values) == coords.xdim2
        @test collect(dataset._variables["ydim1"].values) == coords.ydim1
        @test collect(dataset["x"].coords) == ["xdim1", "xdim2"]
        @test collect(dataset["y"].coords) == ["ydim1"]
        @test collect(dataset._variables["x"].values) == data.x
        @test collect(dataset._variables["y"].values) == data.y
        @test dataset.attrs["prop"] == attrs.prop
        @test dataset.attrs["inference_library"] == library
    end
end

@testset "dict to dataset roundtrip" begin
    rng = MersenneTwister(42)
    J = 8
    K = 6
    nchains = 4
    ndraws = 500
    vars = Dict(
        "a" => randn(rng, nchains, ndraws), "b" => randn(rng, nchains, ndraws, J, K)
    )
    coords = Dict("bi" => 1:J, "bj" => 1:K)
    dims = Dict("b" => ["bi", "bj"])
    attrs = Dict("mykey" => 5)

    ds = ArviZ.dict_to_dataset(vars; library=:MyLib, coords, dims, attrs)
    @test ds isa ArviZ.Dataset
    vars2, kwargs = ArviZ.dataset_to_dict(ds)
    for (k, v) in vars
        @test k ∈ keys(vars2)
        @test vars2[k] ≈ v
    end
    @test kwargs.coords == coords
    @test kwargs.dims == dims
    for (k, v) in attrs
        @test k ∈ keys(kwargs.attrs)
        @test kwargs.attrs[k] == v
    end
    @test "inference_library" ∈ keys(kwargs.attrs)
    @test kwargs.attrs["inference_library"] == "MyLib"
end


text:
#=
multigrid:
- Julia version: 0.6.2
- Author: jerrymei
- Date: 2018-05-22
=#

mutable struct gridWrapper{T}
    u::Array{T,2}
    newu::Array{T,2}
    dx::Array{T,2}
    dy::Array{T,2}
    rhs::Array{T,2}
    f::Array{T,2}
    fft::rFFTWrapper
end

function gridWrapper{T}(f::AbstractArray{T,2})
    fft=rFFTWrapper(f)
    u=zeros(f)
    newu=similar(u)
    dx=similar(u)
    dy=similar(u)
    rhs=similar(u)
    return gridWrapper(u,newu,dx,dy,rhs,f,fft)
end

function restriction!{T}(fine::AbstractArray{T,2}, coarse::AbstractArray{T,2})
    nx,ny=size(coarse)
    @inbounds for x in 1:nx
        mx=2 * x - 1
        lx=mx > 1?mx - 1:2 * nx
        ux=mx + 1
        for y in 1:ny
            my=2 * y - 1
            ly=my > 1?my - 1:2 * ny
            uy=my + 1
            coarse[x,y]=(fine[lx,ly] + 2. * fine[mx,ly] + fine[ux,ly])/16.+ (fine[lx,my] + 2. * fine[mx,my] + fine[ux,my])/8.+ (fine[lx,uy] + 2. * fine[mx,uy] + fine[ux,uy]) / 16.
        end
    end
end

function interpolation!{T}(coarse::AbstractArray{T,2}, fine::AbstractArray{T,2})
    nx,ny =size(coarse)
    @inbounds for x in 1:nx
        mx=2 * x - 1
        ux=x < nx?x + 1:1
        for y in 1:ny
            my=2 * y - 1
            uy=y < ny?y + 1:1
            fine[mx,my]=coarse[x,y]
            fine[mx,my + 1]=(coarse[x,y] + coarse[x,uy]) / 2.
            fine[mx + 1,my]=(coarse[x,y] + coarse[ux,y]) / 2.
            fine[mx + 1,my + 1]=(coarse[x,y] + coarse[x,uy] + coarse[ux,y] + coarse[ux,uy]) / 4.
        end
    end
end

function full_cycle{T}(iteration,grid1,grid2,grid3,grid4,eps::T,tol::T)
    iteration(grid4,eps)
    interpolation!(grid4.u,grid3.u)
    iteration(grid3,eps)
    restriction!(grid3.u,grid4.u)

    iteration(grid4,eps)
    interpolation!(grid4.u,grid3.u)
    iteration(grid3,eps)
    interpolation!(grid3.u,grid2.u)
    iteration(grid2,eps)
    restriction!(grid2.u,grid3.u)
    iteration(grid3,eps)
    restriction!(grid3.u,grid4.u)

    iteration(grid4,eps)
    interpolation!(grid4.u,grid3.u)
    iteration(grid3,eps)
    interpolation!(grid3.u,grid2.u)
    iteration(grid2,eps)
    interpolation!(grid2.u,grid1.u)
end

max_stars_repo_name: JerryLingjieMei/NonlinearPDE

max_issues_repo_name: JerryLingjieMei/NonlinearPDE

max_forks_repo_name: JerryLingjieMei/NonlinearPDE


text:
function assemble_shortest_paths(g::GenericGraph, s, dist_key; distance::Bool=false)

    nodes = [k.index => k for k = vertices(g)]

    Source_nodes = Int64[]
    Predecessors = [k.index => Int64[] for k = vertices(g)] # vymenene {} za []
    sigma = [vertex.index => 0.0 for vertex = vertices(g)]
    Distances = Dict{Int64, Float64}()
    sigma[s] = 1.0
    visited = Dict(s => 0)
    c = 0
    Q = (Array{Any, 1})[]
    Collections.heappush!(Q, [0, c += 1 , s, s],Base.Order.ord((x,y)->isless(x[1], y[1]), identity, false))

    while length(Q) > 0

        heap_element = Collections.heappop!(Q, Base.Order.ord((x,y)->isless(x[1], y[1]), identity, false))
        dist, _, pred, vertex = heap_element[1] , heap_element[2] , heap_element[3] , heap_element[4]

        if haskey(Distances, vertex)
            continue
        end

        sigma[vertex] += sigma[pred]
        push!(Source_nodes,vertex)
        Distances[vertex] = dist

        for edge in out_edges(nodes[vertex], g)

            vw_dist = dist + edge.attributes[dist_key]

            if (haskey(Distances, edge.target.index) == false) && ( haskey(visited, edge.target.index) == false || vw_dist < visited[edge.target.index] )

                visited[edge.target.index] = vw_dist
                Collections.heappush!(Q, [vw_dist, c += 1, vertex, edge.target.index],Base.Order.ord((x,y)->isless(x[1], y[1]), identity, false))
                sigma[edge.target.index] = 0.0
                Predecessors[edge.target.index] = [vertex]

            elseif vw_dist == visited[edge.target.index]

                sigma[edge.target.index] += sigma[vertex]
                push!(Predecessors[edge.target.index],vertex)

            end

        end

    end

    if distance
      return Source_nodes, Predecessors, Distances
    elseif !(distance)
      return Source_nodes, Predecessors, sigma
    end

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
import Random

@enum Fathomed none dominance infeasibility optimality

include("dataStruct.jl")
include("pretrait.jl")
include("branchAndBound.jl")

function firstPhase(prob::Problem, params...)

    @assert length(params)==1 "The parameters must exactly one : ϵ"

    ϵ = params[1]
    XSEm = []
    S = []
    x1 = solve1OKP(weightedScalarRelax(prob, [1., ϵ]))
    x2 = solve1OKP(weightedScalarRelax(prob, [ϵ, 1.]))

    zx1 = evaluate(prob, x1)
    zx2 = evaluate(prob, x2)

    if zx1 == zx2
        push!(XSEm, x1)
    else
        push!(S, (zx1, zx2))
        push!(XSEm, x1, x2)
        while length(S) != 0
            yr, yl = pop!(S)
            λ = [yl[2]-yr[2], yr[1]-yl[1]]
            xe = solve1OKP(weightedScalarRelax(prob, λ))

            zxe = evaluate(prob, xe)

            if sum(λ .* zxe) > sum(λ .* yr)
                push!(XSEm, xe)
                push!(S, (yr, zxe), (zxe, yl))
            end
        end
    end
    return XSEm

end

function secondPhase(prob::Problem, params...)

    @assert length(params) == 1 "There must be 1 parameters to secondPhase : Name of the order"

end

function jules(fname::String = "test.dat", methodName::String = "TwoPhases", params...)
    prob = parser(fname)
    testBandB(prob)
end

function main(fname::String = "1.dat"; verbose = false)

    prob = parser(fname)

    sol = solve1OKPMain(weightedScalarRelax(prob,[1.,0.]), verbose = verbose)

    return sol
end


text:
using Documenter, ExcelFiles

makedocs(modules=[ExcelFiles],
	sitename="ExcelFiles.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ])

deploydocs(repo="github.com/queryverse/ExcelFiles.jl.git")


text:
include("consts.jl")
include("params.jl")
include("equations.jl")
include("GERG2008/GERG2008.jl")
include("EOS_LNG/EOS_LNG.jl")


text:
# NOTE: PlanarAzimuthalSymmetry() does not included all possible wavenumbers
function wavenumbers_bisection_robust(ω::T, medium::PhysicalMedium{T,Dim}, species::Species{T,Dim};
        symmetry::AbstractSetupSymmetry{Dim} = PlanarAzimuthalSymmetry{Dim}(),
        tol::T = 1e-5,
        num_wavenumbers = 3,
        box_k::Vector{Vector{T}} = box_keff(ω, medium, species; tol = tol),
        bisection_mesh_points::Int = Int(round(2 - log(tol))) + 3 * num_wavenumbers,
        bisection_iteration::Int = Int(round(-log(tol) / 3)),
        fixedparameters::Optim.FixedParameters = NelderMeadparameters(),
        optimoptions::Optim.Options{T} = Optim.Options(
            g_tol = tol^T(3), x_tol=tol^T(2),
            iterations = 1000 + Int(round(log(tol)*20))
        ),
        kws...) where {T,Dim}

    # exact low-frequency effective wavenumber is known exactly
    ko = real(ω / effective_medium(medium, species).c)

    disp = dispersion_complex(ω, medium, species, symmetry; kws...)
    # disp = dispersion_complex(ω, medium, species, symmetry; basis_order = basis_order)

    freal(x, y)::T = real(disp(x + y*im))
    fimag(x, y)::T = imag(disp(x + y*im))

    x = LinRange(box_k[1][1],box_k[1][2],bisection_mesh_points);
    y = LinRange(box_k[2][1],box_k[2][2], Int(round(bisection_mesh_points/2)));

    prob_real = MDBM_Problem(freal,[Axis(x,"x"),Axis(y,"y")]);
    prob_imag = MDBM_Problem(fimag,[Axis(x,"x"),Axis(y,"y")]);

    solve!(prob_real,bisection_iteration);
    solve!(prob_imag,bisection_iteration);

    # x_eval,y_eval=getevaluatedpoints(prob_real)
    xr_sol,yr_sol=getinterpolatedsolution(prob_real);
    xi_sol,yi_sol=getinterpolatedsolution(prob_imag)

    x_max = maximum(abs.([xr_sol; xi_sol]))
    y_max = maximum(abs.([yr_sol; yi_sol]))

    x_normalise = sqrt(length(xr_sol)) / x_max
    y_normalise = sqrt(length(yr_sol)) / y_max

    vecs1 = map(xr_sol,yr_sol) do x,y
        round.([x * x_normalise, y * y_normalise])
    end
    vecs2 = map(xi_sol,yi_sol) do x,y
        round.([x * x_normalise, y * y_normalise])
    end

    root_vecs = collect(intersect(Set(vecs1),Set(vecs2)))

    inds1 = [
        findall([v == root for v in vecs1])
    for root in root_vecs]

    inds2 = [
        findall([v == root for v in vecs2])
    for root in root_vecs]

    roots = map(eachindex(root_vecs)) do i
        x = mean([xr_sol[inds1[i]];xi_sol[inds2[i]]])
        y = mean([yr_sol[inds1[i]];yi_sol[inds2[i]]])
        [x, y]
    end
    roots = reduce_kvecs(roots, T(50) * tol * ko)

    # If the mesh is not fine enough,
    roots = sort(roots, by = v -> v[2])
    roots = roots[1:min(length(roots),50*num_wavenumbers)]

    # using Plots
    #points where the function foo was evaluated
    # x_eval,y_eval=getevaluatedpoints(prob_real)
    # scatter(x_eval,y_eval, label="real", markersize = 1.0, markerstrokewidth = 0.0)

    # scatter(xr_sol,yr_sol, label="real", markersize = 1.0, markerstrokewidth = 0.0)
    # scatter!(xi_sol,yi_sol, label="imag", markersize = 1.0, markerstrokewidth = 0.0)
    # scatter!([r[1] for r in roots], [r[2] for r in roots], markersize = 2.0)

    f_vec(x_vec)::T = abs(disp(x_vec[1] + x_vec[2]*im))

    # Used to create the NelderMead simplexer
    dx = T(0.5) * x_max / sqrt(length(roots))
    dy = T(0.5) * y_max / sqrt(length(roots))

    # Here we refine the roots
    roots2 = map(roots) do root

        inner_optimizer = NelderMead(
            initial_simplex =  MySimplexer(dx,dy),
            parameters  = fixedparameters
        )

        res = optimize(f_vec, root, inner_optimizer, optimoptions)

        if res.minimum < T(100)*tol || (Optim.converged(res) && res.minimum < T(1000)*tol)
            res.minimizer
        else
            [zero(T),-one(T)]
        end
    end
    roots2 = reduce_kvecs(roots2, T(50) * tol * ko)

    k_effs = [sol[1] + sol[2]*im for sol in roots2]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end


text:
using GeoScanning
using Test

@testset "GeoScanning.jl" begin
    # Write your own tests here.
end