text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test-4d.dat")
sname = joinpath(dirname(@__FILE__), "..", "data", "tmp", "same.dat")

x, y, z, s, t = read_dat(fname)

@test size(x) == (30,)
@test size(y) == (36,)
@test size(z) == (28,)
@test size(s) == (30,36,28,2)
@test size(t) == (2,)

@test maximum(x) == 72.5
@test maximum(y) == 71.220001
@test maximum(z) == 76.809998
@test maximum(s) == 0.067409396
@test maximum(t) == 0.24

@test minimum(x) == -72.5
@test minimum(y) == -103.779999
@test minimum(z) == -58.189999
@test minimum(s) == 0.0
@test minimum(t) == 0.12

write_dat(sname, x, y, z, s[:,:,:,:], t)

x2, y2, z2, s2, t2 = read_dat(sname)

@test x == x2
@test y == y2
@test z == z2
@test s == s2
@test t == t2

fname = joinpath(dirname(@__FILE__), "..", "data", "test-3d.dat")
x, y, z, s, t = read_dat(fname)

@test size(x) == (30,)
@test size(y) == (36,)
@test size(z) == (28,)
@test size(s) == (30,36,28,1)
@test size(t) == (1,)

@test maximum(x) == 72.5
@test maximum(y) == 71.220001
@test maximum(z) == 76.809998
@test maximum(s) == 33.2692985535
@test maximum(t) == 0

@test minimum(x) == -72.5
@test minimum(y) == -103.779999
@test minimum(z) == -58.189999
@test minimum(s) == -7.5189352036
@test minimum(t) == 0


text:
using DiffEqBayes
using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
using Stan, Mamba, Test

ProjDir = dirname(@__FILE__)

f1 = @ode_def_nohes LotkaVolterraTest4 begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=>1.0 c=>3.0 d=>1.0

cd(ProjDir) do

  isdir("tmp") && rm("tmp", recursive=true)
  isfile("LotkaVolterra2-1.pdf") &&
  rm("LotkaVolterra2-1.pdf")

  println("Four parameter case")
  u0 = [1.0,1.0]
  tspan = (0.0,10.0)
  prob1 = ODEProblem(f1,u0,tspan)
  sol = solve(prob1,Tsit5())
  t = collect(range(1, stop=10, length=10))
  randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
  data = convert(Array,randomized)
  priors = [Truncated(Normal(1.5,0.01),0,2),Truncated(Normal(1.0,0.01),0,1.5),Truncated(Normal(3.0,0.01),0,4),Truncated(Normal(1.0,0.01),0,2)]

  bayesian_result = stan_inference(prob1,t,data,priors;num_samples=100,num_warmup=500,vars =("StanODEData",InverseGamma(4,1)))
  theta1 = bayesian_result.chain_results[:,["theta.1"],:]
  theta2 = bayesian_result.chain_results[:,["theta.2"],:]
  theta3 = bayesian_result.chain_results[:,["theta.3"],:]
  theta4 = bayesian_result.chain_results[:,["theta.4"],:]

  global sim2 = bayesian_result.chain_results[1:end,
    ["theta.1", "theta.2", "theta.3", "theta.4"], :]
  ## Plotting
  p = plot(sim2, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, ncol=4, filename="LotkaVolterra2", fmt=:pdf)

  @test mean(theta1.value[:,:,1]) ≈ 1.5 atol=1e-1
  @test mean(theta2.value[:,:,1]) ≈ 1.0 atol=1e-1
  @test mean(theta3.value[:,:,1]) ≈ 3.0 atol=1e-1
  @test mean(theta4.value[:,:,1]) ≈ 1.0 atol=1e-1
end

ext: jl


text:
mutable struct Image <: GtkBox

    handle::Ptr{Gtk.GObject}
    data
    c::GtkCanvas

    function Image(img)

        data = array_to_rgb(img)
        c = GtkCanvas()
        set_gtk_property!(c, :expand, true)

        @guarded Gtk.ShortNames.draw(c) do widget
            xview, yview = guidata[widget, :xview], guidata[widget, :yview]
            set_coordinates(getgc(widget), xview, yview)

            roi = data[floor(Int, xview.min):ceil(Int, xview.max),
                       floor(Int, yview.min):ceil(Int, yview.max)]

            copy!(widget, roi)
        end

        b = GtkBox(:v)
        push!(b, c)
        # Initialize panning & zooming
        panzoom(c, (1, size(data, 1)), (1, size(data, 2)))
        panzoom_mouse(c, factor=1.0)
        panzoom_key(c)

        i = new(b.handle, data, c)
        Gtk.gobject_move_ref(i, b)

        signal_connect(image_key_press_cb, i, "key-press-event",
        Cint, (Ptr{Gtk.GdkEvent}, ), false, i)
        i
    end
end

array_to_rgb(img::Array{T, 2}) where {T<:Colors.Color} = img

function array_to_rgb(img::Array{T, 2}) where {T<:Number}
    data = Array{Colors.RGB24}(undef, size(img)...)
    img = img .- minimum(img)
    img = img ./ maximum(img)
    for i in eachindex(img)
        data[i] = Colors.RGB(img[i], img[i], img[i])
    end
    data
end
function array_to_rgb(img::Array{T, 3}) where {T<:Number}
    size(img, 3) != 3 && error("The size of the third dimension needs to be equal to 3 (RGB).")
    data = Array{Colors.RGB24}(undef, size(img)...)
    img = img .- minimum(img)
    img = img ./ maximum(img)
    for i = 1:size(img, 1)
        for j = 1:size(img, 2)
            data[i, j] = Colors.RGB(img[i, j, 1], img[i, j, 2], img[i, j, 3])
        end
    end
    data
end

function image(img)
    i = Image(reverse(img, dims=1))
    f = get_tab(fig_ntbook, index(fig_ntbook))
    if typeof(f) == Image
        idx = index(fig_ntbook)
        splice!(fig_ntbook, idx)
        insert!(fig_ntbook, idx, i, "Image")
    else
        idx = length(fig_ntbook)+1
        insert!(fig_ntbook, idx, i, "Image")
    end
    showall(fig_ntbook)
    index(fig_ntbook, idx)
    i
end

function save(f::Image, filename)
    w, h = width(f.c), height(f.c)
    surf = Cairo.CairoPDFSurface(filename, w, h)
    cr = Cairo.CairoContext(surf)
    Cairo.set_source_surface(cr, f.c.backcc.surface)
    Cairo.paint(cr)
    Cairo.finish(surf)
end

@guarded (PROPAGATE) function image_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    i = user_data
    event = convert(Gtk.GdkEvent, eventptr)
    if doing(Action(keyval("r"), ""), event)
        GtkUtilities.PanZoom.zoom_reset(i.c)
    end
    return PROPAGATE
end

@guarded (INTERRUPT) function fig_ntbook_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    ntbook = convert(GtkNotebook, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)

    if doing(Actions["newtab"], event)
        Immerse.figure()
    end
    if doing(Actions["closetab"], event)
        t = get_current_tab(ntbook)
        if typeof(t) == Figure
            Immerse.closefig(t.figno)
        elseif typeof(t) == Image
            close_tab(ntbook)
        end
    end
    return PROPAGATE
end

@guarded (nothing) function fig_ntbook_switch_page_cb(widgetptr::Ptr, pageptr::Ptr, pagenum::Int32, user_data)

    page = convert(Gtk.GtkWidget, pageptr)
    if typeof(page) == Figure
        Immerse.switchfig(_display, page.figno)
    end
    nothing
end

Base.show(io::IO, p::Gadfly.Plot) = write(io, "Gadfly.Plot(...)")

# # I'm not sure why but I need to call display here
function Base.show(io::IO, cc::Compose.Context)
    display(cc)
    write(io, "Compose.Context(...)")
end

function Immerse.figure(;name::AbstractString="Figure $(Immerse.nextfig(Immerse._display))",
                 width::Integer=400,    # TODO: make configurable
                 height::Integer=400)

    i = Immerse.nextfig(Immerse._display)
    f = Immerse.Figure()
    #Gtk.on_signal_destroy((x...)->Immerse.dropfig(Immerse._display, i), f)
    signal_connect(Immerse.on_figure_destroy, f, "destroy", Nothing, (), false, (i, Immerse._display))

    idx = length(fig_ntbook)+1
    insert!(fig_ntbook, idx, f, name)

    showall(fig_ntbook)
    Immerse.initialize_toolbar_callbacks(f)
    Immerse.addfig(Immerse._display, i, f)

    index(fig_ntbook, idx)
    i
end

function Immerse.figure(i::Integer; displayfig::Bool = true)

    Immerse.switchfig(_display, i)
    fig = Immerse.curfig(_display)
    displayfig && display(_display, fig)

    for idx = 1:length(fig_ntbook)
        f = fig_ntbook[idx]
        if typeof(f) == Figure && f.figno == i
            index(fig_ntbook, idx)
        end
    end
    fig
end

function Immerse.closefig(i::Integer)

    fig = Immerse.getfig(Immerse._display, i)
    for idx = 1:length(fig_ntbook)
        f = fig_ntbook[idx]
        if typeof(f) == Figure && f.figno == i

            Immerse.clear_hit(fig)
            close_tab(fig_ntbook, idx)

            destroy(fig)
            return
        end
    end
end


text:
function conditionally_add_symbols!(inputs_to_function::AbstractDict,
                                    lhs_assignments::OrderedCollections.OrderedDict, col)
    # if it's already been assigned at top-level,
    # don't add it to the inputs
    if haskey(lhs_assignments, col)
        return lhs_assignments[col]
    else
        return addkey!(inputs_to_function, col)
    end
end

replace_syms_astable!(inputs_to_function::AbstractDict,
                      lhs_assignments::OrderedCollections.OrderedDict, x) = x
replace_syms_astable!(inputs_to_function::AbstractDict,
                      lhs_assignments::OrderedCollections.OrderedDict, q::QuoteNode) =
    conditionally_add_symbols!(inputs_to_function, lhs_assignments, q)

function replace_syms_astable!(inputs_to_function::AbstractDict,
                               lhs_assignments::OrderedCollections.OrderedDict, e::Expr)
    if onearg(e, :^)
        return e.args[2]
    end

    col = get_column_expr(e)
    if col !== nothing
        return conditionally_add_symbols!(inputs_to_function, lhs_assignments, col)
    elseif e.head == :.
        return replace_dotted_astable!(inputs_to_function, lhs_assignments, e)
    else
        return mapexpr(x -> replace_syms_astable!(inputs_to_function, lhs_assignments, x), e)
    end
end

protect_replace_syms_astable!(inputs_to_function::AbstractDict,
                              lhs_assignments::OrderedCollections.OrderedDict, e) = e
protect_replace_syms_astable!(inputs_to_function::AbstractDict,
                              lhs_assignments::OrderedCollections.OrderedDict, e::Expr) =
    replace_syms!(inputs_to_function, lhs_assignments, e)

function replace_dotted_astable!(inputs_to_function::AbstractDict,
                                 lhs_assignments::OrderedCollections.OrderedDict, e)
    x_new = replace_syms_astable!(inputs_to_function, lhs_assignments, e.args[1])
    y_new = protect_replace_syms_astable!(inputs_to_function, lhs_assignments, e.args[2])
    Expr(:., x_new, y_new)
end

is_column_assigment(ex) = false
function is_column_assigment(ex::Expr)
    ex.head == :(=) && (get_column_expr(ex.args[1]) !== nothing)
end

# Taken from MacroTools.jl
# No docstring so assumed unstable
block(ex) = isexpr(ex, :block) ? ex : :($ex;)

sym_or_str_to_sym(x::Union{AbstractString, Symbol}) = Symbol(x)
sym_or_str_to_sym(x) =
    throw(ArgumentError("New columns created inside @astable must be Symbols or AbstractStrings"))

function get_source_fun_astable(ex; exprflags = deepcopy(DEFAULT_FLAGS))
    inputs_to_function = Dict{Any, Symbol}()
    lhs_assignments = OrderedCollections.OrderedDict{Any, Symbol}()

    # Make sure all top-level assignments are
    # in the args vector
    ex = block(MacroTools.flatten(ex))
    exprs = map(ex.args) do arg
        if is_column_assigment(arg)
            lhs = get_column_expr(arg.args[1])
            rhs = arg.args[2]
            new_ex = replace_syms_astable!(inputs_to_function, lhs_assignments, arg.args[2])
            if haskey(inputs_to_function, lhs)
                new_lhs = inputs_to_function[lhs]
                lhs_assignments[lhs] = new_lhs
            else
                new_lhs = addkey!(lhs_assignments, lhs)
            end

            Expr(:(=), new_lhs, new_ex)
        else
            replace_syms_astable!(inputs_to_function, lhs_assignments, arg)
        end
    end
    source = :(DataFramesMeta.make_source_concrete($(Expr(:vect, keys(inputs_to_function)...))))

    inputargs = Expr(:tuple, values(inputs_to_function)...)
    nt_iterator = (:(DataFramesMeta.sym_or_str_to_sym($k) => $v) for (k, v) in lhs_assignments)
    nt_expr = Expr(:tuple, Expr(:parameters, nt_iterator...))

    body = Expr(:block, Expr(:block, exprs...), nt_expr)

    fun = quote
        $inputargs -> begin
            $body
        end
    end

    # TODO: Add passmissing support by
    # checking if any input arguments missing,
    # and if-so, making a named tuple with
    # missing values
    if exprflags[BYROW_SYM][]
        fun = :(ByRow($fun))
    end

    return source, fun
end


text:
###################################################################################
#                             load dependencies
###################################################################################
cd(@__DIR__)
using Revise, DifferentialEvolutionMCMC, Random, Distributions
Random.seed!(50514)
###################################################################################
#                             define functions
###################################################################################
# returns prior log likelihood
function prior_loglike(μ, σ)
    LL = 0.0
    LL += logpdf(Normal(0, 1), μ)
    LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
    return LL
end

# function for initial values
function sample_prior()
    μ = rand(Normal(0, 1))
    σ = rand(truncated(Cauchy(0, 1), 0, Inf))
    return [μ,σ]
end

# likelihood function
function loglike(data, μ, σ)
    return sum(logpdf.(Normal(μ, σ), data))
end
###################################################################################
#                             generate data
###################################################################################
data = rand(Normal(0.0, 1.0), 50)
###################################################################################
#                             configure sampler
###################################################################################
# parameter names
names = (:μ,:σ)
# parameter bounds
bounds = ((-Inf,Inf),(0.0,Inf))

# model object
model = DEModel(;
    sample_prior,
    prior_loglike,
    loglike,
    data,
    names
)

# DEMCMC sampler object
de = DE(;sample_prior, bounds, burnin = 1000, Np = 6)
# number of interations per particle
n_iter = 2000
###################################################################################
#                             estimate parameters
###################################################################################
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)


text:
module TestDataFrameRow
    using Test, DataFrames
    using DataFrames: columns

    df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
                   b=[2.0, missing, 1.2, 2.0, missing, missing],
                   c=["A", "B", "C", "A", "B", missing],
                   d=CategoricalArray([:A, missing, :C, :A, missing, :C]))
    df2 = DataFrame(a = [1, 2, 3])

    @test names(DataFrameRow(df, 1)) == [:a, :b, :c, :d]

    #
    # Equality
    #
    @test_throws ArgumentError DataFrameRow(df, 1) == DataFrameRow(df2, 1)
    @test DataFrameRow(df, 1) != DataFrameRow(df, 2)
    @test DataFrameRow(df, 1) != DataFrameRow(df, 3)
    @test DataFrameRow(df, 1) == DataFrameRow(df, 4)
    @test DataFrameRow(df, 2) == DataFrameRow(df, 5)
    @test DataFrameRow(df, 2) != DataFrameRow(df, 6)

    # isless()
    df4 = DataFrame(a=[1, 1, 2, 2, 2, 2, missing, missing],
                    b=Union{Float64, Missing}[2.0, 3.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0],
                    c=[:B, missing, :A, :C, :D, :D, :A, :A])
    @test isless(DataFrameRow(df4, 1), DataFrameRow(df4, 2))
    @test !isless(DataFrameRow(df4, 2), DataFrameRow(df4, 1))
    @test !isless(DataFrameRow(df4, 1), DataFrameRow(df4, 1))
    @test isless(DataFrameRow(df4, 1), DataFrameRow(df4, 3))
    @test !isless(DataFrameRow(df4, 3), DataFrameRow(df4, 1))
    @test isless(DataFrameRow(df4, 3), DataFrameRow(df4, 4))
    @test !isless(DataFrameRow(df4, 4), DataFrameRow(df4, 3))
    @test isless(DataFrameRow(df4, 4), DataFrameRow(df4, 5))
    @test !isless(DataFrameRow(df4, 5), DataFrameRow(df4, 4))
    @test !isless(DataFrameRow(df4, 6), DataFrameRow(df4, 5))
    @test !isless(DataFrameRow(df4, 5), DataFrameRow(df4, 6))
    @test isless(DataFrameRow(df4, 7), DataFrameRow(df4, 8))
    @test !isless(DataFrameRow(df4, 8), DataFrameRow(df4, 7))

    # hashing
    @test hash(DataFrameRow(df, 1)) != hash(DataFrameRow(df, 2))
    @test hash(DataFrameRow(df, 1)) != hash(DataFrameRow(df, 3))
    @test hash(DataFrameRow(df, 1)) == hash(DataFrameRow(df, 4))
    @test hash(DataFrameRow(df, 2)) == hash(DataFrameRow(df, 5))
    @test hash(DataFrameRow(df, 2)) != hash(DataFrameRow(df, 6))

    # check that hashrows() function generates the same hashes as DataFrameRow
    df_rowhashes, _ = DataFrames.hashrows(Tuple(columns(df)), false)
    @test df_rowhashes == [hash(dr) for dr in eachrow(df)]

    # test incompatible frames
    @test_throws UndefVarError DataFrameRow(df, 1) == DataFrameRow(df3, 1)

    # test RowGroupDict
    N = 20
    d1 = rand(map(Int64, 1:2), N)
    df5 = DataFrame([d1], [:d1])
    df6 = DataFrame(d1 = [2,3])

    # test_group("group_rows")
    gd = DataFrames.group_rows(df5)
    @test length(unique(gd.groups)) == 2

    # getting groups for the rows of the other frames
    @test length(gd[DataFrameRow(df6, 1)]) > 0
    @test_throws KeyError gd[DataFrameRow(df6, 2)]
    @test isempty(DataFrames.findrows(gd, df6, (gd.df[1],), (df6[1],), 2))

    # grouping empty frame
    gd = DataFrames.group_rows(DataFrame(x=Int[]))
    @test length(unique(gd.groups)) == 0

    # grouping single row
    gd = DataFrames.group_rows(df5[1:1,:])
    @test length(unique(gd.groups)) == 1

    # getproperty, setproperty! and propertynames
    r = DataFrameRow(df, 1)
    @test Base.propertynames(r) == names(df)
    @test r.a === 1
    @test r.b === 2.0
    r.a = 2
    @test r.a === 2
    r.b = 1
    @test r.b === 1.0

    # getindex
    r = DataFrameRow(df, 1)
    @test r[:] == r

    # keys, values and iteration, size
    @test keys(r) == names(df)
    @test values(r) == (df[1, 1], df[1, 2], df[1, 3], df[1, 4])
    @test collect(pairs(r)) == [:a=>df[1, 1], :b=>df[1, 2], :c=>df[1, 3], :d=>df[1, 4]]

    @test haskey(r, :a)
    @test !haskey(r, :zzz)

    @test length(r) == 4
    @test ndims(r) == 1
    @test ndims(typeof(r)) == 1
    @test size(r) == (4,)
    @test size(r, 1) == 4
    @test_throws BoundsError size(r, 2)

    df = DataFrame(a=nothing, b=1)
    io = IOBuffer()
    show(io, DataFrameRow(df, 1))
    @test String(take!(io)) == "DataFrameRow (row 1)\na  \nb  1"

    # convert
    df = DataFrame(a=[1, missing], b=[2.0, 3.0])
    dfr = DataFrameRow(df, 1)
    @test convert(Vector, dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
    @test convert(Vector{Int}, dfr)::Vector{Int} == [1, 2]
    @test Vector(dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
    @test Vector{Int}(dfr)::Vector{Int} == [1, 2]

    # copy
    df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
                   b=[2.0, missing, 1.2, 2.0, missing, missing],
                   c=["A", "B", "C", "A", "B", missing])
    @test copy(DataFrameRow(df, 1)) == (a = 1, b = 2.0, c = "A")
    @test isequal(copy(DataFrameRow(df, 2)), (a = 2, b = missing, c = "B"))

    # parent and parentindices
    @test parent(df[1, :]) === df
    @test parentindices(df[1, :]) == (1, Base.OneTo(3))
    @test parent(df[1, 1:3]) !== df
    @test parentindices(df[1, 1:3]) == (1, Base.OneTo(3))

    # iteration and collect
    ref = ["a", "b", "c"]
    df = DataFrame(permutedims(ref))
    dfr = df[1, :]
    @test collect(dfr) == ref
    @test eltype(collect(dfr)) === String
    for (v1, v2) in zip(ref, dfr)
        @test v1 == v2
    end
    for (i, v) in enumerate(dfr)
        @test v == ref[i]
    end
    dfr = DataFrame(a=1, b=true, c=1.0)[1,:]
    @test eltype(collect(dfr)) === Real
end


text:
@inbounds for i = ix:ix-1+n
    x[i] = $op( y[i], z[i] )
end


text:
# Driver functions for the integrators.

struct MolecularDynamicsParameters
    "Centroid friction (1 / ps)."
    gamma0::Float64
    "Time step (ps)."
    dt::Float64
    "Total duration (ps)."
    duration::Float64
end

"""
    num_steps(mdp::MolecularDynamicsParameters)

Number of molecular dynamics steps for `mdp`.
"""
num_steps(mdp::MolecularDynamicsParameters) = ceil(Int, mdp.duration / mdp.dt)

function run_constrained_equilibration(
        model_T::Type{<:WaterModel}, int_T::Type{<:Integrator}, N::Int,
        T::Float64, R::Float64, P::Int, mdp_equil::MolecularDynamicsParameters;
        qs_init::Maybe{Vector{Float64}}=nothing, output_skip::Int=100)
    is_constrained(int_T) || error("Integrator is not constrained.")

    # Number of atoms in each water molecule.
    A = 3
    # Number of spatial dimensions.
    D = 3

    if !isnothing(qs_init)
        s = State(P, D, A, N; qs_init)
    else
        s = State(P, D, A, N; R)
    end

    model = model_T(N)
    fc = ForceContainer(P, D, A, N)
    int_equil = int_T(model, R, mdp_equil.dt, beta(T), mdp_equil.gamma0, fc, s)

    for i in 1:num_steps(mdp_equil)
        step!(s, int_equil)

        i % output_skip == 0 || continue

        time = i * mdp_equil.dt
        E = energy(int_equil, s)
        E_K = pop!(E, "kinetic")
        E_Vspring = pop!(E, "spring")
        E_Vmodel = pop!(E, "model")
        isempty(E) || error("Unexpected energy contribution.")
        E_total = E_K + E_Vspring + E_Vmodel
        T = get_temperature(int_equil)

        println("$(time) $(E_K) $(E_Vspring) $(E_Vmodel) $(E_total) $(T)")
    end

    s
end

function run_constrained_production(
        model_T::Type{<:WaterModel}, int_T::Type{<:Integrator}, N::Int,
        T::Float64, R::Float64, P::Int, mdp_equil::MolecularDynamicsParameters,
        mdp_prod::MolecularDynamicsParameters;
        qs_init::Maybe{Vector{Float64}}=nothing, aux_output::IO=stderr)
    is_constrained(int_T) || error("Integrator is not constrained.")

    # Number of atoms in each water molecule.
    A = 3
    # Number of spatial dimensions.
    D = 3

    # Number of molecular dynamics steps for equilibration.
    N_equil = num_steps(mdp_equil)
    # Number of molecular dynamics steps for production.
    N_prod = num_steps(mdp_prod)

    if !isnothing(qs_init)
        s = State(P, D, A, N; qs_init)
    else
        s = State(P, D, A, N; R)
    end

    model = model_T(N)
    fc = ForceContainer(P, D, A, N)
    int_equil = int_T(model, R, mdp_equil.dt, beta(T), mdp_equil.gamma0, fc, s)
    int_prod = int_T(model, R, mdp_prod.dt, beta(T), mdp_prod.gamma0, fc, s)

    meter_equil = Progress(N_equil, output=aux_output)

    for _ in ProgressWrapper(1:N_equil, meter_equil)
        step!(s, int_equil)
    end

    constraints = Array{Float64}(undef, N_prod)
    E_Ks = Array{Float64}(undef, N_prod)
    E_Vsprings = Array{Float64}(undef, N_prod)
    E_Vmodels = Array{Float64}(undef, N_prod)
    E_totals = Array{Float64}(undef, N_prod)
    temperatures = Array{Float64}(undef, N_prod)

    forces1 = Array{Float64}(undef, N_prod)
    forces2 = Array{Float64}(undef, N_prod)

    estim1 = estimator1(R, beta(T), int_prod, s)
    estim2 = estimator2(R, beta(T), int_prod, s)

    meter = Progress(N_prod, output=aux_output)

    for i in ProgressWrapper(1:N_prod, meter)
        step!(s, int_prod)

        constraints[i] = com_r(CN1, CN2, CJ, s)
        E = energy(int_prod, s)
        E_Ks[i] = pop!(E, "kinetic")
        E_Vsprings[i] = pop!(E, "spring")
        E_Vmodels[i] = pop!(E, "model")
        isempty(E) || error("Unexpected energy contribution.")
        E_totals[i] = E_Ks[i] + E_Vsprings[i] + E_Vmodels[i]
        temperatures[i] = get_temperature(int_prod)

        forces1[i] = estim1()
        forces2[i] = estim2()
    end

    data_table = []

    for (label, data) in [("xi (nm)", constraints),
                          ("K (kJ / mol)", E_Ks),
                          ("Vspring (kJ / mol)", E_Vsprings),
                          ("Vmodel (kJ / mol)", E_Vmodels),
                          ("E (kJ / mol)", E_totals),
                          ("T (K)", temperatures),
                          ("est1 (1 / nm)", forces1),
                          ("est2 (1 / nm)", forces2)]
        data_mean, data_stderr, data_actime = aggregate(mdp_prod.dt, data)
        push!(data_table, [label, data_mean, data_stderr, data_actime])
    end

    pretty_table(aux_output, permutedims(hcat(data_table...)),
                 ["", "Value", "Error", "Autocorrelation time (ps)"])

    s, aggregate(mdp_prod.dt, forces1), aggregate(mdp_prod.dt, forces2)
end

function run_restrained_equilibration(
        model_T::Type{<:WaterModel}, int_T::Type{<:Integrator}, N::Int,
        T::Float64, R::Float64, P::Int, mdp_equil::MolecularDynamicsParameters,
        bias_force_constant::Float64; qs_init::Maybe{Vector{Float64}}=nothing,
        output_skip::Int=100)
    is_constrained(int_T) && error("Integrator is constrained.")

    # Number of atoms in each water molecule.
    A = 3
    # Number of spatial dimensions.
    D = 3

    if !isnothing(qs_init)
        s = State(P, D, A, N; qs_init)
    else
        s = State(P, D, A, N; R)
    end

    bias = HarmonicRestraint(bias_force_constant, R, CN1, CN2, CJ)
    model = model_T(N)
    fc = ForceContainer(P, D, A, N)
    int_equil = int_T(model, mdp_equil.dt, beta(T), mdp_equil.gamma0, fc, s;
                      extra_forces=Dict("bias" => [bias]))

    for i in 1:num_steps(mdp_equil)
        step!(s, int_equil)

        i % output_skip == 0 || continue

        time = i * mdp_equil.dt
        xi = com_r(CN1, CN2, CJ, s)
        E = energy(int_equil, s)
        E_K = pop!(E, "kinetic")
        E_Vspring = pop!(E, "spring")
        E_Vmodel = pop!(E, "model")
        E_Vbias = pop!(E, "bias")
        isempty(E) || error("Unexpected energy contribution.")
        E_total = E_K + E_Vspring + E_Vmodel + E_Vbias
        T = get_temperature(int_equil)

        print("$(time) $(xi) $(E_K) $(E_Vspring) $(E_Vmodel) $(E_Vbias) ")
        println("$(E_total) $(T)")
    end

    s
end

function run_restrained_production(
        model_T::Type{<:WaterModel}, int_T::Type{<:Integrator}, N::Int,
        T::Float64, R::Float64, P::Int, mdp_equil::MolecularDynamicsParameters,
        mdp_prod::MolecularDynamicsParameters, bias_force_constant::Float64;
        qs_init::Maybe{Vector{Float64}}=nothing, output_skip::Int=100)
    is_constrained(int_T) && error("Integrator is constrained.")

    # Number of atoms in each water molecule.
    A = 3
    # Number of spatial dimensions.
    D = 3

    # Number of molecular dynamics steps for equilibration.
    N_equil = num_steps(mdp_equil)
    # Number of molecular dynamics steps for production.
    N_prod = num_steps(mdp_prod)

    if !isnothing(qs_init)
        s = State(P, D, A, N; qs_init)
    else
        s = State(P, D, A, N; R)
    end

    bias = HarmonicRestraint(bias_force_constant, R, CN1, CN2, CJ)
    model = model_T(N)
    fc = ForceContainer(P, D, A, N)
    int_equil = int_T(model, mdp_equil.dt, beta(T), mdp_equil.gamma0, fc, s;
                      extra_forces=Dict("bias" => [bias]))
    int_prod = int_T(model, mdp_prod.dt, beta(T), mdp_prod.gamma0, fc, s;
                     extra_forces=Dict("bias" => [bias]))

    for _ in 1:N_equil
        step!(s, int_equil)
    end

    for i in 1:N_prod
        step!(s, int_prod)

        i % output_skip == 0 || continue

        println(com_r(CN1, CN2, CJ, s))
    end

    s
end


text:
# Documentation
@testset "Documentation" begin
    using Documenter, Econometrics, Weave, StatsBase

    prefix = ".."

    for file in filter!(file -> endswith(file, ".jmd"), readdir(joinpath(prefix, "docs", "jmd"), join = true))
        weave(file,
              out_path = joinpath(prefix, "docs", "src"),
              doctype = "github",
              )
    end

    DocMeta.setdocmeta!(Econometrics,
                       :DocTestSetup,
                       :(using Econometrics, Documenter, CSV, RDatasets, StatsBase;
                         ENV["COLUMNS"] = 120;
                         ENV["LINES"] = 30;),
                       recursive = true)
    # doctest(Econometrics, fix = true)
    makedocs(sitename = "Econometrics",
             format = Documenter.HTML(assets = [joinpath("assets", "custom.css")]),
             modules = [Econometrics, StatsBase],
             pages = [
                 "Introduction" => "index.md",
                 "Getting Started" => "getting_started.md",
                 "Estimators" => "estimators.md",
                 "API" => "api.md",
                 ],
             source = joinpath(prefix, "docs", "src"),
             build = joinpath(prefix, "docs", "build"),
             )
    @test true
end


text:
""" Schedules

This module will contain the tooling which converts diagrams to DWDs. This will
also include optimization tooling relating to DWD schedules (compressing matrix
operations, pre-computing constant values, moving addition to reduce linear).
"""
module Schedules
using Catlab.Syntax
using Catlab.WiringDiagrams.DirectedWiringDiagrams
using Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.Programs.DiagrammaticPrograms: NamedGraph

using ..Diagrams

export diag2dwd

sp_otimes(expr) = expr isa ObExpr{:otimes} ? expr.args : [expr]
n_args(expr) = length(sp_otimes(expr))

isa_proj(hom) = startswith(string(hom), "proj")
proj_ind(hom) = unsub(split(string(hom), "_")[1][5:end])

super2int = Dict(
  '⁰'=>'0', '¹'=>'1', '²'=>'2', '³'=>'3', '⁴'=>'4', '⁵'=>'5',
  '⁶'=>'6', '⁷'=>'7', '⁸'=>'8', '⁹'=>'9'
)
sub2int = Dict(
  '₀'=>'0', '₁'=>'1', '₂'=>'2', '₃'=>'3', '₄'=>'4', '₅'=>'5',
  '₆'=>'6', '₇'=>'7', '₈'=>'8', '₉'=>'9'
)

unsub(str) = begin
  parse(Int64, join([sub2int[a] for a in str]))
end

function eval_deps!(dwd, graph, el, w2b, el2p, obs; in_els = Dict{Int, Int}(), boundaries = [])
  if el in keys(el2p)
    return el2p[el]
  end

  inputs = incident(graph, el, :tgt)

  labels = graph[:ename]
  el_types = sp_otimes(obs[el])
  outp = if length(inputs) > 1
    sort!(inputs, by = i->proj_ind(labels[i]))
    vcat(map(w -> eval_deps!(dwd, graph, graph[w, :src], w2b, el2p, obs; in_els = in_els), inputs)...)
  elseif length(inputs) == 0
    el ∈ keys(in_els) || error("Element $el has no dependencies, but is not defined in `in_els`")
    [(input_id(dwd), in_els[el])]
  else
    in_ps = eval_deps!(dwd, graph, graph[inputs[1], :src], w2b, el2p, obs; in_els = in_els)
    if isa_proj(graph[inputs[1], :ename])
      [in_ps[proj_ind(graph[inputs[1], :ename])]]
    else
      wires = map(enumerate(in_ps)) do (ip, op)
        Wire(port_value(dwd, Port(op[1], OutputPort, op[2])), op, (w2b[inputs[1]], ip))
      end
      add_wires!(dwd, wires)
      map(1:length(output_ports(dwd, w2b[inputs[1]]))) do p
        (w2b[inputs[1]], p)
      end
    end
  end

  if el in keys(el2p)
    return el2p[el]
  end

  # Evaluate Boundary Conditions
  bcs = incident(graph, el, :src)
  filter!(a->"$(Symbol(graph[a, :ename]))"[1] == '∂', bcs)
  for bc in bcs
    wires = map(enumerate(outp)) do (ip, op)
      Wire(port_value(dwd, Port(op[1], OutputPort, op[2])), op, (w2b[bc], ip))
    end
    add_wires!(dwd, wires)
    outp .= map(1:length(output_ports(dwd, w2b[bc]))) do p
      (w2b[bc], p)
    end
  end
  el2p[el] = outp
  outp
end

name(a::HomExpr) = head(a) == :generator ? args(a)[1] : head(a)

function diag2dwd(diagram; clean = false, calc_states = [], out_vars=[], in_vars=[])
  homs = Vector{Any}(copy(diagram.hom_map))
  obs = Vector{Any}(copy(diagram.ob_map))
  graph = NamedGraph{Any, Any}()
  copy_parts!(graph, dom(diagram).graph)
  graph[:ename] .= homs

  # Expand homs which are composition
  for h in parts(graph, :E)
    h_name = graph[h, :ename]
    elsrc, eltgt = graph[h, :src], graph[h, :tgt]
    if h_name isa HomExpr{:compose}
      args = h_name.args
      rem_part!(graph, :E, h)
      verts = add_parts!(graph, :V, length(args) - 1, vname = [Symbol("anon_", gensym()) for i in 1:(length(args)-1)])
      append!(obs, codom.(args)[1:(end-1)])
      add_parts!(graph, :E, length(args), ename = args,
                          src = vcat([elsrc], verts),
                          tgt = vcat(verts, [eltgt]))
    end
  end

  flipped = fill(false, ne(graph))
  # Flip projections for graph eval
  for (i, h) in enumerate(homs)
    if isa_proj(h) && !flipped[i]
      iw = copy(incident(graph, graph[i, :src], :tgt))
      ow = copy(incident(graph, graph[i, :src], :src))
      if isempty(iw)
        for w in ow
          if isa_proj(homs[w]) && !flipped[w]
            wsrc = graph[w, :src]
            graph[w, :src] = graph[w, :tgt]
            graph[w, :tgt] = wsrc
            flipped[w] = true
          end
        end
      end
    end
  end

  pres = presentation(codom(diagram))

  # Decide on State and Parameter variables
  time_arrs = findall(h-> h isa HomExpr{:∂ₜ}, graph[:ename])
  params = findall(e -> isempty(incident(graph, e, :tgt)), parts(graph, :V))
  state_vals = graph[time_arrs, :src]

  # Mix State and Parameters into single in_els (need to ensure alignment between input and outputs
  in_els = unique(vcat(state_vals, params))
  out_els = copy(graph[time_arrs, :tgt])
  out_names = [Symbol(:∂ₜ, graph[v, :vname]) for v in state_vals]

  # TODO: Establish some kind of consistency expectations here
  if !isempty(in_vars) && !isempty(out_vars)
    in_els = [incident(graph, a, :vname) for a in in_vars]
    out_els = [incident(graph, a, :vname) for a in out_vars]
    out_names = graph[out_els, :vname]
  elseif !isempty(in_vars)
    # This ensures that the output values are in-order with input values
    new_inds = [findfirst(v -> graph[v, :vname] == a, in_els) for a in in_vars]
    non_param_inds = filter(i -> in_els[i] in state_vals , new_inds)
    in_els .= in_els[new_inds]
    out_els .= out_els[non_param_inds]
    out_names .= out_names[non_param_inds]
  elseif !isempty(out_vars)
    out_els = [incident(graph, a, :vname) for a in out_vars]
    out_names = graph[out_els, out_els]
  end
  if !isempty(calc_states)
    calc_inds = findall(v->graph[v, :vname] ∈ calc_states, state_vals)
    out_els = out_els[calc_inds]
    out_names = out_names[calc_inds]
  end

  in_types = [Dict(:name => graph[v, :vname],
                   :type => obs[v])
              for v in in_els]
  out_types = [Dict(:name => out_names[v],
                    :type => obs[out_els[v]])
               for v in 1:length(out_els)]

  # FIXME: Hacky solution to ensuring time_arrs aren't included later
  # need more thoughtful approach to sorting through edges here
  rem_parts!(graph, :E, time_arrs)

  dwd = WiringDiagram(collect(in_types), collect(out_types))

  # Add all necessary boxes to the DWD based on arrows in the diagram
  w2b = map(parts(graph, :E)) do a
    w_type = graph[a, :ename]
    # Fix this if-case. This was meant for times when multiple arguments would
    # go into a single box
    if isa_proj(w_type)
      nothing
    elseif w_type isa HomExpr
      src_el = graph[a, :src]
      tgt_el = graph[a, :tgt]
      in_ports = sp_otimes(obs[src_el])
      out_ports = sp_otimes(obs[tgt_el])
      add_box!(dwd, Box(name(w_type), [Dict(:type=>ip) for ip in in_ports],
                        [Dict(:type=>op) for op in out_ports]))
    end
  end

  # Mapping from elements (vertices in `graph`) to ports in the DWD
  # This ensures that elements which have already been computed are
  # not recomputed
  el2p = Dict{Int, Vector}()
  in_ps = vcat(map(out_els) do el
    eval_deps!(dwd, graph, el, w2b, el2p, obs; in_els=Dict(in_els[i] => i for i in 1:length(in_els)))
  end...)

  wires = map(enumerate(in_ps)) do (ip, op)
    Wire(out_types[ip], op, (output_id(dwd), ip))
  end
  add_wires!(dwd, wires)

  # Remove any boxes without any connecting wires
  if clean
    to_rem = Vector{Int64}()
    for b in 1:nparts(dwd.diagram, :Box)
      if isempty(vcat(incident(dwd.diagram, b, [:tgt, :in_port_box]),
                      incident(dwd.diagram, b, [:in_tgt, :in_port_box]))) &&
         isempty(vcat(incident(dwd.diagram, b, [:src, :out_port_box]),
                      incident(dwd.diagram, b, [:out_src, :out_port_box])))
         push!(to_rem, b)
      end
    end
    rem_boxes!(dwd, to_rem)
  end

  dwd
end

end


text:
using MicroLogging
using Test

import MicroLogging: BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
    shouldlog, handle_message, min_enabled_level, catch_exceptions,
    configure_logging

import Test: TestLogger

#-------------------------------------------------------------------------------
@testset "Logging" begin

include("config.jl")
include("ConsoleLogger.jl")
include("StickyMessages.jl")

end


text:
fm = Z ~ X + Y
df = DataFrame(X = randn(10), Y = randn(10), Z = randn(10));
mf = ModelFrame(Z ~ X + Y, df);


text:
using LowRankIntegrators, SparseArrays

@testset "Burgers equation" begin
    n = 1000 # spatial discretization

    l = π # length of spatial domain
    Δx = l/n # step size
    x_range = Δx/2:Δx:l-Δx/2 # uniform grid

    # boundary conditions
    left(i) = i > 1 ? i - 1 : n
    right(i) = i < n ? i + 1 : 1

    # discretized diff operators
    # laplacian (times viscosity)
    Δ = spzeros(n, n)
    ν = 0.005
    for i in 1:n
        Δ[i,left(i)] = ν/Δx^2
        Δ[i,i] = -2ν/Δx^2
        Δ[i,right(i)] = ν/Δx^2
    end

    # gradient
    ∇ = spzeros(n, n)
    for i in 1:n
        ∇[i,left(i)] = -1/2/Δx
        ∇[i,right(i)] = 1/2/Δx
    end

    function burgers(ρ, (∇,Δ), t)
        return Δ*ρ - (∇*ρ) .* ρ
    end

    # uncertainty range
    m = 20 # parameter realizations scale as m^2
    σ = [0.5,0.5]
    ξ_range = [(ξ_1,ξ_2) for ξ_1 in range(-1,1,length=m), ξ_2 in range(-1,1,length=m)];

    #initial condition
    ub(x) = 0.5*(exp.(cos.(x)) .- 1.5).*sin.(x .+ 2π*0.37) # deterministic initial condition
    uprime(x,ξ,σ) = σ[1]*ξ[1]*sin.(2π*x) .+ σ[2]*ξ[2]*sin.(3π*x) # stochastic fluctuation
    ρ0_mat = hcat([ub(x_range) + uprime(x_range, ξ, σ) for ξ in ξ_range]...) # full rank initial condition
    r = 5 # approximation rank
    lr_ρ0 = truncated_svd(ρ0_mat, r); # intial condition

    dt = 1e-2 # time step
    lr_prob = MatrixDEProblem((ρ,t) -> burgers(ρ, (∇,Δ), t),  lr_ρ0, (0.0, 1.0)) # defines the matrix differential equation problem
    solvers = [UnconventionalAlgorithm(),
               ProjectorSplitting(PrimalLieTrotter()),
               ProjectorSplitting(DualLieTrotter()),
               ProjectorSplitting(Strang()),
               RankAdaptiveUnconventionalAlgorithm(1e-4, rmax=10)]
    for solver in solvers
        lr_sol = LowRankIntegrators.solve(lr_prob, solver, dt) # solves the low rank approximation
    end
end


text:
##### Multi-dimensional interpolations

# TODO: Add implementations that perform interps with a single pass and no alloc

#####
##### CellCenter <-> CellEdge
#####

function interp!(cent::CellCenter, edge::CellEdge{dim}, grid::Grid) where {dim}
    face_orth = CellFace(grid, first(orthogs(dim)))
    interp!(face_orth, edge, grid)
    interp!(cent, face_orth, grid)
end

function interp!(edge::CellEdge{dim}, cent::CellCenter, grid::Grid) where {dim}
    face_orth = CellFace(grid, first(orthogs(dim)))
    interp!(face_orth, cent, grid)
    interp!(edge, face_orth, grid)
end

#####
##### CellCorner <-> CellFace
#####

function interp!(corn::CellCorner, face::CellFace{dim}, grid::Grid) where {dim}
    edge_orth = CellEdge(grid, first(orthogs(dim)))
    interp!(edge_orth, face, grid)
    interp!(corn, edge_orth, grid)
end

function interp!(face::CellFace{dim}, corn::CellCorner, grid::Grid) where {dim}
    edge_orth = CellEdge(grid, first(orthogs(dim)))
    interp!(edge_orth, corn, grid)
    interp!(face, edge_orth, grid)
end

#####
##### CellFace <-> CellEdge
#####

function interp!(face::CellFace{dim}, edge::CellEdge{dim}, grid::Grid) where {dim}
    face_orth = CellFace(grid, first(orthogs(dim)))
    cent = CellCenter(grid)
    interp!(face_orth, edge, grid)
    interp!(cent, face_orth, grid)
    interp!(face, cent, grid)
end

function interp!(edge::CellEdge{dim}, face::CellFace{dim}, grid::Grid) where {dim}
    cent = CellCenter(grid)
    face_orth = CellFace(grid, first(orthogs(dim)))
    interp!(cent, face, grid)
    interp!(face_orth, cent, grid)
    interp!(edge, face_orth, grid)
end

#####
##### CellCenter <-> CellCorner
#####

function interp!(cent::CellCenter, corn::CellCorner, grid::Grid)
    edge = CellEdge(grid, 1)
    interp!(edge, corn, grid)
    interp!(cent, edge, grid)
end

function interp!(corn::CellCorner, cent::CellCenter, grid::Grid)
    edge = CellEdge(grid, 1)
    interp!(edge, cent, grid)
    interp!(corn, edge, grid)
end


text:
include("robot_operations.jl")
pygui(true)

function mark_kross_x(r::Robot)
    for side in ((Nord,Ost),(Sud,Ost),(Sud,West),(Nord,West))
        putmarkers!(r,side)
        move_by_markers!(r,inverse(side))
    end
    putmarker!(r)
end

function putmarkers!(r::Robot,side::NTuple{2,HorizonSide})
    while isborderTulpe(r,side)==false
        moveTulpe(r,side)
        putmarker!(r)
    end
end

function isborderTulpe(r::Robot,side::NTuple{2,HorizonSide})
    return (isborder(r,side[1]) || isborder(r,side[2]))
end

function moveTulpe(r::Robot, side::NTuple{2,HorizonSide})
    for s in side
        move!(r,s)
    end
end

function move_by_markers!(r::Robot,side::NTuple{2,HorizonSide})
    while ismarker(r)
        moveTulpe(r,side)
    end
end

function inverse(side::NTuple{2,HorizonSide})
    return (inverseSide(side[1]), inverseSide(side[2]))
end

r=Robot(animate=true)
mark_kross_x(r)

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
Base.zero(x::arb_series) = arb_series(x.poly.parent(0), x.length)

"""
    bessel_j(ν::arb, z::arb_series[, n = length(z)])

> Compute the Taylor series of the Bessel function.

It's computed using a recursive formula for the Taylor coefficients
of the Bessel functions and then composing it with `z`.
"""
function bessel_j(ν::arb, z::arb_series, n = length(z))
    res = arb_series(parent(z.poly)(), n)

    if n > 0
        x = z[0]
        a0 = bessel_j(ν, x)
        res[0] = a0

        if !isfinite(a0)
            return arb_series(parent(z.poly)(fill(NaN, n)), n)
        end

        if n > 1
            a1 = 1 // 2 * (bessel_j(ν - 1, x) - bessel_j(ν + 1, x))
            res[1] = a1

            if !isfinite(a0)
                return arb_series(parent(z.poly)(fill(NaN, n)), n)
            end

            if n > 2
                a2 = 1 // 8 * (bessel_j(ν - 2, x) + bessel_j(ν + 2, x) - 2a0)
                res[2] = a2

                if !isfinite(a0)
                    return arb_series(parent(z.poly)(fill(NaN, n)), n)
                end

                if n > 3
                    a3 = 1 // 48 * (bessel_j(ν - 3, x) - bessel_j(ν + 3, x) - 6a1)
                    res[3] = a3

                    if !isfinite(a0)
                        return arb_series(parent(z.poly)(fill(NaN, n)), n)
                    end

                    if n > 4
                        x2 = x^2
                        ν2 = ν^2
                        for i = 4:n-1
                            k = i - 4
                            ai =
                                -(
                                    res[k] +
                                    2x * res[k+1] +
                                    (k^2 + 4k + 4 - ν2 + x2) * res[k+2] +
                                    (2k^2 + 11k + 15) * x * res[k+3]
                                ) / (x2 * (k^2 + 7k + 12))
                            res[i] = ai
                        end
                    end
                end
            end
        end
    end
    # Compose the Taylor series for the Bessel function with that of z
    z_tmp = arb_series(deepcopy(z.poly))
    z_tmp[0] = base_ring(parent(z.poly))(0)

    return Nemo.compose(res, z_tmp, n)
end

"""
    bessel_y(ν::arb, z::arb_series[, n = length(z)])

Compute the Taylor series of the Bessel function.

It's computed using a recursive formula for the Taylor coefficients of
the Bessel functions and then composing it with `z`. Note that since
both `bessel_j` and `bessel_y` satisfy the same differential equation
the recursive formula is the same.
"""
function bessel_y(ν::arb, z::arb_series, n = length(z))
    res = arb_series(parent(z.poly)(), n)

    if n > 0
        x = z[0]
        a0 = bessel_y(ν, x)
        res[0] = a0

        if !isfinite(a0)
            return arb_series(parent(z.poly)(fill(NaN, n)), n)
        end

        if n > 1
            a1 = 1 // 2 * (bessel_y(ν - 1, x) - bessel_y(ν + 1, x))
            res[1] = a1

            if !isfinite(a0)
                return arb_series(parent(z.poly)(fill(NaN, n)), n)
            end

            if n > 2
                a2 = 1 // 8 * (bessel_y(ν - 2, x) + bessel_y(ν + 2, x) - 2a0)
                res[2] = a2

                if !isfinite(a0)
                    return arb_series(parent(z.poly)(fill(NaN, n)), n)
                end

                if n > 3
                    a3 = 1 // 48 * (bessel_y(ν - 3, x) - bessel_y(ν + 3, x) - 6a1)
                    res[3] = a3

                    if !isfinite(a0)
                        return arb_series(parent(z.poly)(fill(NaN, n)), n)
                    end

                    if n > 4
                        x2 = x^2
                        ν2 = ν^2
                        for i = 4:n-1
                            k = i - 4
                            ai =
                                -(
                                    res[k] +
                                    2x * res[k+1] +
                                    (k^2 + 4k + 4 - ν2 + x2) * res[k+2] +
                                    (2k^2 + 11k + 15) * x * res[k+3]
                                ) / (x2 * (k^2 + 7k + 12))
                            res[i] = ai
                        end
                    end
                end
            end
        end
    end
    # Compose the Taylor series for the Bessel function with that of z
    z_tmp = arb_series(deepcopy(z.poly))
    z_tmp[0] = base_ring(parent(z.poly))(0)

    return Nemo.compose(res, z_tmp, n)
end

"""
    legendre_p(ν::arb_series, μ::arb, z::arb[, n = length(ν)])

> Compute the Taylor series of the Legendre function with respect to
  the parameter ν.
"""
function legendre_p(ν::arb_series, μ::arb, z::arb, n = length(ν))
    CC = ComplexField(z.parent.prec)
    PP = AcbPolyRing(CC, :x)
    νν = PP(ν.poly)
    μμ = PP(μ)
    zz = PP(z)

    a = -νν
    b = νν + 1
    c = 1 - μμ

    res = PP()

    ccall(
        (:acb_hypgeom_2f1_series_direct, :libarb),
        Cvoid,
        (
            Ref{acb_poly},
            Ref{acb_poly},
            Ref{acb_poly},
            Ref{acb_poly},
            Ref{acb_poly},
            Cint,
            Clong,
            Clong,
        ),
        res,
        a,
        b,
        c,
        0.5 * (1 - zz),
        1,
        2,
        CC.prec,
    )

    realres = arb_series(ν.poly.parent([real(coeff(res, i)) for i = 0:degree(res)]), n)
    ((1 + z) / (1 - z))^(μ / 2) * realres
end

"""
    legendre_p(ν::arb, μ::arb, z::arb_series[, n = length(z)])

> Compute the Taylor series of the Legendre function.

It's computed using a recursive formula for the Taylor coefficients
of the Legendre functions and then composing it with `z`.
"""
function legendre_p(ν::arb, μ::arb, z::arb_series, n = length(z))
    # TODO: Precompute values used several times
    res = arb_series(parent(z.poly)(), n)

    if n > 0
        x = z[0]
        a0 = legendre_p_safe(ν, μ, x)
        res[0] = a0

        if !isfinite(a0)
            return arb_series(parent(z.poly)(fill(NaN, n)), n)
        end

        if n > 1
            x2 = x^2
            onemx2 = 1 - x2

            a1 = ((ν + 1) * x * a0 - (ν - μ + 1) * legendre_p_safe(ν + 1, μ, x)) / (onemx2)
            res[1] = a1

            if !isfinite(a1)
                return arb_series(parent(z.poly)(fill(NaN, n)), n)
            end

            if n > 2
                ν2 = ν^2
                μ2 = μ^2

                a2 = (2x * a1 - (ν2 + ν - μ2 / (onemx2)) * a0) / (2(onemx2))
                res[2] = a2

                if n > 3
                    a3 =
                        (
                            4(onemx2)^2 * x * a2 +
                            2(2μ2 - (ν2 + ν) * (onemx2)) * x * a0 +
                            ((ν2 + ν + 2) * x2 + μ2 - ν2 - ν + 2) * (onemx2) * a1
                        ) / (6 * (onemx2)^3)
                    res[3] = a3

                    if n > 4
                        onemx22 = onemx2^2
                        xonemx2 = x * onemx2
                        ν2pν = ν2 + ν
                        μ2mν2pν = μ2 - ν2pν
                        fourx = 4x
                        for i = 4:n-1
                            k = i - 4
                            ai =
                                (
                                    (k + 1 + ν) * (ν - k) * res[k] -
                                    (0.5 * (2k^2 + (5k + 3) - ν2pν)) * fourx * res[k+1] +
                                    (
                                        ((-6k^2 - 24k - 24) + ν2pν) * x2 +
                                        (2k^2 + 8k + 8) +
                                        μ2mν2pν
                                    ) * res[k+2] +
                                    (2(k + 3) * (2k + 5)) * xonemx2 * res[k+3]
                                ) / (onemx22 * ((k + 4) * (k + 3)))
                            res[i] = ai
                        end
                    end
                end
            end
        end
    end
    # Compose the Taylor series for the Legendre function with that of
    # z
    z_tmp = arb_series(deepcopy(z.poly))
    z_tmp[0] = base_ring(parent(z.poly))(0)

    return Nemo.compose(res, z_tmp, n)
end

function legendre_p_safe(ν::arb, μ::arb, z::arb_series, n = length(z))
    legendre_p(ν, μ, z, n)
end


text:
using Documenter
using Latexify

makedocs(
    modules = [Latexify],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Latexify.jl",
    pages = [
        "index.md",
        # "Functions" => [
        #     "tutorials/latexify.md",
        #     "tutorials/latexinline.md",
        #     "tutorials/latexalign.md",
        #     "tutorials/latexarray.md",
        #     "tutorials/latextabular.md"
        # ],
        "tutorials/recipes.md",
        "Use with other packages" => [
            "tutorials/parameterizedfunctions.md",
            "tutorials/DiffEqBiological.md"
        ],
        "arguments.md",
        "tutorials/inner_workings.md",
    ],
    doctest = false
)

deploydocs(
    #deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/korsbo/Latexify.jl.git",
    target = "build",
    make = nothing,
    deps = nothing
    )

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
function func1big(x)
    return x^2
end

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD
using Flux
using DelimitedFiles
using Plots

function Afunc(system::AbstractSystem,mm::MixedModel,ind::Int64,dt::Float64)
  lam::Float64 = 2.0e3
  pot::Float64 = potential(system,mm.potentials[1])
  ind==1500 && system.x[2]>0.7 && pot<-145.0 ? lam : 0.0
end

### define system ###
dt = 0.0001
steps = 1500
tf = steps*dt

# mueller-brown potential
mb = MullerBrown()
# neural network model
nn = NNModel(2,Chain(Dense(3,10),Dense(10,2)))

# mixed model
model = MixedModel(2,[mb,nn])

# thermostat and system
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(2,0.0,[0.63,0.03],zeros(Float64,2),langevin)

# define value baseline function
Mx = 21
My = 21
Mt = 21
shape = (Mx,My,Mt)
vmux = zeros(Float64,Mx)
vmuy = zeros(Float64,My)
vmut = zeros(Float64,Mt)
for j=1:Mx
  vmux[j] = -1.5+(j-1)*(1.5-(-1.5))/(Mx-1)
end
for j=1:My
  vmuy[j] = -0.5+(j-1)*(2-(-0.5))/(My-1)
end
for j in 1:Mt
  vmut[j] = (j*steps*1.0/(Mt-1)) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  vmut[j] *= dt
end
vsigx = zeros(Float64,Mx)
vsigy = zeros(Float64,My)
vsigt = zeros(Float64,Mt)
vsigx[:] .= ((1.5-(-1.5))/(Mx-1)/2.0)^2 #variance
vsigy[:] .= ((2-(-0.5))/(My-1)/2.0)^2 #variance
vsigt[:] .= ((steps*1.0/(Mt-1))/2.0)^2 #variance
vsigt[:] .*= dt^2
# read in optimal coefficients
v = readdlm("coeff0V.txt")
vcs = zeros(Float64,shape)
for i in 1:21
  vcs[:,:,i] .= reshape(v[:,i],(21,21))'
end
vcs = reshape(vcs, (21*21*21))
vbl = GaussianValueBaseline2D(Mx,My,Mt,vmux,vmuy,vmut,vsigx,vsigy,vsigt,vcs)

# define callback function
mcvb = MCVBCallback(Afunc,nn,vbl)
# set up integrator and run!
integrator = StochasticEuler(system,model)
runtraj!(steps,dt,integrator,mcvb)
println(mcvb.dkl)
writedlm("gradf.txt",mcvb.gradf)
writedlm("gradv.txt",mcvb.gradv*dt)


text:
# # Multi-Class Classification

# We load all the necessary packages
using AbstractGPs
using ApproximateGPs
using ArraysOfArrays
using AugmentedGPLikelihoods
using Distributions
using LinearAlgebra
using SplitApplyCombine
# Plotting libraries
using Plots

# We create some random data (sorted for plotting reasons)
N = 100
Nclass = 3
x = collect(range(-10, 10; length=N));

# We will present simultaneously the bijective and non-bijective likelihood
liks = [
    CategoricalLikelihood(BijectiveSimplexLink(LogisticSoftMaxLink(Nclass))),
    CategoricalLikelihood(LogisticSoftMaxLink(Nclass)),
];

# This is small hack until https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl/pull/68 is merged
AugmentedGPLikelihoods.nlatent(::CategoricalLikelihood{<:BijectiveSimplexLink}) = Nclass - 1
AugmentedGPLikelihoods.nlatent(::CategoricalLikelihood{<:LogisticSoftMaxLink}) = Nclass

SplitApplyCombine.invert(x::ArrayOfSimilarArrays) = nestedview(flatview(x)')
# We define the models
kernel = 5.0 * with_lengthscale(SqExponentialKernel(), 2.0)
gp = GP(kernel)
fz = gp(x, 1e-8);
# We use a multi-output GP to generate our (N-1) latent GPs
gpm = GP(IndependentMOKernel(kernel))
X = MOInput(x, Nclass - 1);

# We sample (N-1) latent GPs to force the setting of the bijective version
# which the non-bijective version should also be able to recover
fs = rand(gpm(X, 1e-6))
fs = nestedview(reduce(hcat, Iterators.partition(fs, N)))
lik_true = liks[1](invert(fs)) # the likelihood
y = rand(lik_true);

# We build the one-hot encoding for each likelihood (different)
Ys = map(liks) do lik
    Y = nestedview(sort(unique(y))[1:nlatent(lik)] .== permutedims(y))
    return Y
end;
# We plot the sampled data
plts = map(["Bijective Logistic-softmax", "Logistic-softmax"]) do title
    plt = plot(; title)
    scatter!(plt, x, y; group=y, label=[1 2 3 4], msw=0.0)
    plot!(plt, x, vcat(fs, [zeros(N)]); color=[1 2 3 4], label="", lw=3.0)
end
plot(plts[1]; title="")
# ## CAVI Updates
# We write our CAVI algorithmm
function u_posterior(fz, m, S)
    return posterior(SparseVariationalApproximation(Centered(), fz, MvNormal(copy(m), S)))
end

function cavi!(fz::AbstractGPs.FiniteGP, lik, x, Y, ms, Ss, qΩ; niter=10)
    K = ApproximateGPs._chol_cov(fz)
    for _ in 1:niter
        posts_u = u_posterior.(Ref(fz), ms, Ss)
        posts_fs = marginals.([p_u(x) for p_u in posts_u])
        aux_posterior!(qΩ, lik, Y, SplitApplyCombine.invert(posts_fs))
        Ss .=
            inv.(
                Symmetric.(Ref(inv(K)) .+ Diagonal.(expected_auglik_precision(lik, qΩ, Y)))
            )
        ms .= Ss .* (expected_auglik_potential(lik, qΩ, Y) .- Ref(K \ mean(fz)))
    end
    return ms, Ss
end
# Now we just initialize the variational parameters and run CAVI
ms_Ss = map(liks, Ys) do lik, Y
    m = nestedview(zeros(N, nlatent(lik)))
    S = [Matrix{Float64}(I(N)) for _ in 1:nlatent(lik)]
    qΩ = init_aux_posterior(lik, N)
    fz = gp(x, 1e-8)
    cavi!(fz, lik, x, Y, m, S, qΩ; niter=20)
    return (; m, S)
end
# And visualize the obtained variational posterior
x_te = -10:0.01:10
for i in 1:2
    for j in 1:nlatent(liks[i])
        plot!(
            plts[i],
            x_te,
            u_posterior(fz, ms_Ss[i].m[j], ms_Ss[i].S[j]);
            color=j,
            alpha=0.3,
            lw=3.0,
            label="",
        )
    end
end
plot(plts...)

p_plts = [plot() for _ in 1:2]
for i in 1:2
    scatter!(p_plts[i], x, y / Nclass; group=y, label=[1 2 3 4], msw=0.0)

    ## vline!(p_plts[i], x, group=fs_ys[i].y, lw=20/length(x) * 20.0, alpha=0.2, ylims=(0,1),title="p(y=k|f)", label="")
    lik_pred =
        liks[i].(
            invert(
                mean.([
                    u_post(x_te) for u_post in u_posterior.(Ref(fz), ms_Ss[i].m, ms_Ss[i].S)
                ]),
            ),
        )
    ps = getproperty.(lik_pred, :p)
    lik_pred_x =
        liks[i].(
            invert(
                mean.([
                    u_post(x) for u_post in u_posterior.(Ref(fz), ms_Ss[i].m, ms_Ss[i].S)
                ]),
            ),
        )
    ps_x = getproperty.(lik_pred_x, :p)
    ps_true = getproperty.(lik_true.v, :p)
    for k in 1:Nclass
        plot!(p_plts[i], x, invert(ps_true)[k]; color=k, lw=2.0, label="")
        plot!(p_plts[i], x_te, invert(ps)[k]; color=k, lw=2.0, label="", ls=:dash)
    end
end
plot(p_plts...)
# ## ELBO
# How can one compute the Augmented ELBO?
# Again AugmentedGPLikelihoods provides helper functions
# to not have to compute everything yourself
function aug_elbo(lik, u_post, x, y)
    qf = marginals(u_post(x))
    qΩ = aux_posterior(lik, y, qf)
    return expected_logtilt(lik, qΩ, y, qf) - aux_kldivergence(lik, qΩ, y) -
           kldivergence(u_post.approx.q, u_post.approx.fz)
end

# aug_elbo(lik, u_posterior(fz, m, S), x, y)
# ## Gibbs Sampling
# We create our Gibbs sampling algorithm (we could do something fancier with
# AbstractMCMC)
function gibbs_sample(fz, lik, Y, fs, Ω; nsamples=200)
    K = ApproximateGPs._chol_cov(fz)
    Σ = [zeros(N, N) for _ in 1:nlatent(lik)]
    μ = [zeros(N) for _ in 1:nlatent(lik)]
    return map(1:nsamples) do _
        aux_sample!(Ω, lik, Y, invert(fs))
        Σ .= inv.(Symmetric.(Ref(inv(K)) .+ Diagonal.(auglik_precision(lik, Ω, Y))))
        μ .= Σ .* (auglik_potential(lik, Ω, Y) .- Ref(K \ mean(fz)))
        rand!.(MvNormal.(μ, Σ), fs)
        return copy(fs)
    end
end;
# We initialize our random variables
samples = map(liks, Ys, plts) do lik, Y, plt
    fs_init = nestedview(randn(N, nlatent(lik)))
    Ω = init_aux_variables(lik, N)
    ## Run the sampling for default number of iterations (200)
    return gibbs_sample(fz, lik, Y, fs_init, Ω)
    ## And visualize the samples overlapped to the variational posterior
    ## that we found earlier.
end;

for i in 1:2
    for fs in samples[i]
        for j in 1:nlatent(liks[i])
            plot!(plts[i], x, fs[j]; color=j, alpha=0.07, label="")
        end
    end
end
plot(plts...)


text:
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: mul!
using Combinatorics: permutations
using LinearAlgebra: rank, svd, Diagonal

"""
`struct SymmetricBasis`

### Constructors

Option 1: pass a `OneParticleBasis`
```julia
SymmetricBasis(φ, symgrp, basis1p, Bsel)
SymmetricBasis(φ, basis1p, Bsel)   # uses default symgrp = O3()
```
will first construct a `PIBasis` from these inputs and then call the second
constructor.

Option 1: pass a `PIBasis`
```julia
SymmetricBasis(φ, symgrp, pibasis)
SymmetricBasis(φ, pibasis)
```
If the PIbasis is already available, this directly constructs a
resulting SymmetricBasis; all possible permutation-invariant basis functions
will be symmetrised and then reduced to a basis (rather than spanning set)
"""
struct SymmetricBasis{BOP, PROP, SYM, REAL, VPROP} <: ACEBasis
   pibasis::PIBasis{BOP}
   A2Bmap::SparseMatrixCSC{PROP, Int}
   symgrp::SYM
   real::REAL
   B_pool::VectorPool{VPROP}
end

Base.length(basis::SymmetricBasis{BOP, PROP}) where {BOP, PROP} =
      size(basis.A2Bmap, 1)

valtype(basis::SymmetricBasis{BOP, PROP}) where {BOP, PROP} = basis.real(PROP)

# TODO: this is not nice, there should be proper promotion
valtype(basis::SymmetricBasis{BOP, PROP}, X::AbstractState) where {BOP, PROP} =
      valtype(basis)

gradtype(basis::SymmetricBasis, X::AbstractState) = gradtype(basis, typeof(X))

function gradtype(basis::SymmetricBasis, cfgorX)
   φ = zero(eltype(basis.A2Bmap))
   dAA = zero(gradtype(basis.pibasis, cfgorX))
   # note up to 0.12.4, basis.real used to be replaced with _real1234
   # a weird hacky thing. not sure why it now works with basis.real
   return typeof( basis.real( coco_o_daa(φ, dAA) ))
end

# -------- FIO

==(B1::SymmetricBasis, B2::SymmetricBasis) =
      ( (B1.pibasis == B2.pibasis) &&
        (B1.A2Bmap == B2.A2Bmap) &&
        (B1.real == B2.real) )

write_dict(B::SymmetricBasis{BOP, PROP}) where {BOP, PROP} =
      Dict( "__id__" => "ACE_SymmetricBasis",
            "pibasis" => write_dict(B.pibasis),
            "A2Bmap" => write_dict(B.A2Bmap),
            "symgrp" => write_dict(B.symgrp),
            "isreal" => (B.real == Base.real) )

read_dict(::Val{:ACE_SymmetricBasis}, D::Dict) =
      SymmetricBasis(read_dict(D["pibasis"]),
                     read_dict(D["A2Bmap"]),
                     read_dict(D["symgrp"]),
                     (D["isreal"] ? Base.real : Base.identity) )
# --------

SymmetricBasis(φ::AbstractProperty,
               basis1p::OneParticleBasis,
               Bsel::AbstractBasisSelector;
               kwargs...) =
      SymmetricBasis(φ, basis1p, O3(), Bsel; kwargs...)

SymmetricBasis(φ::AbstractProperty, pibasis; kwargs...) =
      SymmetricBasis(φ, O3(), pibasis; kwargs...)

SymmetricBasis(φ::AbstractProperty,
               basis1p::OneParticleBasis,
               symgrp::SymmetryGroup,
               Bsel::AbstractBasisSelector;
               isreal=isrealB(φ), kwargs...) =
      SymmetricBasis(φ, symgrp,
                     PIBasis(basis1p, symgrp, Bsel;
                             isreal = isrealAA(φ), kwargs..., property = φ);
                     isreal=isreal)

SymmetricBasis(φ::AbstractProperty, symgrp::SymmetryGroup, pibasis::PIBasis;
               isreal=false) =
      SymmetricBasis(φ, symgrp, pibasis, isreal ? Base.real : Base.identity)

function SymmetricBasis(φ::TP, symgrp::SymmetryGroup, pibasis::PIBasis,
                        _real) where {TP <: AbstractProperty}

   # AA index -> AA spec
   AAspec = get_spec(pibasis)

   # construct the reverse mapping AA spec -> iAA and A spec -> iA
   invAAspec = Dict{Any, Int}()
   for (iAA, AA) in enumerate(AAspec)
      invAAspec[AA] = iAA
   end
   invAspec = Dict{Any, Int}()
   for (iA, A) in enumerate(get_spec(pibasis.basis1p))
      invAspec[A] = iA
   end

   # allocate the datastructure that computes and caches the
   # coupling coefficients
   # TODO: should this be stored with the basis?
   #       or maybe written to a file on disk? and then flushed every time
   #       we finish with a basis construction???
   # TODO: for sure this needs to become a function of the symmetry group?
   rotc = Rot3DCoeffs(φ, real(valtype(pibasis)))
   # allocate triplet format
   TCC = coco_type(TP)
   Irow, Jcol, vals = Int[], Int[], TCC[]
   # count the number of PI basis functions = number of rows
   idxB = 0

   # loop through AA basis, but skip most of them ...
   for (iAA, AA) in enumerate(AAspec)
      # determine whether we need to compute coupling coefficients for this
      # basis function or whether it will be included in a different
      # coco computation?
      if !is_refbasisfcn(symgrp, AA)
         continue
      end
      # compute the cocos
      U, AAcols = coupling_coeffs(symgrp, AA, rotc)

      # loop over the rows of U -> each specifies a basis function which we now
      # need to incorporate into the basis
      for irow = 1:size(U, 1)
         idxB += 1
         # loop over the columns of U / over brows
         for (icol, bcol) in enumerate(AAcols)
            # put bcol into the correct order
            bcol_ordered = _get_ordered(bcol, invAspec)
            # this is a subtle step: bcol and bcol_ordered are equivalent
            # permutation-invariant basis functions. This means we will
            # add the same PI basis function several times, but in the call to
            # `sparse` the values will just be added.
            if !haskey(invAAspec, bcol_ordered)
               @show bcol_ordered
               @show degree(bcol_ordered, NaiveTotalDegree(), pibasis.basis1p)
            end
            idxAA = invAAspec[bcol_ordered]
            push!(Irow, idxB)
            push!(Jcol, idxAA)
            push!(vals, U[irow, icol])
         end
      end
   end
   # TODO: filter and throw out everything that hasn't been used!!
   # create CSC: [   triplet    ]  nrows   ncols
   A2Bmap = sparse(Irow, Jcol, vals, idxB, length(AAspec))
   return SymmetricBasis(pibasis, A2Bmap, symgrp, _real)
end

function SymmetricBasis(pibasis, A2Bmap, symgrp, _real)
   PROP = _real(eltype(A2Bmap))
   B_pool = VectorPool{PROP}()
   return SymmetricBasis(pibasis, A2Bmap, symgrp, _real, B_pool)
end

"""
produce the ordered tuple defining the AA basis function uniquely.
"""
function _get_ordered(bb, invAspec)
   iAs = [invAspec[b] for b in bb]
   return bb[ sortperm(iAs, rev = true) ]
end

# ------------------- exporting the basis spec

# this doesn't provide the "full" specification, just collects
# the n and l but not the m or coupling coefficients.

function get_spec(basis::SymmetricBasis)
   spec_AA = get_spec(basis.pibasis)
   spec_B = []
   for iB = 1:length(basis)
      iAA = findfirst(norm.(basis.A2Bmap[iB, :]) .!= 0)
      push!(spec_B, get_sym_spec(basis.symgrp, spec_AA[iAA]))
   end
   return identity.(spec_B)
end

# ---------------- A modified sparse matmul

# TODO: move this stuff all to aux?

using SparseArrays: AbstractSparseMatrixCSC,
				        nonzeros, rowvals, nzrange

using LinearAlgebra: Transpose

function genmul!(C, A::AbstractSparseMatrixCSC, B, mulop)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    fill!(C, zero(eltype(C)))
    for k in 1:size(C, 2)
        @inbounds for col in 1:size(A, 2)
            αxj = B[col,k]
            for j in nzrange(A, col)
                C[rv[j], k] += mulop(nzv[j], αxj)
            end
        end
    end
    return C
end

function genmul!(C, xA::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B, mulop)
   A = xA.parent
   size(A, 2) == size(C, 1) || throw(DimensionMismatch())
   size(A, 1) == size(B, 1) || throw(DimensionMismatch())
   size(B, 2) == size(C, 2) || throw(DimensionMismatch())
   nzv = nonzeros(A)
   rv = rowvals(A)
   fill!(C, zero(eltype(C)))
   for k in 1:size(C, 2)
       @inbounds for col in 1:size(A, 2)
           tmp = zero(eltype(C))
           for j in nzrange(A, col)
               tmp += mulop(nzv[j], B[rv[j],k])
           end
           C[col,k] += tmp
       end
   end
   return C
end

#dispatching for SVectors and B being a list of matrices rather than a single matrix
#the function below does the same, but return identical coppies for all properties
#only works on things that are parameter independent.
function adjointgenmul!(C::AbstractVector{<: SVector}, A::AbstractSparseMatrixCSC, B, mulop)
   for prop in 1:length(C[1]) #TODO is it worth checking everything?
      size(A, 2) == size(B[prop], 1) || throw(DimensionMismatch())
      size(B[prop], 2) == size(C, 2) || throw(DimensionMismatch())
   end
   size(A, 1) == size(C, 1) || throw(DimensionMismatch())
   nzv = nonzeros(A)
   rv = rowvals(A)
   fill!(C, zero(eltype(C)))
   for k in 1:size(C, 2)
      for prop in 1:length(C[1]) #we add this loop over properties
         @inbounds for col in 1:size(A, 2)
            αxj = B[prop][col,k]
            for j in nzrange(A, col)
                  mop = mulop(nzv[j], αxj) #find the value we need
                  zerotmp = zeros(length(C[1])) #make an array of zeros
                  zerotmp[prop] = 1 #put a one only in the current position
                  tmp = SVector{length(C[1])}(zerotmp) #make it an SVector
                  C[rv[j], k] += mop * tmp #this should add only on the current property
            end
         end
      end
   end
   return C
end

#TODO, is this worth having?
#dispatching for SVectors
#here we simply pass a coppy or a fill() or the mulop to every
#property on the SVector.
function genmul!(C::AbstractVector{<: SVector}, A::AbstractSparseMatrixCSC, B, mulop)
   size(A, 2) == size(B, 1) || throw(DimensionMismatch())
   size(A, 1) == size(C, 1) || throw(DimensionMismatch())
   size(B, 2) == size(C, 2) || throw(DimensionMismatch())
   nzv = nonzeros(A)
   rv = rowvals(A)
   fill!(C, zero(eltype(C)))
   for k in 1:size(C, 2)
       @inbounds for col in 1:size(A, 2)
           αxj = B[col,k]
           for j in nzrange(A, col)
               mop = mulop(nzv[j], αxj)
               C[rv[j], k] += mop * ones(SVector{length(C[1]),eltype(mop)})
           end
       end
   end
   return C
end

# ---------------- Evaluation code

function evaluate!(B, basis::SymmetricBasis, cfg::AbstractConfiguration)
   # compute AA
   AA = acquire_B!(basis.pibasis, cfg)
   evaluate!(AA, basis.pibasis, cfg)
   return evaluate!(B, basis, AA)
end

# this function allows us to attach multiple symmetric bases to a single
#     𝑨 basis
#  TODO: this is extremely inefficient for multiple species But it is simple
#        and clean and will do for now...
function evaluate!(B, basis::SymmetricBasis, AA::AbstractVector{<: Number})
   genmul!(B, basis.A2Bmap, AA, (a, b) -> basis.real(a * b))
end

# ---------------- gradients

function evaluate_d!(dB, basis::SymmetricBasis, cfg::AbstractConfiguration,
                     args...)   # args... could be nothing or sym
   # compute AA
   AA = acquire_B!(basis.pibasis, cfg)
   dAA = acquire_dB!(basis.pibasis, cfg)
   evaluate_ed!(AA, dAA, basis.pibasis, cfg, args...)
   evaluate_d!(dB, basis, AA, dAA)
   return dB
end

function evaluate_d!(dB, basis::SymmetricBasis,
                     AA::AbstractVector{<: Number}, dAA)
   # note up to 0.12.4, basis.real used to be replaced with _real1234
   # a weird hacky thing. not sure why it now works with basis.real
   genmul!(dB, basis.A2Bmap, dAA,
           (a, b) -> basis.real( ACE.coco_o_daa(a, b) ) )
end

function scaling(basis::SymmetricBasis, p)
   wwpi = scaling(basis.pibasis, p)
   wwrpi = abs2.(norm.(basis.A2Bmap)) * abs2.(wwpi)
   return sqrt.(wwrpi)
end


text:
"""
Inversion variables
"""
mutable struct X{T}
	x::Vector{T}
	last_x::Vector{T}
	lower_x::Vector{T}
	upper_x::Vector{T}
	precon::Vector{T}
	preconI::Vector{T}
	weights::Vector{T}
end

function X(n,T)
	return X{T}(zeros(n), randn(n), fill(-Inf,n), fill(Inf,n), ones(n), ones(n), zeros(n))
end


text:
function k(α,β,γ)

			x̂ = cartesian(α)
			ξ = cartesian(β)[1]
			η = cartesian(γ)[1]

			ŷ = [ξ*(1-η), η]

			#jacobian of y--->ŷ = (1-y[2])

			chartI = simplex(p1,x̂,p0)
			chartII = simplex(p2,x̂,p1)
			chartIII = simplex(p0,x̂,p2)

			n1 = neighborhood(chartI, ŷ)
			n2 = neighborhood(chartII, ŷ)
			n3 = neighborhood(chartIII, ŷ)

			yI = cartesian(n1)
			yII = cartesian(n2)
			yIII = cartesian(n3)

			return(Kernel(x̂,yI)*jacobian(n1)*(1-η) +
					Kernel(x̂,yII)*jacobian(n2)*(1-η) +
					Kernel(x̂,yIII)*jacobian(n3)*(1-η))
end

function verifintegral1(sourcechart::CompScienceMeshes.Simplex{3,2,1,3,Float64},
	 testchart::CompScienceMeshes.Simplex{3,2,1,3,Float64}, integrand, accuracy::Int64)

	global Kernel, p0, p1, p2

	Kernel = integrand

	p0 = Sourcechart.vertices[1]
	p1 = Sourcechart.vertices[2]
	p2 = Sourcechart.vertices[3]

	qps1 = quadpoints(Sourcechart, accuracy)

	path = simplex(point(0), point(1))
	qps2 = quadpoints(path, accuracy)

	result = sum(w*w2*w1*k(α,β,γ) for (β,w1) in qps2, (γ,w2) in qps2, (α,w) in qps1)

	return (result)

end

function verifintegral2(sourcechart::CompScienceMeshes.Simplex{3,2,1,3,Float64},
	 testchart::CompScienceMeshes.Simplex{3,2,1,3,Float64}, integrand, accuracy::Int64)

	qps1 = quadpoints(sourcechart, accuracy)
	qps2 = quadpoints(testchart, accuracy)

	result = sum(w2*w1*integrand(cartesian(x),cartesian(y))
				for (x,w2) in qps2, (y,w1) in qps1)

	return (result)

end


text:
export CollectGems

mutable struct CollectGems{R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Gem}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    num_gem_init::Int
    num_gem_current::Int
    gem_reward::Float64
    reward::Float64
    rng::R
end

function CollectGems(;n=8, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT, rng=Random.GLOBAL_RNG)
    objects = (EMPTY, WALL, GEM)
    w = GridWorldBase(objects, n, n)

    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true

    num_gem_init = n - 1
    num_gem_current = num_gem_init

    gem_reward = 1.0
    reward = 0.0

    env = CollectGems(w, agent_start_pos, Agent(dir=agent_start_dir), num_gem_init, num_gem_current, gem_reward, reward, rng)

    reset!(env, agent_start_pos = agent_start_pos, agent_start_dir = agent_start_dir)

    return env

end

function (env::CollectGems)(::MoveForward)
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)
    env.reward = 0.0
    if !env.world[WALL, dest]
        env.agent_pos = dest
        if env.world[GEM, dest]
            env.world[GEM, dest] = false
            env.world[EMPTY, dest] = true
            env.num_gem_current = env.num_gem_current - 1
            env.reward = env.gem_reward
        end
    end
    env
end

RLBase.get_terminal(env::CollectGems) = env.num_gem_current <= 0

function RLBase.reset!(env::CollectGems; agent_start_pos = CartesianIndex(2, 2), agent_start_dir = RIGHT)

    n = size(env.world)[end]
    env.world[EMPTY, 2:n-1, 2:n-1] .= true
    env.world[GEM, 1:n, 1:n] .= false
    env.num_gem_current = env.num_gem_init

    env.reward = 0.0
    env.agent_pos = agent_start_pos
    agent = get_agent(env)
    set_dir!(agent, agent_start_dir)

    gem_placed = 0
    while gem_placed < env.num_gem_init
        gem_pos = CartesianIndex(rand(env.rng, 2:n-1), rand(env.rng, 2:n-1))
        if (gem_pos == env.agent_pos) || (env.world[GEM, gem_pos] == true)
            continue
        else
            env.world[GEM, gem_pos] = true
            env.world[EMPTY, gem_pos] = false
            gem_placed = gem_placed + 1
        end
    end

    return env
end


text:
@with_kw mutable struct Sobol <: GSAMethod
    order::Array{Int}=[0, 1]
    nboot::Int=1
    conf_int::Float64=0.95
end

mutable struct SobolResult{T1, T2, T3, T4}
    S1::T1
    S1_Conf_Int::T2
    S2::T3
    S2_Conf_Int::T4
    ST::T1
    ST_Conf_Int::T2
end

function fuse_designs(A, B)
    d = size(A,1)
    Aᵦ = [copy(A) for i in 1:d]
    for i in 1:d
        Aᵦ[i][i,:] = B[i,:]
    end
    hcat(A,B,reduce(hcat,Aᵦ))
end

function gsa(f, method::Sobol, A::AbstractMatrix, B::AbstractMatrix;
             batch=false, Ei_estimator = :Jansen1999, kwargs...)

    d,n = size(A)
    n = Int(n/method.nboot)
    multioutput = false
    desol = false
    all_points = mapreduce(fuse_designs, hcat, [A[:,n*(i-1)+1:n*(i)] for i in 1:method.nboot], [B[:,n*(i-1)+1:n*(i)] for i in 1:method.nboot])

    if batch
        all_y = f(all_points)
        multioutput = all_y isa AbstractMatrix
    else
        _y = [f(all_points[:, i]) for i in 1:size(all_points, 2)]
        multioutput = !(eltype(_y) <: Number)
        if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
            y_size = size(_y[1])
            _y = vec.(_y)
            desol = true
        end
        all_y = multioutput ? reduce(hcat, _y) : _y
    end

    Eys = multioutput ? Array{eltype(all_y)}[] : eltype(all_y)[]
    Varys = multioutput ? Array{eltype(all_y)}[] : eltype(all_y)[]
    Vᵢs = multioutput ? Array{eltype(all_y)}[] : Array{eltype(all_y)}[]
    Vᵢⱼs = multioutput ? Array{eltype(all_y)}[] : Array{eltype(all_y)}[]
    Eᵢs = multioutput ? Array{eltype(all_y)}[] : Array{eltype(all_y)}[]
    if !multioutput
        for i in 1:d+2:(d+2)*method.nboot
            push!(Eys,mean(all_y[(i-1)*n+1:(i+1)*n]))
            push!(Varys,var(all_y[(i-1)*n+1:(i+1)*n]))

            fA = all_y[(i-1)*n+1:i*n]
            fB = all_y[(i*n+1):(i+1)*n]
            fAⁱ= [all_y[(j*n+1):((j+1)*n)] for j in i+1:(i+d)]

            push!(Vᵢs,[sum(fB.*(fAⁱ[k].-fA)) for k in 1:d]./n)
            if 2 in method.order
                push!(Vᵢⱼs,[sum(abs2, fAⁱ[k] - fAⁱ[j]) for k in 1:d for j in k+1:d]./(2n))
            end
            if Ei_estimator == :Homma1996
                push!(Eᵢs,[Varys[i] .- sum(fA .* fAⁱ[k])./(n) + Eys[i].^2 for k in 1:d])
            elseif Ei_estimator == :Sobol2007
                push!(Eᵢs,[sum(abs2,fA-fAⁱ[k]) for k in 1:d]./(2n))
            elseif Ei_estimator == :Jansen1999
                push!(Eᵢs,[sum(fA.*(fA.-fAⁱ[k])) for k in 1:d]./(n))
            end
        end
    else
        for i in 1:d+2:(d+2)*method.nboot
            push!(Eys,mean(all_y[:, (i-1)*n+1:(i+1)*n],dims=2))
            push!(Varys,var(all_y[:, (i-1)*n+1:(i+1)*n],dims=2))

            fA = all_y[:, (i-1)*n+1:i*n]
            fB = all_y[:, (i*n+1):(i+1)*n]
            fAⁱ= [all_y[:, (j*n+1):((j+1)*n)] for j in i+1:(i+d)]

            push!(Vᵢs,reduce(hcat, [sum(fB.*(fAⁱ[k].-fA), dims=2) for k in 1:d]./n))

            if 2 in method.order
                push!(Vᵢⱼs,reduce(hcat, [sum(abs2, fAⁱ[k] - fAⁱ[j], dims=2) for k in 1:d for j in k+1:d]./(2n)))
            end
            if Ei_estimator == :Homma1996
                push!(Eᵢs,reduce(hcat, [Varys[i] .- sum(fA .* fAⁱ[k], dims=2)./(n) + Eys[i].^2 for k in 1:d]))
            elseif Ei_estimator == :Sobol2007
                push!(Eᵢs,reduce(hcat, [sum(abs2,fA-fAⁱ[k],dims=2) for k in 1:d]./(2n)))
            elseif Ei_estimator == :Jansen1999
                push!(Eᵢs,reduce(hcat, [sum(fA.*(fA.-fAⁱ[k]), dims=2) for k in 1:d]./(n)))
            end
        end
    end
    if 2 in method.order
        Sᵢⱼs = [(Vᵢⱼs[i])./Varys[i] for i in 1:method.nboot]
    end

    Sᵢs = [Vᵢs[i] ./Varys[i] for i in 1:method.nboot]
    Tᵢs = [Eᵢs[i] ./Varys[i] for i in 1:method.nboot]
    if method.nboot > 1
        size_ = size(Sᵢs[1])
        S1 = [[Sᵢ[i] for Sᵢ in Sᵢs] for i in 1:length(Sᵢs[1])]
        ST = [[Tᵢ[i] for Tᵢ in Tᵢs] for i in 1:length(Tᵢs[1])]

        function calc_ci(x)
            alpha = (1 - method.conf_int)
            tstar = quantile(TDist(length(x)-1), 1 - alpha/2)
            std(x)/sqrt(length(x))
        end
        S1_CI = map(calc_ci,S1)
        ST_CI = map(calc_ci,ST)

        if 2 in method.order
            size__= size(Sᵢⱼs[1])
            S2 = [[Sᵢⱼ[i] for Sᵢⱼ in Sᵢⱼs] for i in 1:length(Sᵢⱼs[1])]
            S2_CI = reshape(map(calc_ci,S2),size__...)
            Sᵢⱼ = reshape(mean.(S2),size__...)
        end

        Sᵢ = reshape(mean.(S1),size_...)
        Tᵢ = reshape(mean.(ST),size_...)
    else
        Sᵢ = Sᵢs[1]
        Tᵢ = Tᵢs[1]
        if 2 in method.order
            Sᵢⱼ = Sᵢⱼs[1]
        end
    end
    if desol
        f_shape = x -> [reshape(x[:,i],y_size) for i in 1:size(x,2)]
        Sᵢ = f_shape(Sᵢ)
        if 2 in method.order
            Sᵢⱼ = f_shape(Sᵢⱼ)
        end
        Tᵢ = f_shape(Tᵢ)
    end
    return SobolResult(Sᵢ,
                     method.nboot > 1 ? reshape(S1_CI,size_...) : nothing,
                     2 in method.order ? Sᵢⱼ : nothing,
                     method.nboot > 1 && 2 in method.order ? S2_CI : nothing,
                     Tᵢ ,
                     method.nboot > 1 ? reshape(ST_CI,size_...) : nothing)
end

function gsa(f,method::Sobol,p_range::AbstractVector; N, kwargs...)
    A = QuasiMonteCarlo.generate_design_matrices(N, [i[1] for i in p_range], [i[2] for i in p_range], QuasiMonteCarlo.SobolSample(),2*method.nboot)
    gsa(f, method, hcat(A[1:method.nboot]...), hcat(A[method.nboot+1:end]...); kwargs...)
end


text:
@doc raw"""
```
t = find_threshold(Otsu(), histogram, edges)
```

Under the assumption that the histogram is bimodal the threshold is
set so that the resultant between-class variance is maximal.

# Output

Returns a real number `t` in `edges`. The `edges` parameter represents an
`AbstractRange` which specifies the intervals associated with the histogram bins.

# Details

Let ``f_i`` ``(i=1 \ldots I)`` denote the number of observations in the
``i``th bin of the histogram. Then the probability that an observation
belongs to the ``i``th bin is given by  ``p_i = \frac{f_i}{N}`` (``i = 1,
\ldots, I``), where ``N = \sum_{i=1}^{I}f_i``.

The choice of a threshold ``T`` partitions the data into two categories, ``C_0``
and ``C_1``. Let
```math
P_0(T) = \sum_{i = 1}^T p_i \quad \text{and} \quad P_1(T) = \sum_{i = T+1}^I p_i
```
denote the cumulative probabilities,
```math
\mu_0(T) = \sum_{i = 1}^T i \frac{p_i}{P_0(T)} \quad \text{and} \quad \mu_1(T) = \sum_{i = T+1}^I i \frac{p_i}{P_1(T)}
```
denote the means, and
```math
\sigma_0^2(T) = \sum_{i = 1}^T (i-\mu_0(T))^2 \frac{p_i}{P_0(T)} \quad \text{and} \quad \sigma_1^2(T) = \sum_{i = T+1}^I (i-\mu_1(T))^2 \frac{p_i}{P_1(T)}
```
denote the variances of categories ``C_0`` and ``C_1``, respectively.
Furthermore, let
```math
\mu = P_0(T)\mu_0(T) + P_1(T)\mu_1(T),
```
represent the overall mean,
```math
\sigma_b^2(T) = P_0(T)(\mu_0(T) - \mu)^2 + P_1(T)(\mu_1(T) - \mu)^2,
```
the between-category variance, and
```math
\sigma_w^2(T) = P_0(T) \sigma_0^2(T) +  P_1(T)\sigma_1^2(T)
```
the within-category variance, respectively.

Finding the discrete value ``T`` which maximises the function ``\sigma_b^2(T)``
produces the sought-after threshold value (i.e. the bin which determines the
threshold). As it turns out, that threshold value is equal to the threshold
decided by minimizing the within-category variances criterion ``\sigma_w^2(T)``.
Furthermore, that threshold is also the same as the threshold calculated by
maximizing the ratio of between-category variance to within-category variance.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Otsu(), counts[1:end], edges)
```

# Reference

1. Nobuyuki Otsu (1979). “A threshold selection method from gray-level histograms”. *IEEE Trans. Sys., Man., Cyber.* 9 (1): 62–66. [doi:10.1109/TSMC.1979.4310076](http://dx.doi.org/doi:10.1109/TSMC.1979.4310076)
"""
function find_threshold(algorithm::Otsu,  histogram::AbstractArray, edges::AbstractRange)
  N = sum(histogram)
  pdf = histogram / N
  first_bin = firstindex(pdf)
  last_bin = lastindex(pdf)
  cumulative_zeroth_moment = cumsum(pdf)
  cumulative_first_moment = cumsum(edges .* pdf)
  μ_T = cumulative_first_moment[end]
  σ²_T = sum( ((first_bin:last_bin) .- μ_T).^2  .* pdf )
  maxval = zero(eltype(first(pdf)))

  # Equation (6) for determining the probability of the first class.
  function ω(k)
    let cumulative_zeroth_moment = cumulative_zeroth_moment
      return cumulative_zeroth_moment[k]
    end
  end

  # Equation (7) for determining the mean of the first class.
  function μ(k)
   let cumulative_first_moment = cumulative_first_moment
     return cumulative_first_moment[k]
   end
  end

  # Equation (18) for determining the between-cass variance.
  function σ²(k)
    let μ_T = μ_T
      return (μ_T*ω(k) - μ(k))^2 / (ω(k) * (1-ω(k)))
    end
  end

  t = first_bin
  for k = first_bin:last_bin - 1
    val = σ²(k)
    if val  > maxval
       maxval = val
       t = k
    end
  end
  return edges[t]
end


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

# Tests for /base/stacktraces.jl

let
    @noinline child() = stacktrace()
    @noinline parent() = child()
    @noinline grandparent() = parent()
    line_numbers = @__LINE__() .- [3, 2, 1]
    stack = grandparent()

    # Basic tests.
    @assert length(stack) >= 3 "Compiler has unexpectedly inlined functions"

    @test [:child, :parent, :grandparent] == [f.func for f in stack[1:3]]
    for (line, frame) in zip(line_numbers, stack[1:3])
        @test [Symbol(@__FILE__), line] == [frame.file, frame.line]
    end
    @test [false, false, false] == [f.from_c for f in stack[1:3]]

    # Test remove_frames!
    stack = StackTraces.remove_frames!(grandparent(), :parent)
    @test stack[1] == StackFrame(:grandparent, @__FILE__, line_numbers[3])

    stack = StackTraces.remove_frames!(grandparent(), [:child, :something_nonexistent])
    @test stack[1:2] == [
        StackFrame(:parent, @__FILE__, line_numbers[2]),
        StackFrame(:grandparent, @__FILE__, line_numbers[3])
    ]

    b = PipeBuffer()
    frame = stack[1]
    serialize(b, frame)
    frame2 = deserialize(b)
    @test frame !== frame2
    @test frame == frame2
    @test frame.linfo !== nothing
    @test frame2.linfo === nothing
end

# Test from_c
let (default, with_c, without_c) = (stacktrace(), stacktrace(true), stacktrace(false))
    @test default == without_c
    @test length(with_c) > length(without_c)
    @test !isempty(filter(frame -> frame.from_c, with_c))
    @test isempty(filter(frame -> frame.from_c, without_c))
end

@test StackTraces.lookup(C_NULL) == [StackTraces.UNKNOWN]

let ct = current_task()
    # After a task switch, there should be nothing in catch_backtrace
    yieldto(@task yieldto(ct))
    @test catch_backtrace() == StackFrame[]

    @noinline bad_function() = throw(UndefVarError(:nonexistent))
    @noinline function try_stacktrace()
        try
            bad_function()
        catch
            return stacktrace()
        end
    end
    @noinline function try_catch()
        try
            bad_function()
        catch
            return catch_stacktrace()
        end
    end
    line_numbers = @__LINE__() .- [15, 10, 5]

    # Test try...catch with stacktrace
    @test try_stacktrace()[1] == StackFrame(:try_stacktrace, @__FILE__, line_numbers[2])

    # Test try...catch with catch_stacktrace
    @test try_catch()[1:2] == [
        StackFrame(:bad_function, @__FILE__, line_numbers[1]),
        StackFrame(:try_catch, @__FILE__, line_numbers[3])
    ]
end

module inlined_test
using Test
@inline g(x) = (y = throw("a"); y) # the inliner does not insert the proper markers when inlining a single expression
@inline h(x) = (y = g(x); y)       # this test could be extended to check for that if we switch to linear representation
f(x) = (y = h(x); y)
trace = (try; f(3); catch; catch_stacktrace(); end)[1:3]
can_inline = Bool(Base.JLOptions().can_inline)
for (frame, func, inlined) in zip(trace, [g,h,f], (can_inline, can_inline, false))
    @test frame.func === typeof(func).name.mt.name
    #@test get(frame.linfo).def === which(func, (Any,)).func
    #@test get(frame.linfo).specTypes === Tuple{typeof(func), Int}
    # line
    @test frame.file === Symbol(@__FILE__)
    @test !frame.from_c
    @test frame.inlined === inlined
end
end

let src = Meta.lower(Main, quote let x = 1 end end).args[1]::CodeInfo,
    li = ccall(:jl_new_method_instance_uninit, Ref{Core.MethodInstance}, ()),
    sf

    li.inferred = src
    li.specTypes = Tuple{}
    li.def = @__MODULE__
    sf = StackFrame(:a, :b, 3, li, false, false, 0)
    repr = string(sf)
    @test repr == "Toplevel MethodInstance thunk at b:3"
end
let li = typeof(fieldtype).name.mt.cache.func::Core.MethodInstance,
    sf = StackFrame(:a, :b, 3, li, false, false, 0),
    repr = string(sf)
    @test repr == "fieldtype(...) at b:3"
end

let ctestptr = cglobal((:ctest, "libccalltest")),
    ctest = StackTraces.lookup(ctestptr + 1)

    @test length(ctest) == 1
    @test ctest[1].func === :ctest
    @test ctest[1].linfo === nothing
    @test ctest[1].from_c
    @test ctest[1].pointer === UInt64(ctestptr)
end

# issue #19655
let st = stacktrace(empty!(backtrace()))
    # not in a `catch`, so should return an empty StackTrace
    @test isempty(st)
    @test isa(st, StackTrace)
end

module StackTracesTestMod
    unfiltered_stacktrace() = stacktrace()
    filtered_stacktrace() = StackTraces.remove_frames!(stacktrace(), StackTracesTestMod)
end

# Test that `removes_frames!` can correctly remove frames from within the module
trace = StackTracesTestMod.unfiltered_stacktrace()
@test contains(string(trace), "unfiltered_stacktrace")

trace = StackTracesTestMod.filtered_stacktrace()
@test !contains(string(trace), "filtered_stacktrace")

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using Documenter, Arena

makedocs(
    modules = [Arena],
    format = :html,
    sitename = "Arena.jl",
    pages = Any[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Functions" => "func_ref.md"
    ]
    # html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo   = "github.com/dehann/Arena.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    julia  = "0.6",
    osname = "linux"
)


text:
using Turf

box = [-95, 30 ,-85, 40]
cellside = 60

pg = point_grid(box, cellside) # a grid of points
rg = rectangle_grid(box, cellside, cellside / 2) # a grid of rectangles
tg = triangle_grid(box, cellside) # a grid of triangles
sg = square_grid(box, cellside) # a grid of squares

fc = [pg, rg, tg, sg]

for (i, fcoll) in enumerate(fc)
    open(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/grids/grid$(i).result.geojson", "w") do file
        write(file, geojson(fcoll))
    end
end


text:
using ExcelDataStreams
using Base.Test

# write your own tests here
include("ExcelDataStream.jl")


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

"""
    struct BridgeSampling <: IntegrationAlgorithm

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

BridgeSampling integration algorithm.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct BridgeSampling{TR<:AbstractDensityTransformTarget,ESS<:EffSampleSizeAlgorithm} <: IntegrationAlgorithm
    trafo::TR = PriorToGaussian()
    essalg::ESS = EffSampleSizeFromAC()
    strict::Bool = true
    # ToDo: add argument for proposal density generator
end
export BridgeSampling

function bat_integrate_impl(target::SampledDensity, algorithm::BridgeSampling)
    transformed_target, trafo = bat_transform(algorithm.trafo, target)
    density = unshaped(transformed_target.density)
    samples = unshaped.(transformed_target.samples)

    integral = bridge_sampling_integral(density, samples,algorithm.strict, algorithm.essalg )
    (result = integral,)
end

function bridge_sampling_integral(
    target_density::AbstractDensity,
    target_samples::DensitySampleVector,
    proposal_density::AbstractDensity,
    proposal_samples::DensitySampleVector,
    strict::Bool,
    ess_alg::EffSampleSizeAlgorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), target_samples)
    )

    N1 = Int(sum(target_samples.weight))
    N2 = Int(sum(proposal_samples.weight))

    #####################
    # Evaluate integral #
    #####################
    #calculate elements for iterative determination of marginal likelhood
    l1 = [exp(target_samples.logd[i]-logdensityof(proposal_density,x)) for (i,x) in enumerate(target_samples.v)]
    l2 = [exp(logdensityof(target_density,x)-proposal_samples.logd[i]) for (i,x) in enumerate(proposal_samples.v)]
    s1 = N1/(N2+N1)
    s2 = N2/(N1+N2)

    #calculate marginal likelhood iteratively
    prev_int = 0
    counter = 0
    current_int = 0.1
    while abs(current_int-prev_int)/current_int > 10^(-15)
        prev_int = current_int
        numerator = 0
        for (i, w) in enumerate(proposal_samples.weight)
            numerator += w*(l2[i]/(s1*l2[i]+s2*prev_int))
        end
        numerator = numerator/N2

        denominator = 0
        for (i, w) in enumerate(target_samples.weight)
            denominator += w/(s1*l1[i]+s2*prev_int)
        end
        denominator = denominator/N1

        current_int = numerator/denominator
        if counter == 500
            msg = "The iterative scheme is not converging!!"
            if strict
                throw(ErrorException(msg))
            else
                @warn(msg)
            end
        end
        counter=counter+1
    end

    #################
    #Evaluate error #
    #################
    #pre calculate objects for error estimate
    f1 = [exp(logdensityof(target_density,x))/current_int/(s1*exp(logdensityof(target_density,x))/current_int+s2*exp(proposal_samples.logd[i])) for (i,x) in enumerate(proposal_samples.v)]
    f2 = [[exp(logdensityof(proposal_density,x))/(s1*exp(target_samples.logd[i])/current_int+s2*exp(logdensityof(proposal_density,x)))] for (i,x) in enumerate(target_samples.v)]
    f2_density_vector = DensitySampleVector(f2,target_samples.logd,weight=target_samples.weight)

    mean1, var1 = StatsBase.mean_and_var(f1, FrequencyWeights(proposal_samples.weight))
    mean2, var2 = mean(f2_density_vector)[1],cov(f2_density_vector)[1]

    N1_eff = bat_eff_sample_size(f2_density_vector,ess_alg).result[1]
    # calculate  Root mean squared error
    r_MSE = sqrt(var1/(mean1^2*N2)+(var2/mean2^2)/N1_eff)*current_int

    integral = Measurements.measurement(current_int, r_MSE)

    return integral
end

function bridge_sampling_integral(
    target_density::AbstractDensity,
    target_samples::DensitySampleVector,
    strict::Bool,
    ess_alg::EffSampleSizeAlgorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), target_samples)
    )

    num_samples = size(target_samples.weight)[1]
    n_first = floor(Int,num_samples/2)
    first_batch = target_samples[1:n_first]
    second_batch = target_samples[n_first+1:end]

    #####################
    # proposal function #
    #####################

    #Determine proposal function
    post_mean = vec(mean(first_batch))
    post_cov = Array(cov(first_batch)) #TODO: other covariance approximations
    post_cov_pd = PDMat(cholesky(Positive, post_cov))

    proposal_density = MvNormal(post_mean,post_cov_pd)
    proposal_samples = bat_sample(proposal_density,IIDSampling(nsamples=Int(sum(second_batch.weight)))).result
    proposal_density = convert(DistLikeDensity, proposal_density)

    bridge_sampling_integral(target_density,second_batch,proposal_density,proposal_samples,strict,ess_alg)
end


text:
export spglinecurvy

"""
Use: fNew,xNew,rNew,iter,step,err, nProd = spglinecurvy(x,g,fMax,funForward, funPenalty, b,project,tau, params)\n
Projected backtracking linesearch. On entry, g is the (possibly scaled) steepest
descent direction.
"""
function spglinecurvy(A::TA,
                      x::AbstractArray{ETx},
                      g::AbstractArray{ETg},
                      fMax::Tf,
                      funForward::Function,
                      funPenalty::Function,
                      b::AbstractVector{ETb},
                      project::Function,
                      timeProject::Float64,
                      tau::AbstractFloat,
                      options::spgOptions,
                      params::Dict{String,Any}) where
                        {TA<:Union{joAbstractLinearOperator, AbstractArray},
                         Tf<:Number, ETx<:Number, ETg<:Number, ETb<:Number}

    (options.verbosity > 1) && println("Script entered spglinecurvy for A::AbstractArray")

    nProd = 0
    nProdAt = 0
    EXIT_CONVERGED = 0
    EXIT_ITERATIONS = 1
    EXIT_NODESCENT = 2
    gamma = 1e-4
    maxIts = 10
    step = one(real(ETg))
    sNorm = zero(real(ETx))
    scale = one(real(ETg))
    nSafe = 0
    iter = 0
    debug = false
    n = length(x)

    # Define outputs outside of while loop scope
    xNew = similar(x)
    rNew = Array{ETb,1}()
    fNew = zero(Tf)
    err = -1
    while true

        xNew, tmp_itr = project(x - step*scale*g, tau, timeProject, options, params)
        #rNew = b - funForward(A, xNew, Array{ETx,1}(), params)::Array{ETb,1}
        rNew = b - funForward(A, xNew, Array{ETx,1}(), params)
        nProd += 1
        fNew, dummy_g = funPenalty(rNew, params)

        s = xNew - x
        gts = scale * real(g'*s)

        if gts >= 0
            err = EXIT_NODESCENT
            break
        end

        if fNew < fMax + gamma*step*gts
            err = EXIT_CONVERGED
            break
        elseif iter >= maxIts
            err = EXIT_ITERATIONS
            break
        end

        (options.verbosity > 1) && println("Line Search Curvy Iter: $iter")
        #New linesearch iteration
        iter += 1
        step /= 2

        # Safeguard: If stepMax is huge, then even damped search
        # directions can give exactly the same point after projection.  If
        # we observe this in adjacent iterations, we drastically damp the
        # next search direction.

        sNormOld = sNorm
        sNorm = norm(s) / sqrt(n)

        if abs(sNorm - sNormOld) <= 1e-6*sNorm
            gNorm = norm(g) / sqrt(n)
            scale = sNorm/gNorm/(2^nSafe)
            nSafe+=1
        end

    end

    return fNew, xNew, rNew, iter, step, err, nProd
end

function spglinecurvy(A::TA,
                      x::AbstractArray{ETx},
                      g::AbstractArray{ETg},
                      fMax::Tf,
                      funForward::Function,
                      funPenalty::Function,
                      b::AbstractVector{ETb},
                      project::Function,
                      timeProject::Float64,
                      tau::AbstractFloat,
                      options::spgOptions,
                      params::Dict{String,Any}) where
                        {TA<:Function, Tf<:Number,
                         ETx<:Number,ETg<:Number, ETb<:Number}

    (options.verbosity > 1) && println("Script entered spglinecurvy for A::Function")

    nProd = 0
    nProdAt = 0
    EXIT_CONVERGED = 0
    EXIT_ITERATIONS = 1
    EXIT_NODESCENT = 2
    gamma = 1e-4
    maxIts = 10
    step = one(real(ETg))
    sNorm = zero(real(ETx))
    scale = one(real(ETg))
    nSafe = 0
    iter = 0
    debug = false
    n = length(x)

    # Define outputs outside of while loop scope
    xNew = similar(x)
    rNew = Array{ETb,1}()
    fNew = zero(Tf)
    err = -1

    while true

        xNew, tmp_itr = project(x - step*scale*g, tau, timeProject, options, params)#GOOD
        tmp1 = funForward(A, xNew, Array{ETx,2}(undef,0,0), params)[1]
        rNew = b - tmp1   #different from ml
    if false
        return rNew, xNew, b, x, g, step, scale, tau
    end
        nProd += 1

        fNew, dummy_g = funPenalty(rNew, params)
        s = xNew - x
        gts = scale * real(g'*s)

        if gts >= 0
            err = EXIT_NODESCENT
            break
        end
        if fNew < fMax + gamma*step*gts
            err = EXIT_CONVERGED
            break
        elseif iter >= maxIts
            err = EXIT_ITERATIONS
            break
        end

        #New linesearch iteration
        iter += 1
        step /= 2

        # Safeguard: If stepMax is huge, then even damped search
        # directions can give exactly the same point after projection.  If
        # we observe this in adjacent iterations, we drastically damp the
        # next search direction.

        sNormOld = copy(sNorm)
        sNorm = norm(s) / sqrt(n)

        if abs(sNorm - sNormOld) <= 1e-6*sNorm
            gNorm = norm(g) / sqrt(n)
            scale = sNorm/gNorm/(2^nSafe)
            nSafe+=1
        end

    end

    return fNew, xNew, rNew, iter, step, err, nProd
end