text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(boottestOLSARubin!),StrBootTest{Float32}})   # time: 8.2192955
    Base.precompile(Tuple{typeof(boottestOLSARubin!),StrBootTest{Float64}})   # time: 6.797035
    Base.precompile(Tuple{typeof(boottestWRE!),StrBootTest{Float64}})   # time: 4.739786
    Base.precompile(Tuple{typeof(boottestWRE!),StrBootTest{Float32}})   # time: 4.028054
    Base.precompile(Tuple{Core.kwftype(typeof(wildboottest)),NamedTuple{(:resp, :predexog, :clustid, :auxwttype), Tuple{Vector{Float32}, Matrix{Float64}, Vector{Int32}, Symbol}},typeof(wildboottest),Type,Matrix{Int64},Vector{Float64}})   # time: 1.4766428
    Base.precompile(Tuple{Core.kwftype(typeof(__wildboottest)),NamedTuple{(:resp, :predexog, :predendog, :inst, :R1, :r1, :clustid, :nbootclustvar, :nerrclustvar, :issorted, :hetrobust, :nfe, :feid, :fedfadj, :obswt, :fweights, :maxmatsize, :ptype, :bootstrapc, :liml, :fuller, :kappa, :arubin, :small, :clusteradj, :clustermin, :scorebs, :reps, :imposenull, :auxwttype, :rng, :level, :rtol, :madjtype, :nH0, :ml, :scores, :beta, :A, :gridmin, :gridmax, :gridpoints, :diststat, :getci, :getplot, :getauxweights), Tuple{Vector{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Vector{Float32}, Matrix{Int64}, Int64, Int64, Bool, Bool, Int64, Vector{Int64}, Int64, Vector{Float32}, Bool, Float16, Symbol, Bool, Bool, Float32, Float32, Bool, Bool, Bool, Bool, Bool, Int64, Bool, Symbol, MersenneTwister, Float32, Float32, Symbol, Int16, Bool, Matrix{Float32}, Vector{Float32}, Symmetric{Float32, Matrix{Float32}}, Vector{Float32}, Vector{Float32}, Vector{Float32}, Symbol, Bool, Bool, Bool}},typeof(__wildboottest),Matrix{Float32},Vector{Float32}})   # time: 0.9131017
    Base.precompile(Tuple{Core.kwftype(typeof(__wildboottest)),NamedTuple{(:resp, :predexog, :predendog, :inst, :R1, :r1, :clustid, :nbootclustvar, :nerrclustvar, :issorted, :hetrobust, :nfe, :feid, :fedfadj, :obswt, :fweights, :maxmatsize, :ptype, :bootstrapc, :liml, :fuller, :kappa, :arubin, :small, :clusteradj, :clustermin, :scorebs, :reps, :imposenull, :auxwttype, :rng, :level, :rtol, :madjtype, :nH0, :ml, :scores, :beta, :A, :gridmin, :gridmax, :gridpoints, :diststat, :getci, :getplot, :getauxweights), Tuple{Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}, Matrix{Int64}, Int64, Int64, Bool, Bool, Int64, Vector{Int64}, Int64, Vector{Float64}, Bool, Float16, Symbol, Bool, Bool, Float64, Float64, Bool, Bool, Bool, Bool, Bool, Int64, Bool, Symbol, MersenneTwister, Float64, Float64, Symbol, Int16, Bool, Matrix{Float64}, Vector{Float64}, Symmetric{Float64, Matrix{Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Symbol, Bool, Bool, Bool}},typeof(__wildboottest),Matrix{Float64},Vector{Float64}})   # time: 0.7967473
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple{(:resp, :clustid, :auxwttype, :predexog), Tuple{Vector{Float32}, Vector{Int32}, Symbol, Matrix{Float64}}},typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64}})   # time: 0.2241874
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Int64}})   # time: 0.0757203
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple,typeof(_wildboottest),DataType,Matrix{Float64},Vector{Int64}})   # time: 0.0512666
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple{(:resp, :reps, :scorebs, :predexog, :hetrobust), Tuple{Vector{Float32}, Int64, Bool, Matrix{Float64}, Bool}},typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64}})   # time: 0.0399314
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Vector{Float32},Vector{Float32},SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0370134
    Base.precompile(Tuple{Core.kwftype(typeof(wildboottest)),NamedTuple{(:resp, :predexog, :hetrobust, :scorebs), Tuple{Vector{Float32}, Matrix{Float64}, Bool, Bool}},typeof(wildboottest),Type,Matrix{Int64},Vector{Float64}})   # time: 0.036235
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple{(:resp, :inst, :clustid, :reps, :small, :predendog, :predexog), Tuple{Vector{Float32}, Matrix{Int8}, Vector{Int8}, Int64, Bool, Vector{Float32}, Matrix{Float64}}},typeof(_wildboottest),DataType,Matrix{Int64},Vector{Int64}})   # time: 0.0359046
    Base.precompile(Tuple{Core.kwftype(typeof(wildboottest)),NamedTuple{(:resp, :predexog, :hetrobust, :scorebs, :reps), Tuple{Vector{Float32}, Matrix{Float64}, Bool, Bool, Int64}},typeof(wildboottest),Type,Matrix{Int64},Vector{Float64}})   # time: 0.0341064
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float64},Matrix{Float64},SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0334703
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float32},Matrix{Float32},SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0327042
    Base.precompile(Tuple{Core.kwftype(typeof(wildboottest)),NamedTuple{(:resp, :predexog, :predendog, :inst, :clustid, :small, :reps), Tuple{Vector{Float32}, Matrix{Float64}, Vector{Float32}, Matrix{Int8}, Vector{Int8}, Bool, Int64}},typeof(wildboottest),Type,Matrix{Int64},Vector{Int64}})   # time: 0.0324644
    Base.precompile(Tuple{Core.kwftype(typeof(wildboottest)),NamedTuple{(:resp, :predexog, :predendog, :inst, :clustid, :small, :reps, :imposenull, :getplot, :ptype), Tuple{Vector{Float32}, Matrix{Float64}, Vector{Float32}, Matrix{Int8}, Vector{Int8}, Bool, Int64, Bool, Bool, Symbol}},typeof(wildboottest),Type,Matrix{Int64},Vector{Int64}})   # time: 0.0324147
    let fbody = try __lookup_kwbody__(which(_wildboottest, (DataType,Matrix{Int64},Vector{Float64},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Vector{Float32},Matrix{Float64},Matrix{Float32},Matrix{Float32},Matrix{Float32},Vector{Float32},Bool,Matrix{Int64},Int64,Int64,Bool,Int64,Vector{Int8},Int64,Vector{Float32},Bool,Int64,Symbol,Bool,Bool,Int64,Float64,Bool,Bool,Bool,Bool,Bool,Int64,Bool,Symbol,MersenneTwister,Float64,Float64,Symbol,Int64,Bool,Matrix{Float32},Vector{Float32},Matrix{Float32},Vector{Float32},Vector{Float32},Vector{Int64},Symbol,Bool,Bool,Bool,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64},))
    end
end   # time: 0.0299475
    let fbody = try __lookup_kwbody__(which(_wildboottest, (DataType,Matrix{Int64},Vector{Int64},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Vector{Float32},Matrix{Float64},Vector{Float32},Matrix{Int8},Matrix{Float64},Vector{Float64},Bool,Vector{Int8},Int64,Int64,Bool,Int64,Vector{Int8},Int64,Vector{Float64},Bool,Int64,Symbol,Bool,Bool,Int64,Float64,Bool,Bool,Bool,Bool,Bool,Int64,Bool,Symbol,MersenneTwister,Float64,Float64,Symbol,Int64,Bool,Matrix{Float32},Vector{Float64},Matrix{Float64},Vector{Float64},Vector{Float64},Vector{Int64},Symbol,Bool,Bool,Bool,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Int64},))
    end
end   # time: 0.026748
    let fbody = try __lookup_kwbody__(which(_wildboottest, (DataType,Matrix{Int64},Vector{Int64},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Vector{Float32},Matrix{Float64},Vector{Float32},Matrix{Int8},Matrix{Float32},Vector{Float32},Bool,Vector{Int8},Int64,Int64,Bool,Int64,Vector{Int8},Int64,Vector{Float32},Bool,Int64,Symbol,Bool,Bool,Int64,Float64,Bool,Bool,Bool,Bool,Bool,Int64,Bool,Symbol,MersenneTwister,Float64,Float64,Symbol,Int64,Bool,Matrix{Float32},Vector{Float32},Matrix{Float32},Vector{Float32},Vector{Float32},Vector{Int64},Symbol,Bool,Bool,Bool,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Int64},))
    end
end   # time: 0.0262361
    let fbody = try __lookup_kwbody__(which(_wildboottest, (DataType,Matrix{Int64},Vector{Float64},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Vector{Float32},Matrix{Float64},Matrix{Float64},Matrix{Float64},Matrix{Float64},Vector{Float64},Bool,Vector{Int32},Int64,Int64,Bool,Int64,Vector{Int8},Int64,Vector{Float64},Bool,Int64,Symbol,Bool,Bool,Int64,Float64,Bool,Bool,Bool,Bool,Bool,Int64,Bool,Symbol,MersenneTwister,Float64,Float64,Symbol,Int64,Bool,Matrix{Float32},Vector{Float64},Matrix{Float64},Vector{Float64},Vector{Float64},Vector{Int64},Symbol,Bool,Bool,Bool,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64},))
    end
end   # time: 0.0259131
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true},Matrix{Float64},Vector{Float64},Vector{UnitRange{Int64}}})   # time: 0.0253607
    let fbody = try __lookup_kwbody__(which(_wildboottest, (DataType,Matrix{Int64},Vector{Float64},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Vector{Float32},Matrix{Float64},Matrix{Float32},Matrix{Float32},Matrix{Float32},Vector{Float32},Bool,Vector{Int32},Int64,Int64,Bool,Int64,Vector{Int8},Int64,Vector{Float32},Bool,Int64,Symbol,Bool,Bool,Int64,Float64,Bool,Bool,Bool,Bool,Bool,Int64,Bool,Symbol,MersenneTwister,Float64,Float64,Symbol,Int64,Bool,Matrix{Float32},Vector{Float32},Matrix{Float32},Vector{Float32},Vector{Float32},Vector{Int64},Symbol,Bool,Bool,Bool,typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64},))
    end
end   # time: 0.0253141
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Vector{Float64},Vector{Float64},SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0238532
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),SubArray{Float32, 2, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true},Matrix{Float32},Vector{Float32},Vector{UnitRange{Int64}}})   # time: 0.022691
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true},Matrix{Float64},SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0226257
    Base.precompile(Tuple{typeof(vecconvert),DataType,BitVector})   # time: 0.0223217
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),SubArray{Float32, 2, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true},Matrix{Float32},SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Vector{UnitRange{Int64}}})   # time: 0.0221149
    Base.precompile(Tuple{typeof(matconvert),DataType,Matrix{Int8}})   # time: 0.0216603
    Base.precompile(Tuple{typeof(matconvert),Type{Int64},Matrix{Int16}})   # time: 0.0211393
    Base.precompile(Tuple{typeof(matconvert),DataType,Vector{Int32}})   # time: 0.0209911
    Base.precompile(Tuple{typeof(vecconvert),DataType,Vector{Int16}})   # time: 0.0208953
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Matrix{Float64},Int64,Matrix{Float64},Symmetric{Float64, Matrix{Float64}},Matrix{Float64}})   # time: 0.0207441
    Base.precompile(Tuple{typeof(vecconvert),DataType,Vector{Int8}})   # time: 0.020406
    Base.precompile(Tuple{typeof(matconvert),DataType,Matrix{Int16}})   # time: 0.0191304
    Base.precompile(Tuple{typeof(matconvert),Type{Int64},Matrix{Int8}})   # time: 0.0187837
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Matrix{Float32},Int64,Matrix{Float32},Symmetric{Float32, Matrix{Float32}},Matrix{Float32}})   # time: 0.0172341
    Base.precompile(Tuple{typeof(vecconvert),Type{Int64},Vector{Int16}})   # time: 0.0142607
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),SubArray{Float64, 2, Array{Float64, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true},Int64,Matrix{Float64},Matrix{Float64},Matrix{Float64}})   # time: 0.0126324
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Adjoint{Float64, Vector{Float64}},Int64,Matrix{Float64},Matrix{Float64},Matrix{Float64}})   # time: 0.0094184
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Adjoint{Float32, Vector{Float32}},Int64,Matrix{Float32},Matrix{Float32},Matrix{Float32}})   # time: 0.0093277
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),SubArray{Float32, 2, Array{Float32, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true},Int64,Matrix{Float32},Matrix{Float32},Matrix{Float32}})   # time: 0.0087927
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple{(:resp, :scorebs, :predexog, :hetrobust), Tuple{Vector{Float32}, Bool, Matrix{Float64}, Bool}},typeof(_wildboottest),DataType,Matrix{Int64},Vector{Float64}})   # time: 0.0080437
    Base.precompile(Tuple{Core.kwftype(typeof(_wildboottest)),NamedTuple{(:resp, :inst, :clustid, :reps, :imposenull, :getplot, :small, :predendog, :ptype, :predexog), Tuple{Vector{Float32}, Matrix{Int8}, Vector{Int8}, Int64, Bool, Bool, Bool, Vector{Float32}, Symbol, Matrix{Float64}}},typeof(_wildboottest),DataType,Matrix{Int64},Vector{Int64}})   # time: 0.0079362
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Matrix{Float64},Int64,Matrix{Float64},Matrix{Float64},Matrix{Float64}})   # time: 0.0074214
    Base.precompile(Tuple{typeof(colquadformminus_nonturbo!),Matrix{Float32},Int64,Matrix{Float32},Matrix{Float32},Matrix{Float32}})   # time: 0.0071387
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float32},Matrix{Float32},Vector{Float32},Vector{UnitRange{Int64}}})   # time: 0.0066268
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float64},Matrix{Float64},Vector{Float64},Vector{UnitRange{Int64}}})   # time: 0.0055284
    Base.precompile(Tuple{typeof(matmulplus_nonturbo!),Vector{Float64},Matrix{Float64},Vector{Float64}})   # time: 0.0041138
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float32},Matrix{Float32},Vector{UnitRange{Int64}}})   # time: 0.0040009
    Base.precompile(Tuple{typeof(matmulplus_nonturbo!),Vector{Float32},Matrix{Float32},Vector{Float32}})   # time: 0.0038271
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Matrix{Float64},Matrix{Float64},Vector{UnitRange{Int64}}})   # time: 0.0037496
    Base.precompile(Tuple{Type{BootTestResult{Float64}},Float64,String,Float64,Float64,Int64,Int64,Int64,Int64,Float64,NamedTuple{(:X, :p), Tuple{Tuple{Vector{Float64}, Vector{Float64}}, Vector{Float64}}},NamedTuple{(:X, :p), Tuple{Vector{Float64}, Float64}},Matrix{Float64},Matrix{Float64},Vector{Float64},Matrix{Float64},Nothing})   # time: 0.0034816
    Base.precompile(Tuple{Type{BootTestResult{Float64}},Float64,String,Float64,Float64,Int64,Int64,Int64,Int64,Float64,NamedTuple{(:X, :p), Tuple{Tuple{Vector{Float64}}, Vector{Float64}}},NamedTuple{(:X, :p), Tuple{Vector{Float64}, Float64}},Matrix{Float64},Matrix{Float64},Vector{Float64},Matrix{Float64},Nothing})   # time: 0.003184
    Base.precompile(Tuple{typeof(coldotplus_nonturbo!),SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false},Matrix{Float64},Matrix{Float64}})   # time: 0.0031641
    Base.precompile(Tuple{typeof(coldotplus_nonturbo!),SubArray{Float32, 2, Matrix{Float32}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false},Matrix{Float32},Matrix{Float32}})   # time: 0.0028257
    Base.precompile(Tuple{Type{BootTestResult{Float32}},Float32,String,Float32,Float32,Int64,Int64,Int64,Int64,Float32,NamedTuple{(:X, :p), Tuple{Tuple{Vector{Float32}}, Vector{Float32}}},NamedTuple{(:X, :p), Tuple{Vector{Float32}, Float32}},Matrix{Float32},Matrix{Float32},Vector{Float32},Matrix{Float32},Nothing})   # time: 0.0028149
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Array{Float64, 3},Array{Float64, 3},Vector{UnitRange{Int64}}})   # time: 0.002783
    Base.precompile(Tuple{Type{BootTestResult{Float32}},Float32,String,Float32,Float32,Int64,Int64,Int64,Int64,Float32,NamedTuple{(:X, :p), Tuple{Tuple{Vector{Float32}, Vector{Float32}}, Vector{Float32}}},NamedTuple{(:X, :p), Tuple{Vector{Float32}, Float32}},Matrix{Float32},Matrix{Float32},Vector{Float32},Matrix{Float32},Nothing})   # time: 0.0027813
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Vector{Float32},Vector{Float32},Vector{UnitRange{Int64}}})   # time: 0.0027138
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Vector{Float64},Vector{Float64},Vector{UnitRange{Int64}}})   # time: 0.0026729
    Base.precompile(Tuple{typeof(panelsum_nonturbo!),Array{Float32, 3},Array{Float32, 3},Vector{UnitRange{Int64}}})   # time: 0.0026243
    Base.precompile(Tuple{typeof(coldotplus_nonturbo!),Matrix{Float64},Matrix{Float64},Matrix{Float64}})   # time: 0.0018618
    Base.precompile(Tuple{typeof(coldotplus_nonturbo!),Matrix{Float32},Matrix{Float32},Matrix{Float32}})   # time: 0.0017882
    Base.precompile(Tuple{typeof(matconvert),DataType,Vector{Int8}})   # time: 0.0016576
    Base.precompile(Tuple{Type{BootTestResult{Float64}},Float64,String,Float64,Float64,Int64,Int64,Int64,Int64,Float64,Nothing,Nothing,Nothing,Matrix{Float64},Vector{Float64},Matrix{Float64},Nothing})   # time: 0.0014966
    Base.precompile(Tuple{typeof(matconvert),DataType,Vector{Float32}})   # time: 0.0013353
    Base.precompile(Tuple{Type{BootTestResult{Float32}},Float32,String,Float32,Float32,Int64,Int64,Int64,Int64,Float32,Nothing,NamedTuple{(:X, :p), Tuple{Vector{Float32}, Float32}},Matrix{Float32},Matrix{Float32},Vector{Float32},Matrix{Float32},Nothing})   # time: 0.0013019
    Base.precompile(Tuple{Type{BootTestResult{Float64}},Float64,String,Float64,Float64,Int64,Int64,Int64,Int64,Float64,Nothing,Nothing,Nothing,Matrix{Float64},Vector{Float64},Matrix{Float64},SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}})   # time: 0.0012976
    Base.precompile(Tuple{Type{BootTestResult{Float32}},Float32,String,Float32,Float32,Int64,Int64,Int64,Int64,Float32,Nothing,Nothing,Nothing,Matrix{Float32},Vector{Float32},Matrix{Float32},Nothing})   # time: 0.001277
    Base.precompile(Tuple{Type{BootTestResult{Float32}},Float32,String,Float32,Float32,Int64,Int64,Int64,Int64,Float32,Nothing,Nothing,Nothing,Matrix{Float32},Vector{Float32},Matrix{Float32},SubArray{Float32, 2, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}})   # time: 0.00116
    Base.precompile(Tuple{Type{BootTestResult{Float64}},Float64,String,Float64,Float64,Int64,Int64,Int64,Int64,Float64,Nothing,NamedTuple{(:X, :p), Tuple{Vector{Float64}, Float64}},Matrix{Float64},Matrix{Float64},Vector{Float64},Matrix{Float64},Nothing})   # time: 0.0010974
    Base.precompile(Tuple{typeof(matconvert),Type{Int64},Matrix{Int64}})   # time: 0.0010696
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# merge interval tables
function mergetables(intervals, pars)
    # get column names
    bh, from, to = pars.holeid, pars.from, pars.to

    # read all interval tables
    interv = intervals isa Interval ? [intervals] : intervals
    tabs = [f.file isa String ? CSV.read(f.file, DataFrame) : f.file for f in interv]

    # rename main columns if necessary
    for i in 1:length(tabs)
        t = interv[i]
        rename!(tabs[i], t.holeid => bh, t.from => from, t.to => to)
    end

    # merge all tables and get unique from
    hcols = [bh,from,to]
    out   = vcat(map(x->select(x,hcols),tabs)..., cols=:union)
    out   = unique(vcat(select(out,[bh,from]),rename(select(out,[bh,to]), to => from), cols=:union))
    sort!(out, [bh, from])

    # add unique to and calculate length
    shift = collect(2:size(out,1))
    push!(shift,1) # check if works for every case
    out[!,to]   = out[shift,from]
    out[!,:CHK] = out[shift,bh]
    out = out[(out[!,to] .> out[!,from]) .& (out[!,bh] .== out[!,:CHK]),[bh,from,to]]
    out[!,:LENGTH] = out[!,to] - out[!,from]

    # add all tables values to the merged intervals
    cols = []
    for i in 1:length(tabs)
        push!(cols,setdiff(names(tabs[i]),string.(hcols)))
        tabs[i][!,"_LEN$i"] = tabs[i][!,to]-tabs[i][!,from]
        out = leftjoin(out,select(tabs[i], Not(to)),on=[bh,from],makeunique=true)
        select!(tabs[i], Not("_LEN$i"))
    end

    # leftjoin might affect order of the output after DataFrames 1.0. sort again
    sort!(out, [bh, from])

    # loop all merged intervals and columns
    for i in 1:size(out,1)
        for j in 1:length(cols)
            # check how many merged intervals compose one table interval
            ismissing(out[i,"_LEN$j"]) && continue
            adj = Dict(1 => 0, 2 => out[i,"_LEN$j"]-out[i,"LENGTH"])
            adj[2] ≈ 0 && continue
            while adj[2] > 0 && adj[2] ≉  0
                adj[1] += 1
                adj[2] -= out[i+adj[1],"LENGTH"]
            end

            # repeat table interval values to next i2 merged intervals
            i1, i2 = (i+1), (i+adj[1])
            for col in cols[j]
                crow = out[i,col]
                for k in i1:i2
                    out[k,col] = crow
                end
            end
        end
    end

    # delete auxiliary length variables
    tcols = [Symbol("_LEN$j") for j in 1:length(tabs)]
    select!(out, Not(tcols))
end

ext: jl


text:
let
    method = ME_Z_EXPECTED_G_BOUNDED

    # Feed jobs
    Ch = Channel(1) do ch
        cGLCs = Hd.val("cGLC")
        for (exp, cGLC)  in enumerate(cGLCs)
            put!(ch, (exp, cGLC))
        end
    end

    @threads for thid in 1:nthreads()
        for (exp, cGLC) in Ch
            # handle cache
            datfile = dat_file(string(DAT_FILE_PREFFIX, method); exp)
            if isfile(datfile)
                lock(WLOCK) do

                    @info("Cached loaded (skipping)",
                        exp, cGLC,datfile, thid
                    )
                    println()
                end
                continue
            end

            ## -------------------------------------------------------------------
            # setup
            model = load_model(exp)
            objidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
            M, N = size(model)
            exp_growth = Hd.val("D", exp)

            lock(WLOCK) do
                @info("Doing $(method)", exp, cGLC, thid)
                println()
            end

            ## -------------------------------------------------------------------
            # gradien descent
            epouts = Dict()
            x0 = 0.0
            @assert x0 == 0.0 # Must be zero for ME_Z_OPEN_G_OPEN
            x1 = 10.0
            maxΔ = 5e3
            th = 1e-3
            epout_seed = nothing
            target = exp_growth
            beta_vec = zeros(size(model, 2))

            upfrec_time = 10 # secunds
            last_uptime = time()
            gd_it = 1

            ## -------------------------------------------------------------------
            function upfun(beta)

                if haskey(epouts, beta)
                    epout = epouts[beta]
                else
                    beta_vec[objidx] = beta
                    epout = ChEP.maxent_ep(model;
                        beta_vec,
                        alpha = Inf,
                        damp = 0.9,
                        epsconv = 1e-4,
                        verbose = false,
                        solution = epout_seed,
                        maxiter = 5000
                    )
                    epouts[beta] = epout
                end
                epout_seed = epout
                ep_growth = ChU.av(epout)[objidx]

                update = gd_it == 1 || abs(last_uptime - time()) > upfrec_time ||
                    epout.status != ChEP.CONVERGED_STATUS
                update && lock(WLOCK) do
                    diff = abs.(exp_growth - ep_growth)
                    @info(
                        "Grad descent... ",
                        exp, gd_it,
                        epout.status, epout.iter,
                        ep_growth, exp_growth, diff,
                        beta,
                        thid
                    ); println()
                    gd_it += 1
                    last_uptime = time()
                end
                return ep_growth
            end

            ## -------------------------------------------------------------------
            # maxent
            expβ = UJL.grad_desc(upfun; x0, x1, th, maxΔ,
                target, maxiters = 2000, verbose = false
            )

            ## -------------------------------------------------------------------
            lock(WLOCK) do
                # Storing
                dat = Dict()
                dat[:exp_beta] = expβ
                dat[:epouts] = epouts
                dat[:exp_beta] = maximum(keys(epouts))
                dat[:model] = model |> ChU.compressed_model

                # caching
                serialize(datfile, dat)

                ep_growth = ChU.av(epouts[expβ])[objidx]
                diff = abs.(exp_growth - ep_growth)
                @info("Finished ",
                    exp, expβ,
                    length(epouts),
                    ep_growth, exp_growth, diff,
                    thid
                ); println()
            end

        end # for (exp, cGLC) in Ch
    end # for thid in 1:nthreads()
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ManagedDiskParameters <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    storageAccountType::Any # spec type: Union{ Nothing, StorageAccountType } # spec name: storageAccountType

    function ManagedDiskParameters(;id=nothing, storageAccountType=nothing)
        o = new()
        validate_property(ManagedDiskParameters, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ManagedDiskParameters, Symbol("storageAccountType"), storageAccountType)
        setfield!(o, Symbol("storageAccountType"), storageAccountType)
        o
    end
end # type ManagedDiskParameters

const _property_map_ManagedDiskParameters = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("storageAccountType")=>Symbol("storageAccountType"))
const _property_types_ManagedDiskParameters = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("storageAccountType")=>"StorageAccountType")
Base.propertynames(::Type{ ManagedDiskParameters }) = collect(keys(_property_map_ManagedDiskParameters))
Swagger.property_type(::Type{ ManagedDiskParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ManagedDiskParameters[name]))}
Swagger.field_name(::Type{ ManagedDiskParameters }, property_name::Symbol) =  _property_map_ManagedDiskParameters[property_name]

function check_required(o::ManagedDiskParameters)
    true
end

function validate_property(::Type{ ManagedDiskParameters }, name::Symbol, val)
end


text:
module EbitODEMessages
  using Compat
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("ebit-ode-messages_pb.jl")
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletaction

"""
    ''_leafletaction(;kwargs...)

A LeafletAction component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `mapActionConfig` (optional): . mapActionConfig has the following type: lists containing elements 'type', 'center', 'zoom', 'zoomInOffset', 'zoomOutOffset', 'bounds'.
Those elements have the following types:
  - `type` (a value equal to: 'set-zoom', 'zoom-in', 'zoom-out', 'set-view', 'fly-to', 'fly-to-bounds'; optional)
  - `center` (optional): . center has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
  - `zoom` (Real; optional)
  - `zoomInOffset` (Real; optional)
  - `zoomOutOffset` (Real; optional)
  - `bounds` (optional): . bounds has the following type: lists containing elements 'minx', 'miny', 'maxx', 'maxy'.
Those elements have the following types:
  - `minx` (Real; optional)
  - `miny` (Real; optional)
  - `maxx` (Real; optional)
  - `maxy` (Real; optional)
- `style` (Dict; optional)
"""
function ''_leafletaction(; kwargs...)
        available_props = Symbol[:id, :className, :loading_state, :mapActionConfig, :style]
        wild_props = Symbol[]
        return Component("''_leafletaction", "LeafletAction", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
# Pendulum
model,obj = TrajectoryOptimization.Dynamics.pendulum
n,m = model.n, model.m

u_bound = 5.
Q_min_time = Array(1e-3*Diagonal(I,n))
R_min_time = Array(1e-3*Diagonal(I,m))
Qf_min_time = Array(Diagonal(I,n)*0.0)
tf = obj.tf
x0 = obj.x0
xf = obj.xf

obj_min_time = TrajectoryOptimization.ConstrainedObjective(LQRObjective(Q_min_time, R_min_time, Qf_min_time, tf, x0, xf), tf=:min, u_min=-u_bound, u_max=u_bound)

solver_uncon = TrajectoryOptimization.Solver(model,obj,integration=:rk4,N=31)
solver_min = TrajectoryOptimization.Solver(model,obj_min_time,integration=:rk4,N=31)

solver_min.opts.verbose = false
solver_min.opts.max_dt = 0.15
solver_min.opts.min_dt = 1e-3
solver_min.opts.constraint_tolerance = 0.001 # 0.005
solver_min.opts.R_minimum_time = 15.0 #15.0 #13.5 # 12.0
solver_min.opts.constraint_decrease_ratio = .25
solver_min.opts.penalty_scaling = 2.0
solver_min.opts.outer_loop_update_type = :default
solver_min.opts.iterations = 1000
solver_min.opts.iterations_outerloop = 50 # 20

U = ones(m,solver_min.N-1)
results_uncon,stats_uncon = TrajectoryOptimization.solve(solver_uncon,U)
results_min,stats_min = TrajectoryOptimization.solve(solver_min,U)
# plot(TrajectoryOptimization.to_array(results_min.X)[1:2,:]')
# plot(TrajectoryOptimization.to_array(results_min.U)[1:2,:]')

T_uncon = TrajectoryOptimization.total_time(solver_uncon,results_uncon)
T_min_time = TrajectoryOptimization.total_time(solver_min,results_min)

@test T_min_time < 0.5*T_uncon
@test T_min_time < 1.5
@test norm(results_min.X[end][1:n] - obj.xf) < 1e-3
@test TrajectoryOptimization.max_violation(results_min) < solver_min.opts.constraint_tolerance

#######################
## Box Parallel Park ##
######################

# Solver Options
dt = 0.01
integration = :rk4

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-4
opts.cost_tolerance_intermediate = 1e-3

# Set up model, objective, and solver
N = 51
model, = TrajectoryOptimization.Dynamics.dubinscar
n, m = model.n,model.m

x0 = [0.0;0.0;0.]
xf = [0.0;1.0;0.]
tf =  3.
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]

u_max = 2.0
u_min = -2.0
obj_con_box = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max,u_max=u_max,u_min=u_min)

solver_uncon  = Solver(model, obj, integration=integration, N=N, opts=opts)
solver_con_box = Solver(model, obj_con_box, integration=integration, N=N, opts=opts)

U0 = ones(solver_con_box.model.m,solver_con_box.N-1)
# X0 = line_trajectory(solver_con_box)

# results_uncon, stats_uncon = TrajectoryOptimization.solve(solver_uncon,U0)
results_con_box, stats_con_box = TrajectoryOptimization.solve(solver_con_box,U0)

obj_mintime = update_objective(obj_con_box,tf=:min, u_min=u_min, u_max=u_max)

opts.max_dt = 0.2
opts.min_dt = 1e-3
opts.minimum_time_dt_estimate = tf/(N-1)
opts.constraint_tolerance = 0.001 # 0.005
opts.R_minimum_time = .05 #15.0 #13.5 # 12.0
opts.constraint_decrease_ratio = .25
opts.penalty_scaling = 10.0
opts.outer_loop_update_type = :individual
opts.iterations = 1000
opts.iterations_outerloop = 30 # 20

solver_mintime = Solver(model, obj_mintime, integration=integration, N=N, opts=opts)
solver_mintime.state.penalty_only = false
results_mintime, stats_mintime = solve(solver_mintime,to_array(results_con_box.U))
solver_mintime
T = TrajectoryOptimization.total_time(solver_con_box,results_con_box)
T_min = TrajectoryOptimization.total_time(solver_mintime,results_mintime)
# plot(to_array(results_mintime.U)[1:2,1:solver_mintime.N-1]',labels="")
# plot(to_array(results_mintime.X)[1:3,:]',labels="")
# plot(to_array(results_mintime.X)[1,:],to_array(results_mintime.X)[2,:],width=2,color=:blue,label="Minimum Time")

@test max_violation(results_mintime) < solver_mintime.opts.constraint_tolerance
@test T_min < 0.75*T
@test norm(results_mintime.X[end][1:n] - obj_con_box.xf) < 1e-3
@test T_min < 1.75


text:
#This script illustate how to set up constraints and project a 3D model onto an intersection
# with PARSDMM in serial, parallel or multilevel (serial or parallel)
# Bas Peters, 2017

#THE DATA FOR THIS SCRIPT IS CURRENTLY NOT AVAILABLE ONLINE

using Distributed
using LinearAlgebra
@everywhere using SetIntersectionProjection
using HDF5
ENV["MPLBACKEND"]="qt5agg"
using PyPlot

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

# Load velocity model
#get model at:  ftp://slim.eos.ubc.ca/data/SoftwareRelease/WaveformInversion.jl/3DFWI/overthrust_3D_true_model.h5
#run('wget ftp://slim.eos.ubc.ca/data/SoftwareRelease/WaveformInversion.jl/3DFWI/overthrust_3D_true_model.h5')
n,d,o,m = h5open("overthrust_3D_true_model.h5","r") do file
	read(file, "n", "d", "o", "m")
end
m .= 1000.0 ./ sqrt.(m);
m = m[50:200,50:200,:];

comp_grid = compgrid( (d[1], d[2], d[3]),( size(m,1), size(m,2), size(m,3) ) )
m=vec(m);

#PARSDMM options:
options=PARSDMM_options()
options.FL=Float64
#options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 500

set_zero_subnormals(true)
BLAS.set_num_threads(4)

#select working precision
if options.FL==Float64
  TF = Float64
  TI = Int64
elseif options.FL==Float32
  TF = Float32
  TI = Int32
end

#convert model and computational grid parameters
comp_grid = compgrid( ( convert(TF,comp_grid.d[1]),convert(TF,comp_grid.d[2]),convert(TF,comp_grid.d[3]) ),( convert(TI,comp_grid.n[1]),convert(TI,comp_grid.n[2]),convert(TI,comp_grid.n[3]) ) )
m         = convert(Vector{TF},m)

#define axis limits and colorbar limits
xmax = comp_grid.d[1]*comp_grid.n[1]
ymax = comp_grid.d[2]*comp_grid.n[2]
zmax = comp_grid.d[3]*comp_grid.n[3]
vmi=minimum(m)
vma=maximum(m)

#constraints
constraint = Vector{SetIntersectionProjection.set_definitions}() #initialize

#bounds:
m_min     = minimum(m)
m_max     = maximum(m)-200.0;
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#slope constraints (vertical)
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

options.parallel       = false
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

println("")
println("PARSDMM serial (bounds and bounds on D_z):")
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

#plot PARSDMM logs
figure();
subplot(3, 3, 3);semilogy(log_PARSDMM.r_pri)          ;title("r primal")
subplot(3, 3, 4);semilogy(log_PARSDMM.r_dual)         ;title("r dual")
subplot(3, 3, 1);semilogy(log_PARSDMM.obj)            ;title(L"$ \frac{1}{2} || \mathbf{m}-\mathbf{x} ||_2^2 $")
subplot(3, 3, 2);semilogy(log_PARSDMM.set_feasibility);title("TD feasibility violation")
subplot(3, 3, 5);plot(log_PARSDMM.cg_it)              ;title("nr. of CG iterations")
subplot(3, 3, 6);semilogy(log_PARSDMM.cg_relres)      ;title("CG rel. res.")
subplot(3, 3, 7);semilogy(log_PARSDMM.rho)            ;title("rho")
subplot(3, 3, 8);plot(log_PARSDMM.gamma)              ;title("gamma")
subplot(3, 3, 9);semilogy(log_PARSDMM.evol_x)         ;title("x evolution")

#plot
m_plot = reshape(m,comp_grid.n)
figure();
subplot(3,1,1);imshow(m_plot[:,:,Int64(round(comp_grid.n[3]/2))],cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, ymax, 0]); title("model to project x-y slice")
subplot(3,1,2);imshow(permutedims(m_plot[:,Int64(round(comp_grid.n[2]/2)),:],[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project x-z slice")
subplot(3,1,3);imshow(permutedims(m_plot[Int64(round(comp_grid.n[1]/2)),:,:],[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  ymax, zmax, 0]); title("model to project y-z slice")

x_plot = reshape(x,comp_grid.n)
figure();
subplot(3,1,1);imshow(x_plot[:,:,Int64(round(comp_grid.n[3]/2))],cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, ymax, 0]); title("projected model x-y slice")
subplot(3,1,2);imshow(permutedims(x_plot[:,Int64(round(comp_grid.n[2]/2)),:],[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("projected model x-z slice")
subplot(3,1,3);imshow(permutedims(x_plot[Int64(round(comp_grid.n[1]/2)),:,:],[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  ymax, zmax, 0]); title("projected model y-z slice")

#use multilevel-serial (2-levels)
#2 levels, the gird point spacing at level 2 is 3X that of the original (level 1) grid
options.parallel  = false
n_levels          = 2
coarsening_factor = 3

#set up all required quantities for each level
(TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels) = setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

println("")
println("PARSDMM multilevel-serial (bounds and bounds on D_z):")
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options)
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options)
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options)

#parallel single level
println("")
println("PARSDMM parallel (bounds and bounds on D_z):")
options.parallel       = true
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)

#use multilevel-parallel (2-levels)
options.parallel  = true
#2 levels, the grid point spacing at level 2 is 3X that of the original (level 1) grid
n_levels          = 2
coarsening_factor = 3

#set up all required quantities for each level
(TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels) = setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

println("PARSDMM multilevel-parallel (bounds and bounds on D_z):")
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);


text:
using Random
using LinearAlgebra
const ContactImplicitMPC = Main

include(joinpath(module_dir(), "src", "solver", "interior_point.jl"))
include(joinpath(module_dir(), "src", "solver", "mehrotra.jl"))

################################################################################
# Benchmark Utils
################################################################################

function get_benchmark_initialization(ref_traj::ContactTraj{T, nq},
		model::ContactModel, env::Environment, t::Int;
		z_dis = 1e-2, θ_dis = 1e-2) where {T,nq}
	Random.seed!(10)
	z = zeros(length(ref_traj.z[t]))
	q = deepcopy(ref_traj.θ[t][nq+1:2nq])
	z_initialize!(z, model, env, q)

	θ = deepcopy(ref_traj.θ[t])
	θ[1:nq] += θ_dis * ones(nq)
	return z, θ
end

function get_interior_point_solver(s::Simulation, space::Space, z, θ;
		r_tol = 1e-8, κ_tol = 1e-8, linear = false)
	model = s.model
	env = s.env

	lin = LinearizedStep(s, z, θ, 0.0) # /!| here we do not use κ from ref_traj
	if linear
		ip = interior_point(z, θ,
			s = space,
			oss = OptimizationSpace13(model, env),
			r!  = r!,
			rz! = rz!,
			rθ! = rθ!,
			r  = RLin(s, lin.z, lin.θ, lin.r, lin.rz, lin.rθ),
			rz = RZLin(s, lin.rz),
			rθ = RθLin(s, lin.rθ),
			opts = InteriorPointOptions(
			max_iter = 100,
			max_ls = 1;
			verbose = true,
			r_tol = r_tol,
			κ_tol = κ_tol,
			solver = :empty_solver,
			))
	else
		ip = interior_point(z, θ,
			s = space,
			oss = OptimizationSpace13(model, env),
			r! = s.res.r!,
			rz! = s.res.rz!,
			rθ! = s.res.rθ!,
			rz = s.rz,
			rθ = s.rθ,
			opts = InteriorPointOptions(
				max_iter = 100,
				r_tol = r_tol,
				κ_tol = κ_tol,
			))
	end
	return ip
end

################################################################################
# Benchmark Structures & Methods
################################################################################

# Benchmark
@with_kw mutable struct BenchmarkOptions{T}
	r_tol::T = 1.0e-8            # residual tolerance
	κ_tol::T = 1.0e-8            # bilinear tolerance
	z_dis::T = 1.0e-4            # disturbance on initial z
    θ_dis::T = 1.0e-4            # disturbance on initial θ
    linear::Bool = false         # whether to use or not the linearized residuals
end

mutable struct BenchmarkStatistics{T}
	failure::AbstractVector{Bool}
	iter::AbstractVector{Int}
	r_vio::AbstractVector{T}
	κ_vio::AbstractVector{T}
	cond_schur::AbstractVector{T}
end

function BenchmarkStatistics()
	failure = trues(0)
	iter = zeros(Int, 0)
	r_vio = zeros(0)
	κ_vio = zeros(0)
	cond_schur = zeros(0)
	return BenchmarkStatistics(failure, iter, r_vio, κ_vio, cond_schur)
end

function record!(stats::BenchmarkStatistics, s, i, r, κ, c)
	push!(stats.failure, s)
	push!(stats.iter, i)
	push!(stats.r_vio, r)
	push!(stats.κ_vio, κ)
	push!(stats.cond_schur, c)
	return nothing
end

function benchmark_interior_point!(stats::BenchmarkStatistics, s::Simulation,
	ref_traj::ContactTraj, space::Space; opts::BenchmarkOptions)

	nz = num_var(s.model, s.env)
	for t = 1:ref_traj.H
		z, θ = get_benchmark_initialization(ref_traj, s.model, s.env, t,
			z_dis = opts.z_dis,
			θ_dis = opts.θ_dis)
		ip = get_interior_point_solver(s, space, deepcopy(ref_traj.z[t]), deepcopy(ref_traj.θ[t]),
			r_tol = opts.r_tol,
			κ_tol = opts.κ_tol,
			linear = opts.linear)
		interior_point_solve!(ip, z, θ)
		rv = residual_violation(ip, ip.r)
		κv = bilinear_violation(ip, ip.r)
		fl = !((rv < opts.r_tol) && (κv < opts.κ_tol))
		record!(stats, fl, ip.iterations, rv, κv, NaN)
	end
	return nothing
end

function benchmark!(stats::BenchmarkStatistics, s::Simulation,
		ref_traj::ContactTraj, space::Space; opts::BenchmarkOptions, solver::Symbol)
	if solver == :interior_point
		benchmark_interior_point!(stats, s, ref_traj, space; opts = opts)
	else
		@warn "Unknown solver"
	end
	return nothing
end

function benchmark(s::AbstractVector{Simulation},
		ref_traj::AbstractVector{<:ContactTraj}, space::AbstractVector{<:Space}, dist::AbstractVector{<:Real};
		opts = BenchmarkOptions(), solver::Symbol)

	stats = Vector{BenchmarkStatistics}()
	nsim = length(s)
	ndis = length(dist)
	@assert length(ref_traj) == nsim

	for i = 1:ndis
		sta = BenchmarkStatistics()
		opts.z_dis = dist[i]
		opts.θ_dis = dist[i]
		for j = 1:nsim
			benchmark!(sta, s[j], ref_traj[j], space[j]; opts = opts, solver = solver)
		end
		push!(stats, sta)
	end
	return stats
end

################################################################################
# Loading Problems
################################################################################

s_quadruped = ContactImplicitMPC.get_simulation("quadruped", "flat_2D_lc", "flat")
s_flamingo = ContactImplicitMPC.get_simulation("flamingo", "flat_2D_lc", "flat")
s_hopper = ContactImplicitMPC.get_simulation("hopper_2D", "flat_2D_lc", "flat")
s_particle_2D = ContactImplicitMPC.get_simulation("particle_2D", "flat_2D_nc", "flat_nc")
s_particle = ContactImplicitMPC.get_simulation("particle", "flat_3D_nc", "flat_nc")
s_box_quat_lc = get_simulation("box", "flat_3D_lc", "flat_lc")
s_box_quat_nc = get_simulation("box", "flat_3D_nc", "flat_nc")

ref_traj_quadruped = deepcopy(ContactImplicitMPC.get_trajectory(s_quadruped.model, s_quadruped.env,
    joinpath(ContactImplicitMPC.module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt))
ref_traj_flamingo = deepcopy(ContactImplicitMPC.get_trajectory(s_flamingo.model, s_flamingo.env,
    joinpath(ContactImplicitMPC.module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))
ref_traj_hopper = deepcopy(ContactImplicitMPC.get_trajectory(s_hopper.model, s_hopper.env,
    joinpath(ContactImplicitMPC.module_dir(), "src/dynamics/hopper_2D/gaits/gait_forward.jld2"),
    load_type = :joint_traj))
ref_traj_particle_2D = deepcopy(get_trajectory(s_particle_2D.model, s_particle_2D.env,
	joinpath(module_dir(), "src/dynamics/particle_2D/gaits/gait_NC.jld2"),
	load_type = :joint_traj))
ref_traj_particle = deepcopy(get_trajectory(s_particle.model, s_particle.env,
	joinpath(module_dir(), "src/dynamics/particle/gaits/gait_NC.jld2"),
	load_type = :joint_traj))
ref_traj_box_quat_lc = deepcopy(get_trajectory(s_box_quat_lc.model, s_box_quat_lc.env,
	joinpath(module_dir(), "src/dynamics/box/gaits/gait_LC.jld2"),
	load_type = :joint_traj))
ref_traj_box_quat_nc = deepcopy(get_trajectory(s_box_quat_nc.model, s_box_quat_nc.env,
	joinpath(module_dir(), "src/dynamics/box/gaits/gait_NC.jld2"),
	load_type = :joint_traj))

space_quadruped = Euclidean(num_var(s_quadruped.model, s_quadruped.env))
space_flamingo = Euclidean(num_var(s_flamingo.model, s_flamingo.env))
space_hopper = Euclidean(num_var(s_hopper.model, s_hopper.env))
space_particle_2D = Euclidean(num_var(s_particle_2D.model, s_particle_2D.env))
space_particle = Euclidean(num_var(s_particle.model, s_particle.env))
space_box_quat_lc = rn_quaternion_space(
	num_var(s_box_quat_lc.model, s_box_quat_lc.env) - 1,
	x -> Gz_func(s_box_quat_lc.model, s_box_quat_lc.env, x),
	collect([(1:3)..., (8:num_var(s_box_quat_lc.model, s_box_quat_lc.env))...]),
	collect([(1:3)..., (7:num_var(s_box_quat_lc.model, s_box_quat_lc.env)-1)...]),
	[collect((4:7))],
	[collect((4:6))])
space_box_quat_nc = rn_quaternion_space(
	num_var(s_box_quat_nc.model, s_box_quat_nc.env) - 1,
	x -> Gz_func(s_box_quat_nc.model, s_box_quat_nc.env, x),
	collect([(1:3)..., (8:num_var(s_box_quat_nc.model, s_box_quat_nc.env))...]),
	collect([(1:3)..., (7:num_var(s_box_quat_nc.model, s_box_quat_nc.env)-1)...]),
	[collect((4:7))],
	[collect((4:6))])

################################################################################
# Running Benchmark
################################################################################

s = [
	s_quadruped, s_flamingo, s_hopper,
	# s_particle_2D, s_particle,
	# s_box_quat_lc,
	# s_box_quat_nc
	]
ref_traj = [
	ref_traj_quadruped, ref_traj_flamingo, ref_traj_hopper,
	# ref_traj_particle_2D, ref_traj_particle,
	# ref_traj_box_quat_lc,
	# ref_traj_box_quat_nc
	]
space = [
	space_quadruped, space_flamingo, space_hopper,
	# space_particle_2D, space_particle,
	# space_box_quat_lc,
	# space_box_quat_nc
	]
# s = [s_particle_2D, s_particle]
# ref_traj = [ref_traj_particle_2D, ref_traj_particle]
# s = [s_quadruped, s_flamingo, s_hopper]
# ref_traj = [ref_traj_quadruped, ref_traj_flamingo, ref_traj_hopper]

dist = [1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
tol = 1e-7
opts_nonlin = BenchmarkOptions(r_tol = tol, κ_tol = tol, linear = false)
opts_lin = BenchmarkOptions(r_tol = tol, κ_tol = tol, linear = true)
stats = Dict{Symbol, AbstractVector{BenchmarkStatistics}}()
keys = [
	:interior_point_lin,
	:interior_point_nonlin,
	]

stats[:interior_point_nonlin] = benchmark(
	s,
	deepcopy.(ref_traj),
	space,
	dist,
	opts = opts_nonlin,
	solver = :interior_point,
	)

stats[:interior_point_lin] = benchmark(
	s,
	deepcopy.(ref_traj),
	space,
	dist,
	opts = opts_lin,
	solver = :interior_point,
	)

function display_statistics!(plt, plt_ind, stats::AbstractVector{BenchmarkStatistics}, dist;
		field::Symbol = :iter, solver::String = "solver")
	data = [mean(getfield(s, field)) for s in stats]
	linestyle = occursin("nonlin", solver) ? :solid : :dot
	linewidth = occursin("nonlin", solver) ? 3.0 : 5.0
	marker    = occursin("nonlin", solver) ? :circ : :utriangle
	if occursin("interior_point_base", solver)
		color = :red
	elseif occursin("interior_point", solver)
		color = :orange
	elseif occursin("mehrotra", solver)
		color = :black
	end

	plot!(plt[plt_ind...], dist, data,
		xlabel = "disturbances",
		ylabel = String(field),
		label = false,
		xlims = [minimum(dist), maximum(dist)*100],
		ylims = [0, Inf],
		xaxis = :log,
		legend = :right,
		linestyle = linestyle,
		color = color,
		linewidth = linewidth)
	scatter!(plt[plt_ind...], dist, data,
		label = solver,
		marker = marker,
		markercolor = color,
		markeralpha = 0.6,
		markerstrokecolor = :black,
		markerstrokewidth = 2.0,
		color = color,
		markersize = 8.0,
		)
	display(plt)
	return nothing
end

################################################################################
# Benchmark Results
################################################################################

plt = plot(size = (800, 600), layout = (3, 1), fmt = :svg)
for k in keys
	display_statistics!(plt, [1,1], stats[k], dist, field = :iter, solver = String(k))
end
for k in keys
	display_statistics!(plt, [2,1], stats[k], dist, field = :failure,  solver = String(k))
end
for k in keys
	display_statistics!(plt, [3,1], stats[k], dist, field = :cond_schur,  solver = String(k))
end


text:
using MatrixCompletion

const FLAG_SIMULATION_GAMMA_VARY_RANK_SMALL                 = false
const FLAG_SIMULATION_GAMMA_VARY_RANK_MEDIUM                = true
const FLAG_SIMULATION_GAMMA_VARY_RANK_LARGE                 = false

const FLAG_SIMULATION_GAMMA_VARY_MISSING_PERCENTAGE_SMALL   = false
const FLAG_SIMULATION_GAMMA_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const FLAG_SIMULATION_GAMMA_VARY_MISSING_PERCENTAGE_LARGE   = false

FLAG_SIMULATION_GAMMA_VARY_RANK_SMALL ?
  include("simulation_gamma_vary_rank_small.jl") : nothing

FLAG_SIMULATION_GAMMA_VARY_RANK_MEDIUM ?
  include("simulation_gamma_vary_rank_medium.jl") : nothing

# @testset "$(format("ADMM Algorithm: Small Input Simulation [Gamma 400 x 400]"))" begin
#   if SIMULATION_STATUS_GAMMA_VARY_MISSING_PERCENTAGE_LARGE == false
#     @info("Running Simulation: [Gamma vary rank][SMALL]")
#     let
#       Random.seed!(65536)
#       ROW = 400
#       COL = 400
#       for input_rank in 1:2:400
#         @printf("small case: rank = %d\n", input_rank)
#         dist = Gamma()
#         timer = TimerOutput()
#         Base.Filesystem.mkpath("./test_result/gamma/small(400x400)")
#         io = open("./test_result/gamma/small(400x400)/rank"*string(input_rank)*".log", "w")
#         # io = Base.stdout
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Gamma(10, 0.5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Gamma
#         @timeit timer  "Gamma(400x400)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = nothing,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end

#         log_simulation_result(Gamma(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#   else
#     @info("Already Completed Simualtion[Gamma vary rank][SMALL]")
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input Simulation [Gamma 2000 x 2000]"))" begin
#   if SIMULATION_STATUS_GAMMA_VARY_RANK_MEDIUM == false
#     @info("Running Simulation: [Gamma vary rank][MEDIUM]")
#     let
#       Random.seed!(65536)
#       ROW = 2000
#       COL = 2000
#       for input_rank in collect(400:10:500)
#         @printf("medium case: rank = %d\n", input_rank)
#         dist = Gamma()
#         timer = TimerOutput()
#         Base.Filesystem.mkpath("./test_result/gamma/medium(2000x2000)")
#         io = open("./test_result/gamma/medium(2000x2000)/rank"*string(input_rank)*".log", "w")
#         # io = Base.stdout
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Gamma(10, 0.5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Gamma
#         @timeit timer  "Gamma(2000x2000)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = input_rank,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end

#         log_simulation_result(Gamma(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#     else
#       @info("Already Completed Simualtion[Gamma vary rank][MEDIUM]")
#     end
# end

  # @testset "$(format("ADMM Algorithm: Large Input Simulation [Gamma 4000 x 4000]"))" begin
  #   let
  #     Random.seed!(65536)
  #     ROW = 2000
  #     COL = 2000
  #     for input_rank in union(1, collect(10:20:1000))
  #       @printf("medium case: rank = %d\n", input_rank)
  #       dist = Gamma()
  #       timer = TimerOutput()
  #       Base.Filesystem.mkpath("./test_result/gamma/medium(2000x2000)")
  #       io = open("./test_result/gamma/medium(2000x2000)/rank"*string(input_rank)*".log", "w")
  #       # io = Base.stdout
  #       truth_matrix        = rand([(FixedRankMatrix(Distributions.Gamma(10, 0.5), rank = input_rank), ROW, COL)])
  #       sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
  #       input_matrix        = sample_model.draw(truth_matrix)
  #       manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
  #       manual_type_matrix .= :Gamma
  #       @timeit timer  "Gamma(2000x2000)" * "| rank="*string(input_rank) begin
  #         completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
  #                                                            maxiter         = 200,
  #                                                            ρ               = 0.3,
  #                                                            use_autodiff    = false,
  #                                                            gd_iter         = 3,
  #                                                            debug_mode      = false,
  #                                                            project_rank    = nothing,
  #                                                            io              = io,
  #                                                            type_assignment = manual_type_matrix)
  #       end

  #       log_simulation_result(Gamma(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
  #       show(io, MIME("text/plain"), timer)
  #     end
  #   end
  # end


text:
#!/usr/bin/env julia

using Dates

include("library.jl")

"""
    TemplateBody(m, x, y, z, v)

Create a template of a celestial body. By default, all arguments are `nothing`.

See also: [`Body`](@ref)
"""
Base.@kwdef mutable struct TemplateBody
	m::Union{Float64, Nothing} = nothing

	x::Union{Float64, Nothing} = nothing
	y::Union{Float64, Nothing} = nothing
	z::Union{Float64, Nothing} = nothing

	v::Vector{Union{Float64, Nothing}} = [nothing, nothing, nothing]
end

"""
    Body(tb::TemplateBody)

Create a Body object from a TemplateBody object. Will throw errors if any attribute is `nothing`.
"""
function Body(tb::TemplateBody)::Body
	if isnothing(tb.m); error("Mass must not be nothing"); end
	if isnothing(tb.x); error("x must not be nothing"); end
	if isnothing(tb.y); error("y must not be nothing"); end
	if isnothing(tb.z); error("z must not be nothing"); end
	if isnothing(tb.v[1]); error("Velocity x component must not be nothing"); end
	if isnothing(tb.v[2]); error("Velocity y component must not be nothing"); end
	if isnothing(tb.v[3]); error("Velocity z component must not be nothing"); end

	Body(tb.m, tb.x, tb.y, tb.z, MVector{3, Float64}(tb.v[1], tb.v[2], tb.v[3]))
end

"""
    parsenums(num::String, len::Int)

Parse the string `num` and return a list of all indices to target.

`len` is the length of the bodies list. If any number is larger than it,
then we throw an error. If any number is less than 1, we throw an error.

See also: [`parseargs`](@ref)
"""
function parsenums(num::String, len::Int)::Vector{Int}
	if num == "a"
		return [1:len;]

	# If it's not a range or group
	elseif !occursin("-", num) && !occursin(".", num)
		n = parse(Int, num)
		if n > len; error("You cannot select body $n; there are only $len bodies."); end
		if n < 1; error("You cannot select body number $n; the minimum is 1."); end
		return [n]

	# If it's not a range but it is a group
	elseif !occursin("-", num) && occursin(".", num)
		ns = [parse(Int, i) for i in split(num, ".")]

		if max(ns...) > len; error("You cannot select body $(max(ns...)); there are only $len bodies."); end
		if min(ns...) < 1; error("You cannot select body number $(min(ns...)); the minimum is 1."); end
		return ns

	# If it's a range but not a group
	elseif occursin("-", num) && !occursin(".", num)
		first = parse(Int, split(num, "-")[1])
		second = parse(Int, split(num, "-")[2])
		ns = [first:second;]

		if max(ns...) > len; error("You cannot select body $(max(ns...)); there are only $len bodies."); end
		if min(ns...) < 1; error("You cannot select body number $(min(ns...)); the minimum is 1."); end
		return ns

	# If it's both a range and a group
	else
		splitnum = split(num, ".")

		# We recur down to get the ranges of each part of splitnum
		# We don't need to check that all nums are <= len, because we do that check in the recurrance
		return collect(Iterators.flatten([parsenums(splitnum[i], len) for i in 1:length(splitnum)]))
	end
end

# When we recur, we're calling parsenums with SubString{String}, so we have to convert it to String
parsenums(a::SubString{String}, b::Int) = parsenums(string(a), b)

"""
    parseargs(progname, args)

Parse the arguments and call creategif(). Should be called with PROGRAM_FILE and ARGS.

See also: [`parsenum`](@ref)
"""
function parseargs(progname::String, args::Vector{String})
	# We need to put a space in front to split by " -" and avoid breaking negatives
	connectedargs = " " * join(args, " ")

	if occursin("--help", connectedargs) || occursin("-h", connectedargs) || length(args) < 4 # We need at least 4 args
		println("""
				Usage: $(Sys.iswindows() ? "julia " : "")$progname [--help] -n <number> -f <frames> [-t <seconds>] [-F <filename>] [--cube] [--initial-bounds] [--quiet] [attribute options]

				Options:
				  --help, -h         Display this help text.

				  -n <number>        The number of bodies.
				  -f <frames>        The number of frames to use in the GIF.
				  -t <seconds>       The time step in seconds (default 60).
				  -F <filename>      The name of the gif (will generate a unique name if not specified).

				  --cube             Give the plot cubic bounds.
				  --initial-bounds   Give the plot bounds to only include the initial positions.
				  --quiet            Suppress terminal output; only write command to log file.

				  -m n,m             Set the mass of body number n to be m kg.
				  -p n,x,y,z         Set the initial position of body number n to be (x, y, z).
				  -v n,x,y,z         Set the initial velocity of body number n to be [x, y, z].

				See README.md (read it at https://github.com/DoctorDalek1963/Julia-Gravity) for a full guide on using the options.
				""")
		return
	end

	# We split the args by " -", get rid of the "" in the list, and strip
	# the whitespace from the end of each element
	arglist = strip.(split(connectedargs, " -", keepempty=false))

	n::Union{Int, Nothing} = nothing
	frames::Union{Int, Nothing} = nothing
	Δt::Float64 = 60.0
	filename::Union{String, Nothing} = nothing
	cube = false
	initialbounds = false
	quiet = false

	# We loop through and get all the meta values before getting the body attributes
	# This allows us to make a fixed length StaticArray of TemplateBody objects
	for i in 1:length(arglist)
		arg = arglist[i]

		if startswith(arg, "n")
			if !isnothing(n) # If we're trying to redefine n
				error("-n may only be defined once")
			else
				n = parse(Int, split(arg, " ")[2])
				if n < 1; error("n must be at least 1, not $n."); end
			end

		elseif startswith(arg, "f")
			if !isnothing(frames) # If we're trying to redefine n
				error("-f may only be defined once")
			else
				frames = parse(Int, split(arg, " ")[2])
				if frames < 1; error("frames must be at least 1, not $frames."); end
			end

		elseif startswith(arg, "t")
			Δt = parse(Float64, split(arg, " ")[2])

		elseif startswith(arg, "F")
			filename = split(arg, " ")[2]

		# We need 1 "-" here because we split by " -" and these args are long form
		elseif arg == "-cube"
			cube = true

		elseif arg == "-initial-bounds"
			initialbounds = true

		elseif arg == "-quiet"
			quiet = true
		end
	end

	if isnothing(n); error("-n must be specified"); end
	if isnothing(frames); error("-f must be specified"); end

	templatebodies = MVector{n, TemplateBody}([TemplateBody() for _ in 1:n]...)

	for i in 1:length(arglist)
		arg = arglist[i]

		if startswith(arg, "m")
			datalist = split(split(arg, " ")[2], ",")
			if length(datalist) > 2; error("Mass only accepts one value"); end

			nums = parsenums(datalist[1], n)

			value = parse(Float64, datalist[2])
			if value < 0; error("Mass must not be negative"); end

			for i in 1:n
				# If this is a body that we want to edit
				if in(i, nums)
					templatebodies[i].m = value
				end
			end

		elseif startswith(arg, "p")
			data = split(arg, " ")[2]
			datalist = split(data, ",")

			# We use parsenums() just to check that the user only gave 1 number,
			# so that we can give them a better error message
			nums = parsenums(datalist[1], n)
			if length(nums) > 1
				error("-p only accepts single number selectors; \"$(datalist[1])\" is invalid")
			else
				num = nums[1]
			end

			if num > n; error("You cannot select body $num; there are only $n bodies."); end

			# If we don't have any special selectors and they're just positional
			if !occursin("x", data) && !occursin("y", data) && !occursin("z", data)
				# If we don't have the right number of values
				if length(datalist) != 2 && length(datalist) != 4
					error("-p only accepts 1 or 3 positional values, not $(length(datalist[2:end]))")
				end

				# If we've only got 1 number, use that for all
				if length(datalist[2:end]) == 1
					templatebodies[num].x = parse(Float64, datalist[2])
					templatebodies[num].y = parse(Float64, datalist[2])
					templatebodies[num].z = parse(Float64, datalist[2])
				else
					templatebodies[num].x = parse(Float64, datalist[2])
					templatebodies[num].y = parse(Float64, datalist[3])
					templatebodies[num].z = parse(Float64, datalist[4])
				end

			# If we've got explicit values
			else
				for i in 2:length(datalist)
					value = datalist[i]

					# If we get a value that doesn't start with x, y, or z, then we just ignore its
					# because we can't mix explicit coords with positional coords
					if startswith(value, "x"); templatebodies[num].x = parse(Float64, value[2:end]); end
					if startswith(value, "y"); templatebodies[num].y = parse(Float64, value[2:end]); end
					if startswith(value, "z"); templatebodies[num].z = parse(Float64, value[2:end]); end
				end
			end

		elseif startswith(arg, "v")
			data = split(arg, " ")[2]
			datalist = split(data, ",")

			nums = parsenums(datalist[1], n)

			# If we don't have any special selectors and they're just positional
			if !occursin("x", data) && !occursin("y", data) && !occursin("z", data)
				# If we don't have the right number of values
				if length(datalist) != 2 && length(datalist) != 4
					error("-v only accepts 1 or 3 positional values, not $(length(datalist[2:end]))")
				end

				# If we've only got 1 number, use that for all
				if length(datalist[2:end]) == 1
					values = [parse(Float64, datalist[2]) for _ in 1:3]
				else
					values = [parse(Float64, datalist[i]) for i in 2:4]
				end

				for i in 1:n
					if in(i, nums)
						templatebodies[i].v = values
					end
				end

			# If we've got explicit values
			else
				for i in 2:length(datalist)
					value = datalist[i]

					for i in 1:n
						if in(i, nums)
							# If we get a value that doesn't start with x, y, or z, then we just ignore its
							# because we can't mix explicit values with positional values
							if startswith(value, "x"); templatebodies[i].v[1] = parse(Float64, value[2:end]); end
							if startswith(value, "y"); templatebodies[i].v[2] = parse(Float64, value[2:end]); end
							if startswith(value, "z"); templatebodies[i].v[3] = parse(Float64, value[2:end]); end
						end
					end
				end
			end
		end
	end

	bodies = Body[]

	for i in 1:n
		tb = templatebodies[i]

		# Values in a similar order of magnitude to 1e22
		if isnothing(tb.m); tb.m = 1e22 * rand() * rand(1:100); end

		# Values from 0 to 50e6 with variation, positive and negative
		if isnothing(tb.x); tb.x = randn() * rand(0:50000000); end
		if isnothing(tb.y); tb.y = randn() * rand(0:50000000); end
		if isnothing(tb.z); tb.z = randn() * rand(0:50000000); end

		# Values from 0 to 500 with variation, positive and negative
		if isnothing(tb.v[1]); tb.v[1] = randn() * rand(0:500); end
		if isnothing(tb.v[2]); tb.v[2] = randn() * rand(0:500); end
		if isnothing(tb.v[3]); tb.v[3] = randn() * rand(0:500); end

		# We created a custom Body constuctor earlier to construct a
		# Body from a TemplateBody
		push!(bodies, Body(tb))
	end

	# If desired, we get the initial bounds and use those
	if initialbounds
		xs = [bodies[i].x for i in 1:length(bodies)]
		ys = [bodies[i].y for i in 1:length(bodies)]
		zs = [bodies[i].z for i in 1:length(bodies)]

		bounds = [
				(min(xs...), max(xs...)),
				(min(ys...), max(ys...)),
				(min(zs...), max(zs...))
				]
	else
		bounds = nothing
	end

	# Generate the command and log it
	command = ""
	command *= "$progname -n $n -f $frames -t $Δt"
	if !isnothing(filename); command *= " -F $filename"; end
	if cube; command *= " --cube"; end
	if initialbounds; command *= " --initial-bounds"; end

	for i in 1:n
		b = bodies[i]
		command *= " -m $i,$(b.m) -p $i,$(b.x),$(b.y),$(b.z) -v $i,$(b.v[1]),$(b.v[2]),$(b.v[3])"
	end

	# Write the command to the log file
	open("command_log.txt", "a") do file
		write(file, "[$(Dates.now())] $command\n\n")
	end

	if !quiet
		println("\nThese are the arguments needed to recreate this simulation:\n")
		println(command)
		println()
	end

	println("Simulating...")
	creategif(bodies, frames, Δt, cube, bounds, filename)
	println("Saved gif")
end

if abspath(PROGRAM_FILE) == @__FILE__
	parseargs(PROGRAM_FILE, ARGS)
end


text:
using JLD2, Plots, Printf, Oceananigans

ENV["GKSwstype"] = "100"

#name = "veron_and_melville_Nz128_β1.0e-04"
name = "veron_and_melville_Nz128_Ly0.3_β1.0e-05_unsteady_waves"

yz_file = jldopen(name * "_yz.jld2")
xz_file = jldopen(name * "_xz.jld2")

Nx = yz_file["grid/Nx"]
Ny = yz_file["grid/Ny"]
Nz = yz_file["grid/Nz"]

Lx = yz_file["grid/Lx"]
Ly = yz_file["grid/Ly"]
Lz = yz_file["grid/Lz"]

grid = RegularRectilinearGrid(size = (Nx, Ny, Nz), halo = (3, 3, 3),
                              x = (0, Lx),
                              y = (0, Ly),
                              z = (-Lz, 0),
                              topology = (Periodic, Bounded, Bounded))

xw, yw, zw = nodes((Center, Center, Face), grid)
xu, yu, zu = nodes((Face, Center, Center), grid)
xc, yc, zc = nodes((Center, Center, Center), grid)

iterations = parse.(Int, keys(yz_file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)
    wyz = yz_file["timeseries/w/$iter"][1, :, :]
    uyz = yz_file["timeseries/u/$iter"][1, :, :]
    uxz = xz_file["timeseries/u/$iter"][:, 1, :]
    cxz = xz_file["timeseries/c/$iter"][:, 1, :]

    t = yz_file["timeseries/t/$iter"]

    wmax = maximum(abs, wyz)
    wlim = wmax / 2
    wlevels = range(-wlim, stop=wlim, length=30)
    wlim < wmax && (wlevels = vcat([-wmax], wlevels, [wmax]))

    umax = maximum(abs, uxz)
    ulim = umax / 2
    ulevels = range(-ulim, stop=ulim, length=30)
    ulim < umax && (ulevels = vcat([-umax], ulevels, [umax]))

    cmax = maximum(abs, cxz)
    clim = cmax * 0.8
    clevels = range(0, stop=clim, length=32)
    clim < cmax && (clevels = vcat([-clevels[2]], clevels, [cmax]))

    u_title = @sprintf("u(x, z, t) (m s⁻¹) at t = %s ", prettytime(t))
    wyz_title = "w(y, z, t) (m s⁻¹)"
    uyz_title = "u(y, z, t) (m s⁻¹)"
    cxz_title = "c(x, z, t)"

    @show u_title

    uxzplot = contourf(xu, zu, uxz';
                       linewidth = 0,
                       aspectratio = :equal,
                       xlims = (0, grid.Lx),
                       ylims = (-grid.Lz, 0),
                       color = :balance,
                       levels = ulevels,
                       clims = (-ulim, ulim))

    cxzplot = contourf(xc, zc, cxz';
                       linewidth = 0,
                       aspectratio = :equal,
                       xlims = (0, grid.Lx),
                       ylims = (-grid.Lz, 0),
                       color = :thermal,
                       levels = clevels,
                       clims = (0, clim))

    wyzplot = contourf(yw, zw, wyz';
                       linewidth = 0,
                       aspectratio = :equal,
                       xlims = (0, grid.Ly),
                       ylims = (-grid.Lz, 0),
                       color = :balance,
                       levels = wlevels,
                       clims = (-wlim, wlim))

    uyzplot = contourf(yu, zu, uyz';
                       linewidth = 0,
                       aspectratio = :equal,
                       xlims = (0, grid.Ly),
                       ylims = (-grid.Lz, 0),
                       color = :balance,
                       levels = ulevels,
                       clims = (-ulim, ulim))

    layout = @layout [
        a
        b
        Plots.grid(1, 2)
    ]

    plot(uxzplot, cxzplot, wyzplot, uyzplot,
         layout = layout,
         size = (1200, 1200),
         title = [u_title cxz_title wyz_title uyz_title])
end

close(yz_file)
close(xz_file)

mp4(anim, name * ".mp4", fps=8)


text:
"""
The `Cursor` type is embeded into Branches of a TTree such that when
we need to read the content of a Branch, we don't need to go through
the Directory and find the TKey and then seek to where the Branch is.

!!! note
    The `io` inside a `Cursor` is in fact only a buffer, it is NOT
    a `io` that refers to the whole file's stream.
"""
struct Cursor
    start::Int64
    io::IO
    tkey
    refs::Dict{Int32, Any}
end

Base.position(c::Cursor) = position(c.io)

function unpack() end
packedsizeof(T::Type) = sum(sizeof.(fieldtypes(T)))

@inline readtype(io, ::Type{T}) where T<:Union{Integer, AbstractFloat} = ntoh(read(io, T))
@inline readtype(io, ::Type{T}) where T<:Bool = read(io, T)
@inline readtype(io, v::Type{T}) where T<:AbstractVector{UInt8} = read(io, length(v))

# Non-C strings in .root are preceeded by 1 or more bytes signifying the length
# of the string that follows.
function readtype(io, ::Type{T}) where T<:AbstractString
    length = readtype(io, UInt8)

    if length == 255
        # first byte 0xff is useless now
        # https://github.com/scikit-hep/uproot3/blob/54f5151fb7c686c3a161fbe44b9f299e482f346b/uproot3/source/cursor.py#L91
        length = readtype(io, UInt32)
    end

    T(read(io, length))
end

struct CString
    value::String
end

function readtype(io, ::Type{T}) where {T<:CString}
    out = Char[]
    char = read(io, Char)
    while char != '\0'
        push!(out, char)
        char = read(io, Char)
    end
    String(out)
end

macro io(data)
    struct_name = data.args[2]
    types = []
    for f in data.args[3].args
        isa(f, LineNumberNode) && continue
        isa(f, Symbol) && error("Untyped field")
        Meta.isexpr(f, :(::)) || error("")
        push!(types, f.args[2])
    end

    # TODO: Need to figure out how to deal with Strings, probably dynamically
    # create the sizes instead of defining it just based on sizeof.
    # struct_size = sum([sizeof(eval(t)) for t in types])

    quote
        $(esc(data))  # executing the code to create the actual struct
        # Base.sizeof(::Type{$(esc(struct_name))}) = $struct_size

        function $(@__MODULE__).unpack(io, ::Type{T}) where {T<:$(esc(struct_name))}
            $(esc(struct_name))($([:(readtype(io, $t)) for t in types]...))
        end

        nothing  # supress REPL output
    end
end

struct Preamble
    start::Int64
    cnt::Union{UInt32, Missing}
    version::UInt16
    type::Type
end

"""
Reads the preamble of an object.

The cursor will be put into the right place depending on the data.
"""
function Preamble(io, ::Type{T}) where {T}
    start = position(io)
    cnt = readtype(io, UInt32)
    version = readtype(io, UInt16)
    if Int64(cnt) & Const.kByteCountMask > 0
        cnt = Int64(cnt) & ~Const.kByteCountMask
        return Preamble(start, cnt + 4, version, T)
    else
        seek(io, start)
        version = readtype(io, UInt16)
        return Preamble(start, missing, version, T)
    end
end

function Preamble(io::Cursor, ::Type{T}) where {T}
    Preamble(io.io, T)
end

"""
    function skiptobj(io)

Skips a TOBject.
"""
function skiptobj(io)
    version = readtype(io, Int16)
    if Int64(version) & Const.kByteCountMask > 0
        skip(io, 4)
    end
    fUniqueID = readtype(io, UInt32)
    fBits = readtype(io, UInt32)
    if fBits & Const.kIsReferenced > 0
        skip(io, 2)
    end
end

"""
    function endcheck(io, preamble::Preamble)

Checks if everything went well after parsing a TOBject. Used in conjuction
with `Preamble`.
"""
function endcheck(io, preamble::T) where {T<:Preamble}
    if !ismissing(preamble.cnt)
        observed = position(io) - preamble.start
        if observed != preamble.cnt
            error("Object '$(preamble.type)' has $(observed) bytes; expected $(preamble.cnt)")
        end
    end
    nothing
end

function nametitle(io)
    preamble = Preamble(io, Missing)
    skiptobj(io)
    name = readtype(io, String)
    title = readtype(io, String)
    endcheck(io, preamble)
    name, title
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :(gelu(0.1097699546178803 + 0.5720816873107069*gelu(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.11808627755560508*gelu(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.12730937872071246*gelu(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.6166642606025916*gelu(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.45996993866072433*gelu(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + gelu(-0.8288330741680832 + -0.5580861219922526*gelu(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.8786358511296881*gelu(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.6915444540383735*gelu(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.5554450106282389*gelu(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.20937059223327026*gelu(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + gelu(0.13615042538554878 + 0.08502572267265318*gelu(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.9980047609363689*gelu(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.5678523898140919*gelu(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.29618367003382495*gelu(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.7256587322261439*gelu(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + gelu(-0.16139551949919406 + -0.5886735107938796*gelu(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.02630736114123211*gelu(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.03441572289730743*gelu(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.3935924617025486*gelu(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.3425056157482871*gelu(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + gelu(-0.5329713534130147 + -0.3183971019485696*gelu(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.5847814320797045*gelu(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.7330877632484647*gelu(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + -0.37731199331975374*gelu(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.08817423127621415*gelu(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using SampleChains: SampleChains
using SampleChains.TupleVectors: TupleVectors
using SampleChainsDynamicHMC: SampleChainsDynamicHMC
using SampleChainsDynamicHMC.TransformVariables

# minimal AbstractChain implementation
struct TestChain{T} <: SampleChains.AbstractChain{T}
    samples::TupleVectors.TupleVector{T}
    info::AbstractVector
end
SampleChains.samples(chain::TestChain) = getfield(chain, :samples)
SampleChains.info(chain::TestChain) = getfield(chain, :info)

function samplechains_dynamichmc_sample(nchains, ndraws)
    # μ ~ Normal(0, 1) |> iid(2, 3), σ ~ HalfNormal(0, 1), y[i] ~ Normal(μ[i], σ)
    y = [0.74 0.15 -1.08; -0.42 1.08 -0.52]
    function ℓ(nt)
        μ = nt.μ
        σ = nt.σ
        return -(sum(abs2, μ) + σ^2 + sum((y .- μ) .^ 2) / σ^2) / 2 - length(y) * log(σ)
    end
    t = as((μ=as(Array, 2, 3), σ=asℝ₊))
    chain = SampleChains.newchain(nchains, SampleChainsDynamicHMC.dynamichmc(), ℓ, t)
    return SampleChains.sample!(chain, ndraws - 1)
end

@testset "SampleChains" begin
    @testset "TestChain with $nchains chains" for nchains in (1, 4)
        ndraws = 10
        dims = Dict(:y => [:a], :z => [:b, :c])
        coords = Dict(:a => ["a1", "a2"], :b => ["b1", "b2"], :c => ["c1", "c2", "c3"])
        tvs = map(1:nchains) do _
            init = (x=randn(), y=randn(2), z=randn(2, 3))
            tv = TupleVectors.TupleVector(undef, init, ndraws)
            copyto!(tv.x, randn(ndraws))
            copyto!(tv.y, [randn(2) for _ in 1:ndraws])
            copyto!(tv.z, [randn(2, 3) for _ in 1:ndraws])
            return tv
        end
        nts = map(collect, tvs)
        info = [(x=3, y=4) for _ in 1:ndraws]
        chains = [TestChain(tv, info) for tv in tvs]
        multichain = SampleChains.MultiChain(chains)

        kwargs = (dims=dims, coords=coords, library="MyLib")
        data = Dict(
            "Vector{AbstractChain}" => chains,
            "NTuple{N,AbstractChain}" => Tuple(chains),
            "MultiChain" => multichain,
        )
        if nchains === 1
            data["AbstractChain"] = only(chains)
        end

        @testset "$k" for (k, chaindata) in data
            @testset "$group" for group in (:posterior, :prior)
                if group === :posterior
                    idata = from_samplechains(chaindata; kwargs...)
                    idata_nt = from_namedtuple(nts; kwargs...)
                else
                    idata = from_samplechains(; group => chaindata, kwargs...)
                    idata_nt = from_namedtuple(; group => nts, kwargs...)
                end
                idata_conv = convert_to_inference_data(chaindata; group, kwargs...)
                map((idata, idata_conv)) do _idata
                    test_namedtuple_data(
                        _idata, group, propertynames(multichain), nchains, ndraws; kwargs...
                    )
                    @test only(ArviZ.groupnames(_idata)) === group
                end
                @testset for var_name in propertynames(multichain)
                    @test getproperty(getproperty(idata, group), var_name).values ==
                        getproperty(getproperty(idata_nt, group), var_name).values
                    @test getproperty(getproperty(idata, group), var_name).values ==
                        getproperty(getproperty(idata_conv, group), var_name).values
                end
            end
        end
    end
    @testset "SampleChainsDynamicHMC" begin
        expected_stats_vars = (
            :acceptance_rate, :n_steps, :diverging, :energy, :tree_depth, :turning
        )

        multichain = samplechains_dynamichmc_sample(4, 10)
        idata = convert_to_inference_data(multichain)
        @test sort([:posterior, :sample_stats]) == ArviZ.groupnames(idata)
        stats_vars = setdiff(map(Symbol, idata.sample_stats.variables), (:chain, :draw))
        missing_vars = setdiff(expected_stats_vars, stats_vars)
        @test isempty(missing_vars)
        post_vars = setdiff(map(Symbol, idata.posterior.variables), (:chain, :draw))
        @test issubset((:μ, :σ), post_vars)
        @test size(idata.posterior.μ.values) == (4, 10, 2, 3)
        @test size(idata.posterior.σ.values) == (4, 10)

        idata = convert_to_inference_data(multichain; group=:prior)
        @test sort([:prior, :sample_stats_prior]) == ArviZ.groupnames(idata)
        stats_vars = setdiff(
            map(Symbol, idata.sample_stats_prior.variables), (:chain, :draw)
        )
        missing_vars = setdiff(expected_stats_vars, stats_vars)
        @test isempty(missing_vars)
    end
end


text:
using MixedVolume
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "CayleyConfiguration" begin
    A1 = [0 1 0; 0 0 1]
    A2 = [1 0; 0 1]

    CC = CayleyConfiguration(A1, A2)
    @test CC.A == [0 1 0 1 0; 0 0 1 0 1; 1 1 1 0 0; 0 0 0 1 1]
    @test CC.offsets == [0, 3]
    @test size(CC) == (4, 5)
end

@testset "CircuitTable" begin
    A1 = [0 0 1 1; 0 2 0 1]
    A2 = [0 0 1 2; 0 1 1 0]
    A = CayleyConfiguration(A1, A2)

    M = [(2, 3), (1, 3)]

    T = CircuitTable(A, M)
    @test T isa CircuitTable{Int}

    w = [0, 0, 0, -1, 0, -3, -4, -8]

    @test dot(T, 1, w) == 8
    @test dot(T, 2, w) == 0
    @test dot(T, 3, w) == 0
    @test dot(T, 4, w) == -1
    @test dot(T, 4, [0, 0, 0, -2, 0, -3, -4, -8]) == 2
    @test dot(T, 5, w) == 0
    @test dot(T, 7, w) == 0
end

@testset "MixedCell" begin
    A1 = [0 0 1 1; 0 2 0 1]
    A2 = [0 0 1 2; 0 1 1 0]
    A = CayleyConfiguration(A1, A2)

    M = MixedCell(A, [(2, 3), (1, 3)])
    @test M isa MixedCell{Int}
end


text:
struct SmolyakDR{S<:Grid,T<:SmolyakGrid,nx} <: AbstractDecisionRule{S,T,nx}
    grid_exo::S
    grid_endo::T
    coefs::Vector{Matrix{Float64}}
end

#####
##### 1-argument decision rule
#####

function SmolyakDR(
        grid_exo::S, grid_endo::T, ::Union{Val{nx},Type{Val{nx}}}
    ) where S <: EmptyGrid where T <: SmolyakGrid where nx
    coefs = [Array{Float64}(undef,n_nodes(grid_endo), nx)]
    SmolyakDR{S,T,nx}(grid_exo, grid_endo, coefs)
end

function set_values!(
        dr::SmolyakDR{<:G}, values::Vector{Matrix{Float64}}
    ) where G <: Union{EmptyGrid,UnstructuredGrid}
    qnodes = qr(dr.grid_endo.B_nodes, Val(true))
    for i in 1:length(values)
        ldiv!(dr.coefs[i], qnodes , values[i])
    end
end

function set_values!(
        dr::SmolyakDR{<:G,<:SmolyakGrid,nx},
        values::Vector{<:Array{Value{nx}}}
    ) where G <: Union{EmptyGrid,UnstructuredGrid} where nx
    if length(values) != length(dr.coefs)
        msg = "The length of values ($(length(values))) is not the same "
        msg *= "as the length of the coefficient Vector ($(length(dr.coefs)))"
        error(msg)
    end
    qnodes = qr(dr.grid_endo.B_nodes,Val(true))
    for i in 1:length(values)
        N = length(values[i])
        data = reshape(reinterpret(Float64, vec(values[i])), (nx, N))'
        ldiv!(dr.coefs[i], qnodes, data)
    end
end

function evaluate(dr::SmolyakDR{<:EmptyGrid}, z::AbstractMatrix)
    B = BM.evalbase(dr.grid_endo.smol_params, z)
    B*dr.coefs[1]
end

function evaluate(dr::SmolyakDR{<:EmptyGrid,SmolyakGrid{d}}, points::AbstractVector{Point{d}}) where d
    N = length(points)
    mat = reshape(reinterpret(Float64, (points)), (d, N))'
    evaluate(dr, mat)
end

function evaluate(dr::SmolyakDR{<:EmptyGrid,SmolyakGrid{d}}, points::Point{d}) where d
    Point{d}(evaluate(dr, reshape([points...], (1,d)))[1,:]...)
end

####
#### UnstructuredGrid × CartesianGrid 2 continous arguments d.r.
####

function SmolyakDR(grid_exo::UnstructuredGrid, grid_endo::SmolyakGrid, ::Union{Val{nx},Type{Val{nx}}}) where nx
    coefs = [Array{Float64}(undef, n_nodes(grid_endo), nx) for i in 1:n_nodes(grid_exo)]
    SmolyakDR{typeof(grid_exo),typeof(grid_endo),nx}(grid_exo, grid_endo, coefs)
end

function evaluate(dr::SmolyakDR{<:UnstructuredGrid}, i::Int, z::AbstractMatrix)
    @boundscheck begin
        n_funcs = length(dr.coefs)
        if i > n_funcs
            msg = "Only $n_funcs are known, but function $i was requested"
            throw(BoundsError(msg))
        end
    end
    BM.evalbase(dr.grid_endo.smol_params, z) * dr.coefs[i]
end

function evaluate(dr::SmolyakDR{<:UnstructuredGrid}, i::Int, z::Vector{Point{d}}) where d
    N = length(z)
    mat = reshape(reinterpret(Float64, vec(z)), (d, N))'
    evaluate(dr, i, mat)
end


text:
export RBF2D, interp1, RBF3D

@doc raw"""
    function RBF2D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}};
        c::Union{PyObject, Array{Float64, 1}, Missing} = missing,
        eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
        d::Union{PyObject, Array{Float64, 1}} = zeros(0),
        kind::Int64 = 0)

Constructs a radial basis function representation on a 2D domain

$$f(x, y) = \sum_{i=1}^N c_i \phi(r; \epsilon_i) + d_0 + d_1 x + d_2 y$$

Here `d` can be either 0, 1 (only $d_0$ is present), or 3 ($d_0$, $d_1$, and $d_2$ are all present).

`kind` determines the type of radial basis functions

* 0:Gaussian

$$\phi(r; \epsilon) = e^{-(\epsilon r)^2}$$

* 1:Multiquadric

$$\phi(r; \epsilon) = \sqrt{1+(\epsilon r)^2}$$

* 2:Inverse quadratic

$$\phi(r; \epsilon) = \frac{1}{1+(\epsilon r)^2}$$

* 3:Inverse multiquadric

$$\phi(r; \epsilon) = \frac{1}{\sqrt{1+(\epsilon r)^2}}$$

Returns a callable struct, i.e. to evaluates the function at locations $(x, y)$ (`x` and `y` are both vectors), run
```julia
rbf(x, y)
```
"""
mutable struct RBF2D
    xc::PyObject
    yc::PyObject
    eps::PyObject
    c::PyObject
    d::PyObject
    kind::Int64

    function RBF2D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}};
            c::Union{PyObject, Array{Float64, 1}, Missing} = missing,
            eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
            d::Union{PyObject, Array{Float64, 1}} = zeros(0),
            kind::Int64 = 0)
        if isa(eps, Real)
            eps = eps * ones(length(xc))
        end
        c = coalesce(c, Variable(zeros(length(xc))))
        eps = coalesce(eps, ones(length(xc)))
        @assert length(xc)==length(yc)==length(c)==length(eps)
        @assert length(d) in [0,1,3]
        new(constant(xc), constant(yc), constant(eps), constant(c), constant(d), kind)
    end
end

function (o::RBF2D)(x,y)
    radial_basis_function_ = load_op_and_grad(libadcme,"radial_basis_function")
    x,y,xc,yc,eps,c,d,kind = convert_to_tensor(Any[x,y,o.xc,o.yc,o.eps,o.c,o.d,o.kind], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64])
    out = radial_basis_function_(x,y,xc,yc,eps,c,d,kind)
    set_shape(out, (length(x,)))
end

function Base.:show(io::IO, o::RBF2D)
    ct = hasproperty(o.c, :trainable) && o.c.trainable
    et = hasproperty(o.eps, :trainable) && o.eps.trainable
    dt = hasproperty(o.d, :trainable) && o.d.trainable
    print("RadialBasisFunction(NumOfCenters=$(length(o.xc)),NumOfAdditionalTerm=$(length(o.d)),CoeffIsVariable=$(ct),ShapeIsVariable=$(et),AdditionalTermIsVariable=$dt)")
end

"""
    interp1(x::Union{Array{Float64, 1}, PyObject},v::Union{Array{Float64, 1}, PyObject},xq::Union{Array{Float64, 1}, PyObject})

returns interpolated values of a 1-D function at specific query points using linear interpolation.
Vector x contains the sample points, and v contains the corresponding values, v(x).
Vector xq contains the coordinates of the query points.

!!! info
    `x` should be sorted in ascending order.

# Example
```julia
x = sort(rand(10))
y = constant(@. x^2 + 1.0)
z = [x[1]; x[2]; rand(5) * (x[end]-x[1]) .+ x[1]; x[end]]
u = interp1(x,y,z)
```
"""
function interp1(x::Union{Array{Float64, 1}, PyObject},v::Union{Array{Float64, 1}, PyObject},xq::Union{Array{Float64, 1}, PyObject})
    interp_dim_one_ = load_system_op("interp_dim_one")
    x,v,xq = convert_to_tensor(Any[x,v,xq], [Float64,Float64,Float64])
    out = interp_dim_one_(x,v,xq)
end

@doc raw"""
    RBF3D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}},
        zc::Union{PyObject, Array{Float64, 1}};
        c::Union{PyObject, Array{Float64, 1}, Missing} = missing,
        eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
        d::Union{PyObject, Array{Float64, 1}} = zeros(0),
        kind::Int64 = 0)

Constructs a radial basis function representation on a 3D domain

$$f(x, y, z) = \sum_{i=1}^N c_i \phi(r; \epsilon_i) + d_0 + d_1 x + d_2 y + d_3 z$$

Here `d` can be either 0, 1 (only $d_0$ is present), or 4 ($d_0$, $d_1$, $d_2$, and $d_3$ are all present).

`kind` determines the type of radial basis functions

* 0:Gaussian

$$\phi(r; \epsilon) = e^{-(\epsilon r)^2}$$

* 1:Multiquadric

$$\phi(r; \epsilon) = \sqrt{1+(\epsilon r)^2}$$

* 2:Inverse quadratic

$$\phi(r; \epsilon) = \frac{1}{1+(\epsilon r)^2}$$

* 3:Inverse multiquadric

$$\phi(r; \epsilon) = \frac{1}{\sqrt{1+(\epsilon r)^2}}$$

Returns a callable struct, i.e. to evaluates the function at locations $(x, y, z)$ (`x`, `y`, and `z` are all vectors), run
```julia
rbf(x, y, z)
```
"""
mutable struct RBF3D
    xc::PyObject
    yc::PyObject
    zc::PyObject
    eps::PyObject
    c::PyObject
    d::PyObject
    kind::Int64

    function RBF3D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}},
            zc::Union{PyObject, Array{Float64, 1}};
            c::Union{PyObject, Array{Float64, 1}, Missing} = missing,
            eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
            d::Union{PyObject, Array{Float64, 1}} = zeros(0),
            kind::Int64 = 0)
        if isa(eps, Real)
            eps = eps * ones(length(xc))
        end
        c = coalesce(c, Variable(zeros(length(xc))))
        eps = coalesce(eps, ones(length(xc)))
        @assert length(xc)==length(yc)==length(zc)==length(c)==length(eps)
        @assert length(d) in [0,1,4]
        new(constant(xc), constant(yc), constant(zc), constant(eps), constant(c), constant(d), kind)
    end
end

function (o::RBF3D)(x,y,z)
    radial_basis_function_ = load_op_and_grad(libadcme,"radial_basis_function_three_d")
    x,y,z,xc,yc,zc,eps,c,d,kind = convert_to_tensor(Any[x,y,z,o.xc,o.yc,o.zc, o.eps,o.c,o.d,o.kind], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64])
    out = radial_basis_function_(x,y,z,xc,yc,zc,eps,c,d,kind)
    set_shape(out, (length(x,)))
end

function Base.:show(io::IO, o::RBF3D)
    ct = hasproperty(o.c, :trainable) && o.c.trainable
    et = hasproperty(o.eps, :trainable) && o.eps.trainable
    dt = hasproperty(o.d, :trainable) && o.d.trainable
    print("RadialBasisFunction(NumOfCenters=$(length(o.xc)),NumOfAdditionalTerm=$(length(o.d)),CoeffIsVariable=$(ct),ShapeIsVariable=$(et),AdditionalTermIsVariable=$dt)")
end