text:
function RandomPercentile()
    return RandomPercentile(())
end

function RandomPercentile(arg0::RandomGenerator)
    return RandomPercentile((RandomGenerator,), arg0)
end

function RandomPercentile(arg0::RandomPercentile)
    return RandomPercentile((RandomPercentile,), arg0)
end

function RandomPercentile(arg0::jdouble)
    return RandomPercentile((jdouble,), arg0)
end

function RandomPercentile(arg0::jdouble, arg1::RandomGenerator)
    return RandomPercentile((jdouble, RandomGenerator), arg0, arg1)
end

function aggregate(obj::RandomPercentile, arg0::RandomPercentile)
    return jcall(obj, "aggregate", void, (RandomPercentile,), arg0)
end

function clear(obj::RandomPercentile)
    return jcall(obj, "clear", void, ())
end

function copy(obj::RandomPercentile)
    return jcall(obj, "copy", RandomPercentile, ())
end

function evaluate(obj::RandomPercentile, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function evaluate(obj::RandomPercentile, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "evaluate", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end

function evaluate(obj::RandomPercentile, arg0::jdouble, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (jdouble, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end

function get_aggregate_n(obj::RandomPercentile, arg0::Collection)
    return jcall(obj, "getAggregateN", jdouble, (Collection,), arg0)
end

function get_aggregate_quantile_rank(obj::RandomPercentile, arg0::jdouble, arg1::Collection)
    return jcall(obj, "getAggregateQuantileRank", jdouble, (jdouble, Collection), arg0, arg1)
end

function get_aggregate_rank(obj::RandomPercentile, arg0::jdouble, arg1::Collection)
    return jcall(obj, "getAggregateRank", jdouble, (jdouble, Collection), arg0, arg1)
end

function get_n(obj::RandomPercentile)
    return jcall(obj, "getN", jlong, ())
end

function get_quantile_rank(obj::RandomPercentile, arg0::jdouble)
    return jcall(obj, "getQuantileRank", jdouble, (jdouble,), arg0)
end

function get_rank(obj::RandomPercentile, arg0::jdouble)
    return jcall(obj, "getRank", jdouble, (jdouble,), arg0)
end

function get_result(obj::RandomPercentile)
    return jcall(obj, "getResult", jdouble, ())
end

function get_result(obj::RandomPercentile, arg0::jdouble)
    return jcall(obj, "getResult", jdouble, (jdouble,), arg0)
end

function increment(obj::RandomPercentile, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end

function max_values_retained(::Type{RandomPercentile}, arg0::jdouble)
    return jcall(RandomPercentile, "maxValuesRetained", jlong, (jdouble,), arg0)
end

function reduce(obj::RandomPercentile, arg0::jdouble, arg1::Collection)
    return jcall(obj, "reduce", jdouble, (jdouble, Collection), arg0, arg1)
end


text:
struct NoBackend <: AbstractBackend end

const _backendType = Dict{Symbol, DataType}(:none => NoBackend)
const _backendSymbol = Dict{DataType, Symbol}(NoBackend => :none)
const _backends = Symbol[]
const _initialized_backends = Set{Symbol}()

"Returns a list of supported backends"
backends() = _backends

"Returns the name of the current backend"
backend_name() = CURRENT_BACKEND.sym
_backend_instance(sym::Symbol) = haskey(_backendType, sym) ? _backendType[sym]() : error("Unsupported backend $sym")

macro init_backend(s)
    str = lowercase(string(s))
    sym = Symbol(str)
    T = Symbol(string(s) * "Backend")
    esc(quote
        struct $T <: AbstractBackend end
        export $sym
        $sym(; kw...) = (default(; kw...); backend(Symbol($str)))
        backend_name(::$T) = Symbol($str)
        push!(_backends, Symbol($str))
        _backendType[Symbol($str)] = $T
        _backendSymbol[$T] = Symbol($str)
        include("backends/" * $str * ".jl")
    end)
end

include("backends/web.jl")
# include("backends/supported.jl")

# ---------------------------------------------------------

function add_backend(pkg::Symbol)
    info("To do a standard install of $pkg, copy and run this:\n\n")
    println(add_backend_string(_backend_instance(pkg)))
    println()
end
add_backend_string(b::AbstractBackend) = warn("No custom install defined for $(backend_name(b))")

# don't do anything as a default
_create_backend_figure(plt::Plot) = nothing
_prepare_plot_object(plt::Plot) = nothing
_initialize_subplot(plt::Plot, sp::Subplot) = nothing

_series_added(plt::Plot, series::Series) = nothing
_series_updated(plt::Plot, series::Series) = nothing

_before_layout_calcs(plt::Plot) = nothing

title_padding(sp::Subplot) = sp[:title] == "" ? 0mm : sp[:titlefontsize] * pt
guide_padding(axis::Axis) = axis[:guide] == "" ? 0mm : axis[:guidefontsize] * pt

"Returns the (width,height) of a text label."
function text_size(lablen::Int, sz::Number, rot::Number = 0)
    # we need to compute the size of the ticks generically
    # this means computing the bounding box and then getting the width/height
    # note:
    ptsz = sz * pt
    width = 0.8lablen * ptsz

    # now compute the generalized "height" after rotation as the "opposite+adjacent" of 2 triangles
    height = abs(sind(rot)) * width + abs(cosd(rot)) * ptsz
    width = abs(sind(rot+90)) * width + abs(cosd(rot+90)) * ptsz
    width, height
end
text_size(lab::AbstractString, sz::Number, rot::Number = 0) = text_size(length(lab), sz, rot)

# account for the size/length/rotation of tick labels
function tick_padding(axis::Axis)
    ticks = get_ticks(axis)
    if ticks == nothing
        0mm
    else
        vals, labs = ticks
        isempty(labs) && return 0mm
        # ptsz = axis[:tickfont].pointsize * pt
        longest_label = maximum(length(lab) for lab in labs)

        # generalize by "rotating" y labels
        rot = axis[:rotation] + (axis[:letter] == :y ? 90 : 0)

        # # we need to compute the size of the ticks generically
        # # this means computing the bounding box and then getting the width/height
        # labelwidth = 0.8longest_label * ptsz
        #
        #
        # # now compute the generalized "height" after rotation as the "opposite+adjacent" of 2 triangles
        # hgt = abs(sind(rot)) * labelwidth + abs(cosd(rot)) * ptsz + 1mm
        # hgt

        # get the height of the rotated label
        text_size(longest_label, axis[:tickfontsize], rot)[2]
    end
end

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot)
    # TODO: something different when `is3d(sp) == true`
    leftpad   = tick_padding(sp[:yaxis]) + sp[:left_margin]   + guide_padding(sp[:yaxis])
    toppad    = sp[:top_margin]    + title_padding(sp)
    rightpad  = sp[:right_margin]
    bottompad = tick_padding(sp[:xaxis]) + sp[:bottom_margin] + guide_padding(sp[:xaxis])

    # switch them?
    if sp[:xaxis][:mirror]
        bottompad, toppad = toppad, bottompad
    end
    if sp[:yaxis][:mirror]
        leftpad, rightpad = rightpad, leftpad
    end

    # @show (leftpad, toppad, rightpad, bottompad)
    sp.minpad = (leftpad, toppad, rightpad, bottompad)
end

_update_plot_object(plt::Plot) = nothing

# ---------------------------------------------------------

mutable struct CurrentBackend
  sym::Symbol
  pkg::AbstractBackend
end
CurrentBackend(sym::Symbol) = CurrentBackend(sym, _backend_instance(sym))

# ---------------------------------------------------------

function pickDefaultBackend()
    env_default = get(ENV, "PLOTS_DEFAULT_BACKEND", "")
    if env_default != ""
        try
            Pkg.installed(env_default)  # this will error if not installed
            sym = Symbol(lowercase(env_default))
            if haskey(_backendType, sym)
                return backend(sym)
            else
                warn("You have set PLOTS_DEFAULT_BACKEND=$env_default but it is not a valid backend package.  Choose from:\n\t",
                     join(sort(_backends), "\n\t"))
            end
        catch
            warn("You have set PLOTS_DEFAULT_BACKEND=$env_default but it is not installed.")
        end
    end

    # the ordering/inclusion of this package list is my semi-arbitrary guess at
    # which one someone will want to use if they have the package installed...accounting for
    # features, speed, and robustness
    for pkgstr in ("GR", "PyPlot", "PlotlyJS", "PGFPlots", "UnicodePlots", "InspectDR", "GLVisualize")
        if Pkg.installed(pkgstr) != nothing
            return backend(Symbol(lowercase(pkgstr)))
        end
    end

    # the default if nothing else is installed
    backend(:plotly)
end

# ---------------------------------------------------------

"""
Returns the current plotting package name.  Initializes package on first call.
"""
function backend()

  global CURRENT_BACKEND
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end

  sym = CURRENT_BACKEND.sym
  if !(sym in _initialized_backends)

    # # initialize
    # println("[Plots.jl] Initializing backend: ", sym)

    inst = _backend_instance(sym)
    try
      _initialize_backend(inst)
    catch err
      warn("Couldn't initialize $sym.  (might need to install it?)")
      add_backend(sym)
      rethrow(err)
    end

    push!(_initialized_backends, sym)

  end
  CURRENT_BACKEND.pkg
end

"""
Set the plot backend.
"""
function backend(pkg::AbstractBackend)
    CURRENT_BACKEND.sym = backend_name(pkg)
    warn_on_deprecated_backend(CURRENT_BACKEND.sym)
    CURRENT_BACKEND.pkg = pkg
    backend()
end

function backend(modname::Symbol)
    warn_on_deprecated_backend(modname)
    CURRENT_BACKEND.sym = modname
    CURRENT_BACKEND.pkg = _backend_instance(modname)
    backend()
end

const _deprecated_backends = [:qwt, :winston, :bokeh, :gadfly, :immerse]

function warn_on_deprecated_backend(bsym::Symbol)
    if bsym in _deprecated_backends
        warn("Backend $bsym has been deprecated.  It may not work as originally intended.")
    end
end

# ---------------------------------------------------------

# these are args which every backend supports because they're not used in the backend code
const _base_supported_args = [
    :color_palette,
    :background_color, :background_color_subplot,
    :foreground_color, :foreground_color_subplot,
    :group,
    :seriestype,
    :seriescolor, :seriesalpha,
    :smooth,
    :xerror, :yerror,
    :subplot,
    :x, :y, :z,
    :show, :size,
    :margin,
    :left_margin,
    :right_margin,
    :top_margin,
    :bottom_margin,
    :html_output_format,
    :layout,
    :link,
    :primary,
    :series_annotations,
    :subplot_index,
    :discrete_values,
    :projection,

]

function merge_with_base_supported(v::AVec)
    v = vcat(v, _base_supported_args)
    for vi in v
        if haskey(_axis_defaults, vi)
            for letter in (:x,:y,:z)
                push!(v, Symbol(letter,vi))
            end
        end
    end
    Set(v)
end

# @init_backend Immerse
# @init_backend Gadfly
@init_backend PyPlot
# @init_backend Qwt
@init_backend UnicodePlots
# @init_backend Winston
# @init_backend Bokeh
@init_backend Plotly
@init_backend PlotlyJS
@init_backend GR
@init_backend GLVisualize
@init_backend PGFPlots
@init_backend InspectDR
@init_backend HDF5

# ---------------------------------------------------------

# create the various `is_xxx_supported` and `supported_xxxs` methods
# by default they pass through to checking membership in `_gr_xxx`
for s in (:attr, :seriestype, :marker, :style, :scale)
    f = Symbol("is_", s, "_supported")
    f2 = Symbol("supported_", s, "s")
    @eval begin
        $f(::AbstractBackend, $s) = false
        $f(bend::AbstractBackend, $s::AbstractVector) = all(v -> $f(bend, v), $s)
        $f($s) = $f(backend(), $s)
        $f2() = $f2(backend())
    end

    for bend in backends()
        bend_type = typeof(_backend_instance(bend))
        v = Symbol("_", bend, "_", s)
        @eval begin
            $f(::$bend_type, $s::Symbol) = $s in $v
            $f2(::$bend_type) = $v
        end
    end
end

# is_subplot_supported(::AbstractBackend) = false
# is_subplot_supported() = is_subplot_supported(backend())

ext: jl


text:
using DimensionalData, Tables, Test, DataFrames

using DimensionalData.LookupArrays, DimensionalData.Dimensions
using DimensionalData: DimTable, DimColumn, DimArrayColumn, dimstride

x = X([:a, :b, :c])
y = Y([10.0, 20.0])
d = Dim{:test}(1.0:1.0:3.0)
dimz = x, y, d
da = DimArray(ones(3, 2, 3), dimz; name=:data)
da2 = DimArray(fill(2, (3, 2, 3)), dimz; name=:data2)

@testset "DimArray Tables interface" begin
    ds = DimStack(da)
    t = Tables.columns(ds)
    @test t isa DimTable
    @test dims(t) == dims(da)

    @test Tables.columns(t) === t
    @test t[:X] isa DimColumn
    @test t[:data] isa DimArrayColumn
    @test length(t[:X]) == length(t[:Y]) == length(t[:test]) == 18

    @test Tables.istable(t) == Tables.istable(da) == Tables.istable(ds) == true
    @test Tables.columnaccess(t) == Tables.columnaccess(da) == Tables.columnaccess(ds) == true
    @test Tables.rowaccess(t) == Tables.rowaccess(ds) == Tables.rowaccess(ds) == false
    @test Tables.columnnames(t) == Tables.columnnames(da) == Tables.columnnames(ds) == (:X, :Y, :test, :data)

    sa = Tables.schema(da)
    sds = Tables.schema(ds)
    st = Tables.schema(t)
    @test sa.names == sds.names == st.names == (:X, :Y, :test, :data)
    @test sa.types == sds.types == st.types == (Symbol, Float64, Float64, Float64)

    @test Tables.getcolumn(t, 1) == Tables.getcolumn(t, :X) == Tables.getcolumn(t, X) ==
          Tables.getcolumn(ds, 1) == Tables.getcolumn(ds, :X) == Tables.getcolumn(ds, X) ==
          Tables.getcolumn(da, 1) == Tables.getcolumn(da, :X) == Tables.getcolumn(da, X) ==
          Tables.getcolumn(da, 1)[:] == repeat([:a, :b, :c], 6)
    @test Tables.getcolumn(t, 2) == Tables.getcolumn(t, :Y) ==
          Tables.getcolumn(da, 2) == Tables.getcolumn(da, :Y) ==
          Tables.getcolumn(ds, 2) == Tables.getcolumn(ds, :Y) ==
          Tables.getcolumn(ds, 2)[:] == repeat([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], 3)
    @test Tables.getcolumn(t, 3) == Tables.getcolumn(t, :test) ==
          Tables.getcolumn(da, 3) == Tables.getcolumn(da, :test) ==
          Tables.getcolumn(ds, 3) == Tables.getcolumn(ds, :test) ==
          Tables.getcolumn(ds, 3)[:] == vcat(repeat([1.0], 6), repeat([2.0], 6), repeat([3.0], 6))
    @test Tables.getcolumn(t, 4) == Tables.getcolumn(t, :data) ==
          Tables.getcolumn(da, 4) == Tables.getcolumn(da, :data) ==
          Tables.getcolumn(ds, 4) == Tables.getcolumn(ds, :data) ==
          Tables.getcolumn(ds, 4)[:] == ones(3 * 2 * 3)
    @test Tables.getcolumn(t, Float64, 4, :data) == ones(3 * 2 * 3)
    @test Tables.getcolumn(t, Float64, 2, :Y) == Tables.getcolumn(da, Float64, 2, :Y) ==
          Tables.getcolumn(ds, Float64, 2, :Y) ==
          Tables.getcolumn(ds, Float64, 2, :Y)[:] == repeat([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], 3)
    @test_throws ArgumentError Tables.getcolumn(t, :NotAColumn)
    @test_throws BoundsError Tables.getcolumn(t, 5)
end

@testset "DimColumn" begin
    c = DimColumn(dims(da, Y), dims(da))
    @test length(c) == length(da)
    @test size(c) == (length(da),)
    @test axes(c) == (Base.OneTo(length(da)),)
    @test vec(c) == Array(c) == Vector(c) ==
        repeat([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], 3)
    @test c[1] == 10.0
    @test c[4] == 20.0
    @test c[7] == 10.0
    @test c[18] == 20.0
    @test c[1:5] == [10.0, 10.0, 10.0, 20.0, 20.0]
    @test_throws BoundsError c[-1]
    @test_throws BoundsError c[19]

    cX = DimColumn(dims(da, X), dims(da))
    @test vec(cX) == Array(cX) == Vector(cX) == repeat([:a, :b, :c], 6)
end

@testset "DataFrame conversion" begin
    ds = DimStack(da, da2)
    @time t = DimTable(ds)
    @time df = DataFrame(t; copycols=true)
    @test names(df) == ["X", "Y", "test", "data", "data2"]
    @test Tables.columntype(df, :X) == Symbol
    @test Tables.columntype(df, :data) == Float64
    @test Tables.columntype(df, :data2) == Int

    @test Tables.getcolumn(df, 1)[:] == Tables.getcolumn(df, :X)[1:18] ==
        repeat([:a, :b, :c], 6)
    @test Tables.getcolumn(t, 2) == Tables.getcolumn(df, :Y) ==
        repeat([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], 3)
    @test Tables.getcolumn(df, 3) == Tables.getcolumn(df, :test) ==
        vcat(repeat([1.0], 6), repeat([2.0], 6), repeat([3.0], 6))
    @test Tables.getcolumn(t, 4) == Tables.getcolumn(t, :data) ==
        ones(3 * 2 * 3)
    @test Tables.getcolumn(t, 5) == Tables.getcolumn(t, :data2) ==
        fill(2, 3 * 2 * 3)
end

@testset "Mixed size" begin
    da1 = DimArray(reshape(11:28, (3, 2, 3)), (x, y, d); name=:data1)
    da2 = DimArray(reshape(1.0:6.0, (2, 3)), (y, d); name=:data2)
    ds = DimStack(da1, da2)
    @time t = DimTable(ds)
    @time df = DataFrame(t; copycols=true)
    @test names(df) == ["X", "Y", "test", "data1", "data2"]
    @test Tables.columntype(df, :X) == Symbol
    @test Tables.columntype(df, :data1) == Int
    @test Tables.columntype(df, :data2) == Float64

    @test Tables.getcolumn(df, 1)[:] == Tables.getcolumn(df, :X)[1:18] ==
        repeat([:a, :b, :c], 6)
    @test Tables.getcolumn(t, 2) == Tables.getcolumn(df, :Y) ==
        repeat([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], 3)
    @test Tables.getcolumn(df, 3) == Tables.getcolumn(df, :test) ==
        vcat(repeat([1.0], 6), repeat([2.0], 6), repeat([3.0], 6))
    @test Tables.getcolumn(t, 4) == Tables.getcolumn(t, :data1) == 11:28
    @test Tables.getcolumn(t, 5) == Tables.getcolumn(t, :data2) == vcat(([x, x, x] for x in 1.0:6.0)...)
end

@testset "dim methods" begin
    ds = DimStack(da)
    @test dims(ds) == dims(da)
    @test lookup(ds) == lookup(dims(da))
end


text:
"""
    getdata_enzstepsize([enz_stepsizetotest, seed]) -> DataFrame

Run controlled thick and thin enz experiments at different enz step sizes
as defined by an array input.

Other simulation controls include:
    set_MAX_STEPS_PER_WALK()
    set_NUMBER_OF_WALKS()
Lab data shows thin enzyme has ~ double the yield of thick enzyme, so
target when the sensor receives half as many peroxides at the thicker parameter
"""
function getdata_enzstepsize(enz_stepsizetotest::Array{Float64, 1}=[0.005], seed::Int=randseed())
    enz_print(seed)
    ss_arr, sensor_arr, escape_arr, unresv_arr, thick_arr = ([], [], [], [], [])
    for ss in enz_stepsizetotest
        for thickorthin in [true, false]
            set_THICK_ENZ(thickorthin)
            set_ENZ_STEP_SIZE(ss)
            append!(ss_arr, ss)

            data = Dict{String,Integer}()
            data["sensor"]  = 0
            data["escape"] = 0
            data["particles unresolved"] = 0
            data = run_sim!(data, seed)
            append!(sensor_arr, data["sensor"])
            append!(escape_arr, data["escape"])
            append!(unresv_arr, data["particles unresolved"])
            if thickorthin
                push!(thick_arr, "thick, $ENZYME_MAX_Y um")
            else
                push!(thick_arr, "thin, $ENZYME_MAX_Y um")
            end
        end
    end
    df = DataFrame(
        thick_or_thin = thick_arr,
        enz_step_size = ss_arr,
        sensor_yield = sensor_arr,
        escaped = escape_arr,
        unresolved = unresv_arr)
    current_df(df)
    current_seed(seed)
    current_path("out/enzdata/")
    return df
end

"Helper function, not exported in Enz module."
function enz_print(seed::Int)
    println("############################")
    println("   Compare Enz Step Sizes   ")
    println("############################")
    println("_________Parameters_________")
    println("# of trials:\t\t", string(NUMBER_OF_WALKS))
    println("# of steps (max):\t", string(MAX_STEPS_PER_WALK))
    println("Step size, water:\t", string(WATER_STEP_SIZE))
    println("Random seed:\t\t$seed")
    nothing
end

"""
    enz_sim([seed])

Run one simulation of the enz experiment with one sensor and no walls.

Simulation controls:
    set_THICK_ENZ() -> 2.0 um when true, 0.2 um when false
    set_ENZ_STEP_SIZE()
    set_NUMBER_OF_WALKS()
    set_MAX_STEPS_PER_WALK()
"""
function enz_sim(seed::Int=randseed())
    println("############################")
    println("       Enz Experiment       ")
    println("############################")
    println("_________Parameters_________")
    println("Enzyme Thickness:\t", string(ENZYME_MAX_Y))
    println("# of trials:\t\t", string(NUMBER_OF_WALKS))
    println("# of steps (max):\t", string(MAX_STEPS_PER_WALK))
    println("Step size, water:\t", string(WATER_STEP_SIZE))
    println("Step size, enz:\t\t", string(ENZ_STEP_SIZE))
    println("Random seed:\t\t$seed")

    data = Dict{String,Integer}()
    data["sensor"] = 0
    data["escape"] = 0
    data["particles unresolved"] = 0
    output = run_sim!(data, seed)

    println("_________Results____________")
    println("# in sensor:\t\t", string(data["sensor"]))
    println("# of escaped:\t\t", string(data["escape"]))
    println("# unresolved:\t\t", string(data["particles unresolved"]))
    nothing
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
"""
Add variables to the OptimizationContainer for any component.
"""
function add_variables!(
    optimization_container::OptimizationContainer,
    ::Type{T},
    devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}},
    formulation::Union{AbstractDeviceFormulation, AbstractServiceFormulation},
) where {T <: AuxVariableType, U <: PSY.Component}
    add_variable!(optimization_container, T(), devices, formulation)
end

@doc raw"""
Default implementation of adding auxiliary variable to the model.
"""
function add_variable!(
    optimization_container::OptimizationContainer,
    ::T,
    devices::U,
    formulation,
) where {
    T <: AuxVariableType,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    @assert !isempty(devices)
    time_steps = model_time_steps(optimization_container)
    add_aux_var_container!(
        optimization_container,
        AuxVarKey(T, D),
        [PSY.get_name(d) for d in devices],
        time_steps,
    )
    return
end


text:
using Distributions
using Random
using LinearAlgebra
using PDMats
using Statistics
using Test
import JSON
import Distributions: _univariate, _multivariate, _rand_params

@testset "matrixvariates" begin
#=
    1. baseline tests
    2. compare 1 x 1 matrix-variate with univariate
    3. compare row/column matrix-variate with multivariate
    4. compare density evaluation against archived Stan output
    5. special, distribution-specific tests
    6. main testing method
    7. run matrix-variate unit tests
=#

#  =============================================================================
#  1. baseline test
#  =============================================================================

#  --------------------------------------------------
#  Check that a random draw from d has the right properties
#  --------------------------------------------------

function test_draw(d::MatrixDistribution, X::AbstractMatrix)
    @test size(d) == size(X)
    @test size(d, 1) == size(X, 1)
    @test size(d, 2) == size(X, 2)
    @test length(d) == length(X)
    @test rank(d) == rank(X)
    @test insupport(d, X)
    @test logpdf(d, X) ≈ log(pdf(d, X))
    @test logpdf(d, [X, X]) ≈ log.(pdf(d, [X, X]))
    @test loglikelihood(d, X) ≈ logpdf(d, X)
    @test loglikelihood(d, [X, X]) ≈ 2 * logpdf(d, X)
    if d isa MatrixFDist
        # Broken since `pdadd` is not defined for SubArray
        @test_broken loglikelihood(d, cat(X, X; dims=3)) ≈ 2 * logpdf(d, X)
    else
        @test loglikelihood(d, cat(X, X; dims=3)) ≈ 2 * logpdf(d, X)
    end
    nothing
end

test_draw(d::MatrixDistribution) = test_draw(d, rand(d))

#  --------------------------------------------------
#  Check that sample quantities are close to population quantities
#  --------------------------------------------------

function test_draws(d::MatrixDistribution, draws::AbstractArray)
    @test isapprox(mean(draws), mean(d), atol = 0.1)
    @test isapprox(cov(hcat(vec.(draws)...)'), cov(d) , atol = 0.1)
    nothing
end

function test_draws(d::LKJ, draws::AbstractArray)
    @test isapprox(mean(draws), mean(d), atol = 0.1)
    @test isapprox(var(draws), var(d) , atol = 0.1)
    nothing
end

function test_draws(d::MatrixDistribution, M::Integer)
    rng = MersenneTwister(123)
    @testset "Testing matrix-variates with $key" for (key, func) in
        Dict("rand(...)" => [rand, rand],
             "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])
        test_draws(d, func[2](d, M))
    end
    @testset "Testing matrix-variates with $key" for (key, func) in
        Dict("rand!(..., false)" => [(dist, X) -> rand!(dist, X, false), rand!],
             "rand!(rng, ..., false)" => [(dist, X) -> rand!(rng, dist, X, false), (dist, X) -> rand!(rng, dist, X)])
        m = [Matrix{partype(d)}(undef, size(d)) for _ in Base.OneTo(M)]
        x = func[1](d, m)
        @test x ≡ m
        @test isapprox(mean(x) , mean(d) , atol = 0.1)
        m3 = Array{partype(d), 3}(undef, size(d)..., M)
        x = func[2](d, m3)
        @test x ≡ m3
        @test isapprox(mean(x) , mean(mean(d)) , atol = 0.1)
    end
    @testset "Testing matrix-variates with $key" for (key, func) in
        Dict("rand!(..., true)" => (dist, X) -> rand!(dist, X, true),
             "rand!(rng, true)" => (dist, X) -> rand!(rng, dist, X, true))
        m = Vector{Matrix{partype(d)}}(undef, M)
        x = func(d, m)
        @test x ≡ m
        @test isapprox(mean(x), mean(d) , atol = 0.1)
    end
    repeats = 10
    m = Vector{Matrix{partype(d)}}(undef, repeats)
    rand!(d, m)
    @test isassigned(m, 1)
    m1=m[1]
    rand!(d, m)
    @test m1 ≡ m[1]
    rand!(d, m, true)
    @test m1 ≢ m[1]
    m1 = m[1]
    rand!(d, m, false)
    @test m1 ≡ m[1]
    nothing
end

#  --------------------------------------------------
#  Check that the convert and partype methods work
#  --------------------------------------------------

function test_convert(d::MatrixDistribution)
    distname = getproperty(parentmodule(typeof(d)), nameof(typeof(d)))
    @test distname(params(d)...) == d
    @test d == deepcopy(d)
    for elty in (Float32, Float64, BigFloat)
        del1 = convert(distname{elty}, d)
        del2 = convert(distname{elty}, getfield.(Ref(d), fieldnames(typeof(d)))...)
        @test del1 isa distname{elty}
        @test del2 isa distname{elty}
        @test partype(del1) == elty
        @test partype(del2) == elty
    end
    nothing
end

#  --------------------------------------------------
#  Check that cov and var agree
#  --------------------------------------------------

function test_cov(d::MatrixDistribution)
    @test vec(var(d)) ≈ diag(cov(d))
    @test reshape(cov(d), size(d)..., size(d)...) ≈ cov(d, Val(false))
end

test_cov(d::LKJ) = nothing

#  --------------------------------------------------
#  Check dim
#  --------------------------------------------------

function test_dim(d::MatrixDistribution)
    @test dim(d) == size(d, 1)
    @test dim(d) == size(d, 2)
    @test dim(d) == size(mean(d), 1)
    @test dim(d) == size(mean(d), 2)
end

test_dim(d::Union{MatrixNormal, MatrixTDist}) = nothing

#  --------------------------------------------------
#  RUN EVERYTHING
#  --------------------------------------------------

function test_distr(d::MatrixDistribution, M::Integer)
    test_draw(d)
    test_draws(d, M)
    test_cov(d)
    test_convert(d)
    test_dim(d)
    nothing
end

function test_distr(dist::Type{<:MatrixDistribution},
                    n::Integer,
                    p::Integer,
                    M::Integer)
    d = dist(_rand_params(dist, Float64, n, p)...)
    test_distr(d, M)
    nothing
end

#  =============================================================================
#  2. test matrix-variate against the univariate it collapses to in the 1 x 1 case
#  =============================================================================

function test_against_univariate(D::MatrixDistribution, d::UnivariateDistribution)
    X = rand(D)
    x = X[1]
    @test logpdf(D, X) ≈ logpdf(d, x)
    @test pdf(D, X) ≈ pdf(d, x)
    @test mean(D)[1] ≈ mean(d)
    @test var(D)[1] ≈ var(d)
    nothing
end

function test_draws_against_univariate_cdf(D::MatrixDistribution, d::UnivariateDistribution)
    α = 0.025
    M = 100000
    matvardraws = [rand(D)[1] for m in 1:M]
    @test pvalue_kolmogorovsmirnoff(matvardraws, d) >= α
    nothing
end

test_draws_against_univariate_cdf(D::LKJ, d::UnivariateDistribution) = nothing

function test_against_univariate(dist::Type{<:MatrixDistribution})
    D = dist(_rand_params(dist, Float64, 1, 1)...)
    d = _univariate(D)
    test_against_univariate(D, d)
    test_draws_against_univariate_cdf(D, d)
    nothing
end

#  =============================================================================
#  3. test matrix-variate against the multivariate it collapses to in the row/column case
#  =============================================================================

function test_against_multivariate(D::MatrixDistribution, d::MultivariateDistribution)
    X = rand(D)
    x = vec(X)
    @test logpdf(D, X) ≈ logpdf(d, x)
    @test pdf(D, X) ≈ pdf(d, x)
    @test vec(mean(D)) ≈ mean(d)
    @test cov(D) ≈ cov(d)
    nothing
end

function test_against_multivariate(dist::Union{Type{MatrixNormal}, Type{MatrixTDist}}, n::Integer, p::Integer)
    D = dist(_rand_params(dist, Float64, n, 1)...)
    d = _multivariate(D)
    test_against_multivariate(D, d)
    D = dist(_rand_params(dist, Float64, 1, p)...)
    d = _multivariate(D)
    test_against_multivariate(D, d)
    nothing
end

test_against_multivariate(dist::Type{<:MatrixDistribution}, n::Integer, p::Integer) = nothing

#  =============================================================================
#  4. test density evaluation against archived output from Stan
#  =============================================================================

function jsonparams2dist(dist::Type{Wishart}, dict)
    ν = dict["params"][1][1]
    S = zeros(Float64, dict["dims"]...)
    S[:] = Vector{Float64}( dict["params"][2] )
    return Wishart(ν, S)
end

function jsonparams2dist(dist::Type{InverseWishart}, dict)
    ν = dict["params"][1][1]
    S = zeros(Float64, dict["dims"]...)
    S[:] = Vector{Float64}( dict["params"][2] )
    return InverseWishart(ν, S)
end

function jsonparams2dist(dist::Type{LKJ}, dict)
    d = dict["params"][1][1]
    η = dict["params"][2][1]
    return LKJ(d, η)
end

function unpack_matvar_json_dict(dist::Type{<:MatrixDistribution}, dict)
    d = jsonparams2dist(dist, dict)
    X = zeros(Float64, dict["dims"]...)
    X[:] = Vector{Float64}(dict["X"])
    lpdf = dict["lpdf"][1]
    return d, X, lpdf
end

function test_against_stan(dist::Type{<:MatrixDistribution})
    filename = joinpath(@__DIR__, "ref", "matrixvariates", "jsonfiles", "$(dist)_stan_output.json")
    stan_output = JSON.parsefile(filename)
    K = length(stan_output)
    for k in 1:K
        d, X, lpdf = unpack_matvar_json_dict(dist, stan_output[k])
        @test isapprox(logpdf(d, X), lpdf, atol = 1e-10)
        @test isapprox(logpdf(d, [X, X]), [lpdf, lpdf], atol=1e-8)
        @test isapprox(pdf(d, X), exp(lpdf), atol = 1e-6)
        @test isapprox(pdf(d, [X, X]), exp.([lpdf, lpdf]), atol=1e-6)
    end
    nothing
end

function test_against_stan(dist::Union{Type{MatrixNormal}, Type{MatrixTDist}, Type{MatrixBeta}, Type{MatrixFDist}})
    nothing
end

#  =============================================================================
#  5. special, distribution-specific tests
#  =============================================================================

test_special(dist::Type{<:MatrixDistribution}) = nothing

function test_special(dist::Type{MatrixNormal})
    D = MatrixNormal(_rand_params(MatrixNormal, Float64, 2, 2)...)
    @testset "X ~ MN(M, U, V) ⟺ vec(X) ~ N(vec(M), V ⊗ U)" begin
        d = vec(D)
        test_against_multivariate(D, d)
    end
    @testset "MatrixNormal mode" begin
        @test mode(D) == D.M
    end
    @testset "PDMat mixing and matching" begin
        n = 3
        p = 4
        M = randn(n, p)
        u = rand()
        U_scale = ScalMat(n, u)
        U_dense = Matrix(U_scale)
        U_pd    = PDMat(U_dense)
        U_pdiag = PDiagMat(u*ones(n))
        v = rand(p)
        V_pdiag = PDiagMat(v)
        V_dense = Matrix(V_pdiag)
        V_pd    = PDMat(V_dense)
        UV = kron(V_dense, U_dense)
        baseeval = logpdf(MatrixNormal(M, U_dense, V_dense), M)
        for U in [U_scale, U_dense, U_pd, U_pdiag]
            for V in [V_pdiag, V_dense, V_pd]
                d = MatrixNormal(M, U, V)
                @test cov(d) ≈ UV
                @test logpdf(d, M) ≈ baseeval
            end
        end
    end
    nothing
end

function test_special(dist::Type{Wishart})
    n = 3
    M = 5000
    α = 0.05
    ν, Σ = _rand_params(Wishart, Float64, n, n)
    d = Wishart(ν, Σ)
    H = rand(d, M)
    @testset "deprecations" begin
        for warn in (true, false)
            @test @test_deprecated(Wishart(n - 1, Σ, warn)) == Wishart(n - 1, Σ)
        end
    end
    @testset "meanlogdet" begin
        @test isapprox(Distributions.meanlogdet(d), mean(logdet.(H)), atol = 0.1)
    end
    @testset "H ~ W(ν, Σ), a ~ p(a) independent ⟹ a'Ha / a'Σa ~ χ²(ν)" begin
        q = MvTDist(10, randn(n), rand(d))
        ρ = Chisq(ν)
        A = rand(q, M)
        z = [A[:, m]'*H[m]*A[:, m] / (A[:, m]'*Σ*A[:, m]) for m in 1:M]
        @test pvalue_kolmogorovsmirnoff(z, ρ) >= α
    end
    @testset "H ~ W(ν, I) ⟹ H[i, i] ~ χ²(ν)" begin
        κ = n + 1
        ρ = Chisq(κ)
        g = Wishart(κ, ScalMat(n, 1))
        mymats = zeros(n, n, M)
        for m in 1:M
            mymats[:, :, m] = rand(g)
        end
        for i in 1:n
            @test pvalue_kolmogorovsmirnoff(mymats[i, i, :], ρ) >= α / n
        end
    end
    @testset "Check Singular Branch" begin
        # Check that no warnings are shown: #1410
        rank1 = @test_logs Wishart(n - 2, Σ)
        rank2 = @test_logs Wishart(n - 1, Σ)
        test_draw(rank1)
        test_draw(rank2)
        test_draws(rank1, rand(rank1, 10^6))
        test_draws(rank2, rand(rank2, 10^6))
        test_cov(rank1)
        test_cov(rank2)

        X = H[1]
        @test Distributions.singular_wishart_logkernel(d, X) ≈ Distributions.nonsingular_wishart_logkernel(d, X)
        @test Distributions.singular_wishart_logc0(n, ν, d.S, rank(d)) ≈ Distributions.nonsingular_wishart_logc0(n, ν, d.S)
        @test logpdf(d, X) ≈ Distributions.singular_wishart_logkernel(d, X) + Distributions.singular_wishart_logc0(n, ν, d.S, rank(d))
    end
    nothing
end

function test_special(dist::Type{MatrixTDist})
    @testset "MT(v, M, vΣ, Ω) → MN(M, Σ, Ω) as v → ∞" begin
        n, p = (6, 3)
        M, Σ, Ω = _rand_params(MatrixNormal, Float64, n, p)
        MT = MatrixTDist(1000, M, 1000Σ, Ω)
        MN = MatrixNormal(M, Σ, Ω)
        A = rand(MN)
        @test isapprox(logpdf(MT, A), logpdf(MN, A), atol = 0.1)
    end
    @testset "PDMat mixing and matching" begin
        n = 3
        p = 4
        ν = max(n, p) + 1
        M = randn(n, p)
        u = rand()
        U_scale = ScalMat(n, u)
        U_dense = Matrix(U_scale)
        U_pd    = PDMat(U_dense)
        U_pdiag = PDiagMat(u*ones(n))
        v = rand(p)
        V_pdiag = PDiagMat(v)
        V_dense = Matrix(V_pdiag)
        V_pd    = PDMat(V_dense)
        UV = kron(V_dense, U_dense) ./ (ν - 2)
        baseeval = logpdf(MatrixTDist(ν, M, U_dense, V_dense), M)
        for U in [U_scale, U_dense, U_pd, U_pdiag]
            for V in [V_pdiag, V_dense, V_pd]
                d = MatrixTDist(ν, M, U, V)
                @test cov(d) ≈ UV
                @test logpdf(d, M) ≈ baseeval
            end
        end
    end
    nothing
end

function test_special(dist::Type{LKJ})
    @testset "LKJ mode" begin
        @test mode(LKJ(5, 1.5)) == mean(LKJ(5, 1.5))
        @test_throws ArgumentError mode( LKJ(5, 0.5) )
    end
    @testset "LKJ marginals" begin
        d = 4
        η = abs(3randn())
        G = LKJ(d, η)
        M = 10000
        α = 0.05
        L = sum(1:(d - 1))
        ρ = Distributions._marginal(G)
        mymats = zeros(d, d, M)
        for m in 1:M
            mymats[:, :, m] = rand(G)
        end
        for i in 1:d
            for j in 1:i-1
                @test pvalue_kolmogorovsmirnoff(mymats[i, j, :], ρ) >= α / L
            end
        end
    end
    @testset "LKJ integrating constant" begin
        #  =============
        #  odd non-uniform
        #  =============
        d = 5
        η = 2.3
        lkj = LKJ(d, η)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
        @test lkj.logc0 == -Distributions.lkj_onion_loginvconst(d, η)
        #  =============
        #  odd uniform
        #  =============
        d = 5
        η = 1.0
        lkj = LKJ(d, η)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst_uniform_odd(d)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_vine_loginvconst_uniform(d)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.corr_logvolume(d)
        @test lkj.logc0 == -Distributions.lkj_onion_loginvconst_uniform_odd(d)
        #  =============
        #  even non-uniform
        #  =============
        d = 6
        η = 2.3
        lkj = LKJ(d, η)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
        @test lkj.logc0 == -Distributions.lkj_onion_loginvconst(d, η)
        #  =============
        #  even uniform
        #  =============
        d = 6
        η = 1.0
        lkj = LKJ(d, η)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst_uniform_even(d)
        @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_vine_loginvconst_uniform(d)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
        @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.corr_logvolume(d)
        @test lkj.logc0 == -Distributions.lkj_onion_loginvconst_uniform_even(d)
    end
    @testset "check integrating constant as a volume" begin
        #  d = 2: Lebesgue measure of the set of correlation matrices is 2.
        volume2D = 2
        @test volume2D ≈ exp( Distributions.lkj_onion_loginvconst(2, 1) )
        @test 1 / volume2D ≈ exp( LKJ(2, 1).logc0 )
        #  d = 3: Lebesgue measure of the set of correlation matrices is π²/2.
        #  See here: https://www.jstor.org/stable/2684832
        volume3D = 0.5π^2
        @test volume3D ≈ exp( Distributions.lkj_onion_loginvconst(3, 1) )
        @test 1 / volume3D ≈ exp( LKJ(3, 1).logc0 )
        #  d = 4: Lebesgue measure of the set of correlation matrices is (32/27)π².
        #  See here: https://doi.org/10.4169/amer.math.monthly.123.9.909
        volume4D = (32 / 27)*π^2
        @test volume4D ≈ exp( Distributions.lkj_onion_loginvconst(4, 1) )
        @test 1 / volume4D ≈ exp( LKJ(4, 1).logc0 )
    end
    nothing
end

#  =============================================================================
#  6. main method
#  =============================================================================

function test_matrixvariate(dist::Type{<:MatrixDistribution},
                            n::Integer,
                            p::Integer,
                            M::Integer)
    test_distr(dist, n, p, M)
    test_against_univariate(dist)
    test_against_multivariate(dist, n, p)
    test_against_stan(dist)
    test_special(dist)
    nothing
end

#  =============================================================================
#  7. run unit tests for matrix-variate distributions
#  =============================================================================

matrixvariates = [(MatrixNormal, 2, 4, 10^5),
                  (Wishart, 2, 2, 10^6),
                  (InverseWishart, 2, 2, 10^6),
                  (MatrixTDist, 2, 4, 10^5),
                  (MatrixBeta, 3, 3, 10^5),
                  (MatrixFDist, 3, 3, 10^5),
                  (LKJ, 3, 3, 10^5)]

for distribution in matrixvariates
    dist, n, p, M = distribution
    println("    testing $(dist)")
    @testset "$(dist)" begin
        test_matrixvariate(dist, n, p, M)
    end
end
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ da4b6f2a-923e-11eb-3f34-cd40c2f87cd5
begin
	using Distributions
	using PlutoUI
	using Statistics
	using StatsPlots
end;

# ╔═╡ c032f824-923e-11eb-2e29-29847cbf1063
md"# A Central Limit Theorem"

# ╔═╡ 0f1865b4-93e3-11eb-38a1-a11e4048d827
TableOfContents()

# ╔═╡ 1c8c2488-93e3-11eb-1864-a7f2a7454e1b
md"## Words"

# ╔═╡ 1a896920-9307-11eb-3a80-63564275f568
md"The Central Limit Theorem discussed here describes the distribution of the random variable, the sample mean.  We don't usually think about a sample mean as a random variable, but that's the first major hurdle to understanding any Central Limit Theorem, functions of random variables are random variables.  Any statistic can be thought of as a random variable."

# ╔═╡ 80a3f330-9308-11eb-0ca5-77b0627d968f
md"Since statistics are random variables, they often follow some distribution.  Most times, we don't know the distribution associated with a statistic.  When the statistic is the mean, then this Central Limit Theorem tells us that the distribution of the sample mean is approximately Normal."

# ╔═╡ b64df9c2-9308-11eb-27f5-3d230dc92312
md"According to this Central Limit Theorem, **the distribution of the sample mean is approximately Normal**.  Therefore the density function of the sample mean describes which values any random mean is most likely to take on."

# ╔═╡ debdea5e-9308-11eb-29f7-71208f83e317
md"Note that we will have to massage the math behind this Central Limit Theorem to better visualize what is going on."

# ╔═╡ 4e184eaa-9307-11eb-1af5-4b2fdc6bbc6b
md"We usually think about one dataset being randomly sampled from a distribution.  Here, we'll consider a sample mean as coming from a distribution.  To help us imagine this abstraction, consider you and $$R - 1$$ of your friends each going out to collect your own random sample from the same process.  By the time each of you collects your own sample of size $$N$$ and then calculates your own mean (from your own data), there will be $$R$$ means."

# ╔═╡ 3e9630c6-9309-11eb-3cf5-73b109c4b412
md"This Central Limit Theorem says that the shape of these $$R$$ means will be approximately Normal, and more over that the width of this Normal distribution depends on $$N$$ (and the variance of the theoretical distribution from which you and your friends sampled from) and the approximation will improve as $$R$$ goes off to infinity."

# ╔═╡ 75dab232-9309-11eb-0f2a-ed982699b343
md"This last bit is where the math and our vizualization below differ.  Mathematically, the approximation and the width of the distribution of the sample mean both only depend on the sample size $$N$$.  The mathematical theory says that as $$N \to \infty$$, then the Normal approximation improves and the width of the distribution will decrease with the square root of the sample size $$N$$."

# ╔═╡ dc33209e-9309-11eb-099b-0f6bc414032a
md"To see visualize this process, pick a distribution below that is supposed to represent the distribution you and your $$R-1$$ friends are sampling from.  Then adjust the sample size and the number of friends $$R$$ you have to see the effect of the sampling distribution of the sample mean.  **Notice that this Central Limit Theorem is indifferent to the shape of the theoretical distribution; whatever shape it has, the sampling distribution is still approximately Normal.**"

# ╔═╡ 36c19130-93e3-11eb-0e54-f7c547f3b438
md"## Pictures"

# ╔═╡ 22658db8-923f-11eb-21ed-0bd2f08bf338
md"""Pick a distribution
$(@bind distribution Select([
	"Gamma" => "Gamma",
	"Normal" => "Normal",
	"Uniform" => "Uniform"]))
"""

# ╔═╡ f2d7090a-9301-11eb-06c8-0fa5fc97ce84
md"""
 $$N$$ = $(@bind N Slider(1:301, default = 27, show_value = true)) $$\quad$$
 $$R$$ = $(@bind R Slider(1:1001, default = 101, show_value = true))
"""

# ╔═╡ f69147d2-930a-11eb-3757-07178c0dadc8
md"## Math"

# ╔═╡ 7b0e370a-930a-11eb-2e33-0bed472a2026
html"

"

# ╔═╡ 1f93e9d0-930b-11eb-2220-332796c4523e
md"Assume you and your $$R-1$$ friends are randomly sampling from distribution $$F$$ with mean $$\mu$$ and variance $$\sigma^2$$.  In a simple version of the Central Limit Theorem of sample means, we require $$\sigma^2 < \infty$$ and the original data to be independent and identically distributed, $$X_1, \ldots, X_N \sim_{iid} F$$."

# ╔═╡ ff9123be-93e3-11eb-3fc5-076772232b87
md"Denote the sample mean as $$\hat{\mu} = N^{-1} \sum_{n = 1}^N X_N$$.  Then by this Central Limit Theorem we have that

$\frac{\hat{\mu} - \mu}{\sigma / \sqrt{N}} \sim N(0, 1)$

as $$N \rightarrow \infty$$."

# ╔═╡ bb1c25f4-9301-11eb-26ce-158f5bb94eae
D =
if distribution == "Gamma"
    Gamma(1, 10)
elseif distribution == "Normal"
	Normal(0, 10)
elseif distribution == "Uniform"
	Uniform(1, 6)
end;

# ╔═╡ aa79d4ac-923f-11eb-111d-3dceaf63f63d
begin
	p = plot(D, label = "", linewidth = 2, title = "Theoretical distribution: $(distribution)")
	plot!(p, [mean(D)], [0], t = :scatter, color = :black, label = "mean", markersize = 15, markershape = :utriangle)
end

# ╔═╡ a44bae34-9302-11eb-3367-cb12648d16be
samplemeans = [mean(rand(D, N)) for _ in 1:R]

# ╔═╡ 1966e0fe-9302-11eb-3ec2-0515c69711a1
begin
	h = histogram(samplemeans, legend = false,
		title = "Sampling distribution \n R = $(R) sample means \n each with N = $(N) observations")
	plot!(h, [mean(D)], [0], t = :scatter, color = :black, label = "mean", markersize = 15, markershape = :utriangle)
end

# ╔═╡ 86fde510-9300-11eb-2dd9-c1dc68850b33
onedataset = rand(D, N)

# ╔═╡ Cell order:
# ╟─c032f824-923e-11eb-2e29-29847cbf1063
# ╟─0f1865b4-93e3-11eb-38a1-a11e4048d827
# ╟─1c8c2488-93e3-11eb-1864-a7f2a7454e1b
# ╟─1a896920-9307-11eb-3a80-63564275f568
# ╟─80a3f330-9308-11eb-0ca5-77b0627d968f
# ╟─b64df9c2-9308-11eb-27f5-3d230dc92312
# ╟─debdea5e-9308-11eb-29f7-71208f83e317
# ╟─4e184eaa-9307-11eb-1af5-4b2fdc6bbc6b
# ╟─3e9630c6-9309-11eb-3cf5-73b109c4b412
# ╟─75dab232-9309-11eb-0f2a-ed982699b343
# ╟─dc33209e-9309-11eb-099b-0f6bc414032a
# ╟─36c19130-93e3-11eb-0e54-f7c547f3b438
# ╟─22658db8-923f-11eb-21ed-0bd2f08bf338
# ╟─aa79d4ac-923f-11eb-111d-3dceaf63f63d
# ╟─f2d7090a-9301-11eb-06c8-0fa5fc97ce84
# ╟─1966e0fe-9302-11eb-3ec2-0515c69711a1
# ╟─a44bae34-9302-11eb-3367-cb12648d16be
# ╟─86fde510-9300-11eb-2dd9-c1dc68850b33
# ╟─f69147d2-930a-11eb-3757-07178c0dadc8
# ╟─7b0e370a-930a-11eb-2e33-0bed472a2026
# ╟─1f93e9d0-930b-11eb-2220-332796c4523e
# ╟─ff9123be-93e3-11eb-3fc5-076772232b87
# ╟─bb1c25f4-9301-11eb-26ce-158f5bb94eae
# ╟─da4b6f2a-923e-11eb-3f34-cd40c2f87cd5

max_stars_repo_name: tim-au/pluto

max_issues_repo_name: tim-au/pluto

max_forks_repo_name: tim-au/pluto


text:
using Test
using DelimitedFiles
using SpectralEnvelope
using Random

cd(@__DIR__)
test = readdlm("DNA_data.txt",',')[1,:]
x,y,e = spectral_envelope(test; m=0)
@test round(spectral_envelope(test)[2][5]; digits = 3) == round(0.003; digits = 3)
@test round(get_mappings(test, 0.33)["A"]; digits = 2) == round(0.54; digits = 2)


text:
function gen_current(props,myprop::Prop)

    if issubtype(typeof(props),Tuple)
        confields=Array(Expr,length(props))
        curfields=Array(Expr,length(props))
        initfields=Array(Expr,length(props))
        mylengths=0

        for i=1:length(props)
            confields[i]=speed_con(props[i],mylengths)
            curfields[i]=speed_cur(props[i],mylengths)
            initfields[i]=speed_init(props[i],mylengths)
            mylengths=length(fieldnames(props[i]))-1
        end
    else
        confields=[speed_con(props,0)]
        curfields=[speed_cur(props,0)]
        initfields=[speed_init(props,0)]
    end

    @eval begin

        function con!(p::$(typeof(myprop)),v::Array{Float64,1},internal::Array{Int64,1},dt=.025)

            @fastmath for i=1:length(internal)
                k=internal[i]
                $([confields[j] for j=1:length(confields)]...)
            end

            nothing
        end

        function cur!(p::$(typeof(myprop)),v::Array{Float64,1},im::Array{Float64,1},internal::Array{Int64,1})

            @fastmath for i=1:length(internal)
                k=internal[i]
                $([curfields[j] for j=1:length(curfields)]...)
            end

            nothing
        end

        function init!(p::$(typeof(myprop)),v::Array{Float64,1},internal::Array{Int64,1})
            @fastmath for i=1:length(internal)
                k=internal[i]
                $([initfields[j] for j=1:length(initfields)]...)
            end
        end
    end
end

const ena=50.0
const ek=-77.0

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
using ArgParse

s=ArgParseSettings()
@add_arg_table s begin
    "--part2"
        help="Do part 2"
        action=:store_true
end
parsed_args = parse_args(ARGS, s)

function part1(numbers)
    lastSeen = typemax(Int64)
    increases = 0
    for value in numbers
        if value > lastSeen
            increases +=1
        end
        lastSeen = value
    end
    increases
end

function part2(numbers)
    increases = 0
    for idx in 1:(Base.length(numbers)-3)
        if sum(numbers[idx:idx+2]) < sum(numbers[idx+1:idx+3])
            increases += 1
        end
    end
    increases
end

numbers = [parse(Int64, line) for line in eachline("input/day_1.txt")]
if parsed_args["part2"]
    answer = part2(numbers)
else
    answer = part1(numbers)
end
println("The answer is: $answer")


text:
using jWinReg
using Base.Test

#set test values
parentkey = "HKCU\\Julia"
subkey1 = "test1"
subkey2 = "test2"
dwName = "dwTest"
dwValue = convert(Int32, 2205)
strName = "strTest"
strValue = "Hello!"
mstrName = "mstrTest"
mstrValue = ["Some", "Vectorized", "Words"]
defaValue = "def"

info("regcreatekey - test #1: recursive creation of keys")
res = regcreatekey(parentkey * "\\" * subkey1)
@test res == 1
res = regcreatekey(parentkey * "\\" * subkey2)
@test res == 1
res = regcreatekey(parentkey * "\\" * subkey1)
@test res == 0

info("regwrite - test #2.1: write values")
res = regwrite(parentkey, dwName, dwValue)
@test res
res = regwrite(parentkey, strName, strValue)
@test res
res = regwrite(parentkey, mstrName, mstrValue)
@test res

info("regwrite - test #2.2: write default value")
res = regwrite(parentkey * "\\" * subkey1, "", defaValue)
@test res

info("... dumping values on file ...")
path = joinpath(Pkg.dir("WinReg"), "test", "regdump.reg")
cmd = `cmd /c regedit /e ""$path"" HKEY_CURRENT_USER\\Julia`
#run(cmd)

info("regread - test #3.1: read existent values")
res = regread(parentkey, dwName)
@test res == dwValue
res = regread(parentkey, strName)
@test res == strValue
res = regread(parentkey, mstrName)
@test length(res) == length(mstrValue)
for w in res
	@test (w in mstrValue)
end

info("regread - test #3.2: read default value")
res = regread(parentkey * "\\" * subkey1, "")
@test res == defaValue

info("regread - test #3.3: read not-existent value")
res = regread(parentkey, "noval")
@test res == nothing

info("regkeylist - test #4.1: list existent subkeys of existent key")
res = regkeylist(parentkey)
@test subkey1 in res
@test subkey2 in res

info("regkeylist - test #4.2: list not-existent subkeys of existent key")
res = regkeylist(parentkey * "\\" * subkey2)
@test length(res) == 0

info("regkeylist - test #4.3: list subkeys of not-existent key")
res = regkeylist(parentkey * "\\nokey")
@test res == nothing

info("regvaluelist - test #5.1: list existent values of existent key")
res = regvaluelist(parentkey)
@test length(res) == 3
@test (dwName, UInt32) in res
@test (strName, AbstractString) in res
@test (mstrName, Vector{AbstractString}) in res

info("regvaluelist - test #5.2: list default value of existent key")
res = regvaluelist(parentkey * "\\" * subkey1)
@test ("", AbstractString) in res

info("regvaluelist - test #5.3: list not-existent values of existent key")
res = regvaluelist(parentkey * "\\" * subkey2)
@test length(res) == 0

info("regvaluelist - test #5.4: list values of not-existent key")
res = regvaluelist(parentkey * "\\nokey")
@test res == nothing

info("regdelete - test #6.1: delete key")
res = regdelete(parentkey * "\\" * subkey2)
@test res

info("regdelete - test #6.2: delete not-existent value")
res = regdelete(parentkey * "\\" * subkey1, dwName)
@test !res

info("regdelete - test #6.3: delete value")
res = regdelete(parentkey, dwName)
@test res

info("regdelete - test #6.4: delete default value")
res = regdelete(parentkey * "\\" * subkey1, "")
@test res

info("regdelete - test #6.5: delete keys recursively (error)")
res = regdelete(parentkey)
@test !res

info("regdelete - test #6.6: delete test subkey")
res = regdelete(parentkey * "\\" * subkey1)
@test res

info("regdelete - test #6.7: delete test key with subvalues")
res = regdelete(parentkey)
@test res


text:
using Documenter, ImmersedLayers

ENV["GKSwstype"] = "nul" # removes GKS warnings during plotting

makedocs(
    sitename = "ImmersedLayers.jl",
    doctest = true,
    clean = true,
    modules = [ImmersedLayers],
    pages = [
        "Home" => "index.md",
        "Basics" => ["manual/caches.md",
                     "manual/surfaceops.md",
                     "manual/multbodies.md",
                     "manual/gridops.md",
                     "manual/helmholtz.md",
                     "manual/matrices.md",
                     "manual/dirichlet.md",
                     "manual/utilities.md"
                     ],
         "Setting up PDEs" => ["manual/problems.md",
                                "manual/neumann.md",
                                "manual/stokes.md",
                                "manual/heatconduction.md",
                                "manual/heatconduction-unbounded.md"
                                ],
        #"Internals" => [ "internals/properties.md"]
    ],
    #format = Documenter.HTML(assets = ["assets/custom.css"])
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict()
            )
        ))
    ),
    #assets = ["assets/custom.css"],
    #strict = true
)

#if "DOCUMENTER_KEY" in keys(ENV)
deploydocs(
     repo = "github.com/JuliaIBPM/ImmersedLayers.jl.git",
     target = "build",
     deps = nothing,
     make = nothing
     #versions = "v^"
)
#end


text:
include("../src/PBCM.jl")

cell_general = cell()
cell_design_op = cylindrical_cell_designer(cell_general)
cost = cost_default()
cell_general, cost = convert_all(cell_general, cost, mult)

using DiffEqSensitivity
using PyPlot

per = 30

no_cells_yr = 10
no_cells_yr_intrvl = [no_cells_yr*(1 - (per/100)) , no_cells_yr*(1 + (per/100))]

cost_pos_AM = 18.0
cost_pos_AM_intrvl = [cost_pos_AM*(1 - (per/100)) , cost_pos_AM*(1 + (per/100))]

cost_neg_AM = 12.5
cost_neg_AM_intrvl = [cost_neg_AM*(1 - (per/100)) , cost_neg_AM*(1 + (per/100))]

cost_sep = 1.1
cost_sep_intrvl    = [cost_sep  * (1 - (per/100)) , cost_sep * (1  + (per/100))]

function gsa_cost(arr)

    cost.general_costs.no_units_mfg   = converter(arr[1], mult.units_mfg)
    cost.cell_costs.cathode.AM[1]     = arr[2]
    cost.cell_costs.anode.AM[1]       = arr[3]
    cost.cell_costs.seperator_cost[1] = arr[4]

    dollars_kWh = cost_calc(cell_general, cost, system="Cell", cost_verbosity=0)[1]

    return dollars_kWh
end

# gsa_cost([0.25, 150.0, 130.0, 25.0, 19.0])

gsa_result = gsa(gsa_cost, Sobol(order = [2]), [no_cells_yr_intrvl, cost_pos_AM_intrvl, cost_neg_AM_intrvl, cost_sep_intrvl], N = 10000)

x = [0, 1, 2, 3, 1, 2, 3, 2, 3, 3]
y = [0, 0, 0, 0, 1, 1, 1, 2, 2, 3]

print("\n")
S1 = gsa_result.S1
S2 = gsa_result.S2
print(S1)
print(S2)

z = [S1[1],S2[1],S2[2],S2[3],S1[2],S2[4],S2[5],S1[3],S2[6],S1[4]]
clf()
r = bar3D(x, y, z*0, 0.4, 0.4, z)
xticks([0,1,2,3], ["No Cells","Cost Pos AM", "Cost Neg AM", "Cost Sep"])
yticks([1,2,3,4], ["No Cells","Cost Pos AM", "Cost Neg AM", "Cost Sep"])

title(string("GSA: 2nd Order: ", per, "%"))
zlabel("Sobol Indices")
figure(2)

# savefig("GSA.png", dpi=500)

max_stars_repo_name: BattModels/BatteryCost

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdline

"""
    ''_antdline(;kwargs...)

An AntdLine component.

Keyword arguments:
- `id` (String; optional)
- `annotations` (optional)
- `appendPadding` (Real | Array of Reals; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `color` (optional): . color has the following type: String | Array of Strings | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `connectNulls` (Bool; optional)
- `data` (Array of Dicts; required)
- `height` (Real; optional)
- `isStack` (Bool; optional)
- `label` (optional)
- `legend` (optional)
- `lineStyle` (optional): . lineStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-CN', 'en-US'; optional)
- `meta` (optional)
- `padding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `point` (optional): . point has the following type: lists containing elements 'color', 'shape', 'style'.
Those elements have the following types:
  - `color` (optional): . color has the following type: String | Array of Strings | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
  - `shape` (optional): . shape has the following type: String | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
  - `style` (optional): . style has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `recentlyPointClickRecord` (optional): . recentlyPointClickRecord has the following type: lists containing elements 'timestamp', 'data'.
Those elements have the following types:
  - `timestamp` (Real; optional)
  - `data` (Dict; optional)
- `recentlyTooltipChangeRecord` (optional): . recentlyTooltipChangeRecord has the following type: lists containing elements 'timestamp', 'data'.
Those elements have the following types:
  - `timestamp` (Real; optional)
  - `data` (Array of Dicts; optional)
- `renderer` (a value equal to: 'canvas', 'svg'; optional)
- `seriesField` (String; optional)
- `slider` (optional)
- `smooth` (Bool; optional)
- `stepType` (String; optional)
- `style` (Dict; optional)
- `tooltip` (optional)
- `width` (Real; optional)
- `xAxis` (optional)
- `xField` (String; required)
- `yAxis` (optional)
- `yField` (String; required)
"""
function ''_antdline(; kwargs...)
        available_props = Symbol[:id, :annotations, :appendPadding, :autoFit, :className, :color, :connectNulls, :data, :height, :isStack, :label, :legend, :lineStyle, :loading_state, :locale, :meta, :padding, :point, :recentlyPointClickRecord, :recentlyTooltipChangeRecord, :renderer, :seriesField, :slider, :smooth, :stepType, :style, :tooltip, :width, :xAxis, :xField, :yAxis, :yField]
        wild_props = Symbol[]
        return Component("''_antdline", "AntdLine", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
module TestIssorted

using ParallelRadixSort
using Test

testdata_integers = [
    ("$T[...]", () -> rand(T, 2^20)) for
    T in [UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64]
]

testdata_floats = [("$T[...]", () -> rand(T, 2^20)) for T in [Float32, Float64]]

raw_testdata = """
rand('a':'z', 100)
[rand('a':'z', 10) for _ in 1:100]
[join(rand('a':'z', 10)) for _ in 1:100]
[rand(1:10, 10) for _ in 1:100]
"""

# A vector of `(datalabel, function)`
testdata_misc = map(split(raw_testdata, "\n", keepempty = false)) do code
    fn = "() -> $code"
    @debug "Evaluating: $fn"
    (code, Base.include_string(@__MODULE__, fn))
end

testdata = vcat(testdata_integers, testdata_floats, testdata_misc)

@testset "$alglabel" for (alglabel, alg) in [
    (:MSDRadixSort, ParallelRadixSort.MSDRadixSort),
    (:StableMSDRadixSort, ParallelRadixSort.StableMSDRadixSort),
    (:ParallelMSDRadixSort, ParallelRadixSort.ParallelMSDRadixSort),
    (:ParallelStableMSDRadixSort, ParallelRadixSort.ParallelStableMSDRadixSort),
    (
        "ParallelMSDRadixSort(basesize = 2^10)",
        ParallelRadixSort.ParallelMSDRadixSort(basesize = 2^10),
    ),
]
    @testset "$datalabel" for (datalabel, fn) in testdata
        vector = fn()
        @test issorted(sort(vector; alg = alg))
        eltype(vector) isa Number || continue
        @test issorted(sort(vector; alg = alg, rev = true); rev = true)
        @test issorted(sort!(copy(vector); alg = alg, rev = true); rev = true)
    end
end

end  # module


text:
"""
Things this needs to abstract:

  abstract type LinearSolver

Preconditioning:
  abstract type AbstractPC
  abstract type AbstractPetscMatPC
  abstract type AbstractPetscMatFreePC
  union of those two

Linear operators:
  abstract type AbstractLO
  abstract type AbstractDenseLO
  abstract type AbstractSparseDirectLO
  abstract type AbstractPetscMatLO
  abstract type AbstractPetscMatFreeLO
  union of those two

functions:
  calcLinearOperator  - calculate Jacobian, or if matfree, calc some stuff
  applyLinearOperator  - does Ax=b
  applyLinearOperatorTranspose  - does A'x = b
  needParallelData
  free

Todo:
  which implements the others (type -> function)
  signatures for functions

"""

max_stars_repo_name: OptimalDesignLab/CSEInterfaces

max_issues_repo_name: OptimalDesignLab/CSEInterfaces

max_forks_repo_name: OptimalDesignLab/CSEInterfaces


text:
recipe = Recipe{Assembler}("piercing-rounds-magazine", 3)
@RIN "firearm-magazine" 1
@RIN "steel-plate" 1
@RIN "copper-plate" 5

recipe = Recipe{Assembler}("uranium-rounds-magazine", 10)
@RIN "piercing-rounds-magazine" 1
@RIN "uranium-238" 1

recipe = Recipe{Assembler}("rocket", 8)
@RIN "electronic-circuit" 1
@RIN "explosives" 1
@RIN "iron-plate" 2

recipe = Recipe{Assembler}("explosive-rocket", 8)
@RIN "rocket" 1
@RIN "explosives" 2

recipe = Recipe{Assembler}("atomic-bomb", 50)
@RIN "processing-unit" 20
@RIN "explosives" 10
@RIN "uranium-235" 30

recipe = Recipe{Assembler}("shotgun-shell", 3)
@RIN "copper-plate" 2
@RIN "iron-plate" 2

recipe = Recipe{Assembler}("piercing-shotgun-shell", 8)
@RIN "shotgun-shell" 2
@RIN "copper-plate" 5
@RIN "steel-plate" 2

recipe = Recipe{Assembler}("railgun-dart", 8)
@RIN "steel-plate" 5
@RIN "electronic-circuit" 5

recipe = Recipe{Assembler}("cannon-shell", 0.5)

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
import Simplex: nrows, ncols

@testset "Cplex.jl" begin
    """
    Original LP:

    min  -10 x1 - 12 x2 - 12 x3
    s.t. -20 <= - x1 - 2 x2 - 2 x3
            2 x1 +   x2 + 2 x3 <= 20
        10 <= 2 x1 + 2 x2 +   x3 <= 20
        -10 <= - x1 -   x2
        x1, x2 >= 0
        -1 <= x3 <= 10
    """
    A = sparse([1; 2; 3; 4; 1; 2; 3; 4; 1; 2; 3], [1; 1; 1; 1; 2; 2; 2; 2; 3; 3; 3], [-1.; 2.; 2.; -1; -2; 1.; 2.; -1.; -2.; 2.; 1.])
    bl = [-20.; -Inf; 10.; -10.]
    bu = [Inf; 20.; 20.; Inf]
    c = [-10.; -12.; -12.]
    xl = [0.; 0.; -1.]; xu = [Inf; Inf; 10.]
    @testset "cplex_basis" begin
        lp = MatOI.LPForm(MOI.MIN_SENSE, c, A, bl, bu, xl, xu)
        Simplex.presolve(lp)
        Simplex.scaling!(lp)
        B = Simplex.PhaseOne.Cplex.cplex_basis(lp)
        @test B == [4,5,6,7]

        """
        canonical:
            minimize -10 x1 - 12 x2 - 12 x3
            subject to
                -0.5 x1     - x2     - x3 - x4 == -10
                     x1 + 0.5 x2     + x3 + x5 == 10
                     x1     + x2 + 0.5 x3 - x6 == 5
                    -x1     - x2          - x7 == -10
                     x1, x2, x4, x5, x7 >= 0
                     -1 <= x3 <= 10
                     0 <= x6 <= 5
        """
        canonical = Simplex.canonical_form(lp)
        @test canonical.A == sparse([
            -0.5 -1.0 -1.0 -1.0 0.0 0.0 0.0;
            1.0 0.5 1.0 0.0 1.0 0.0 0.0;
            1.0 1.0 0.5 0.0 0.0 -1.0 0.0;
            -1.0 -1.0 0.0 0.0 0.0 0.0 -1.0])
        @test canonical.v_lb ≈ [0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0]
        @test canonical.v_ub ≈ [Inf, Inf, 10.0, Inf, Inf, 5.0, Inf]
        @test canonical.b ≈ [-10.0, 10.0, 5.0, -10.0]

        """
        phase-one:
            minimize 0
            subject to
                c1:-0.5 x1     - x2     - x3 - x4                == -10
                c2:     x1 + 0.5 x2     + x3      + x4           == 10
                c3:     x1     + x2 + 0.5 x3           - x6      == 5
                c4:   - x1     - x2                         - x7 == -10
                x1, x2, x4, x5, x7 >= 0
                -1 <= x3 <= 10
                0 <= x6 <= 5

        """
        p1lp, artif_col_idx = Simplex.PhaseOne.Artificial.reformulate(canonical, basis = B)
        @test p1lp.A == canonical.A
        @test p1lp.v_lb == canonical.v_lb
        @test p1lp.v_ub == canonical.v_ub
        @test p1lp.b == canonical.b
        @test p1lp.c == zeros(ncols(p1lp))

        x, is_feasible = Simplex.PhaseOne.Cplex.compute_x(p1lp, B)
        @test x ≈ Float64[0, 0, -1, 11, 11, -5.5, 10]
        @test is_feasible == false

        num_xvars = ncols(canonical)
        num_artif = ncols(p1lp) - num_xvars
        @test num_xvars == 7
        @test num_artif == 0

        """
        CPLEX basis LP:
            minimize x8
            subject to
                c1:-0.5 x1     - x2     - x3 - x4                == -10
                c2:     x1 + 0.5 x2     + x3      + x5           == 10
                c3:     x1     + x2 + 0.5 x3           - x6      == 5
                c4:   - x1     - x2                         - x7 == -10
                c5:                                      x6 + x8 - x9 == 0
                        x8..x9 >= 0
        """
        cpxlp, newB, artif_col_idx = Simplex.PhaseOne.Cplex.reformulate(p1lp, x, num_xvars, num_artif, B)
        @test cpxlp.A == sparse(
            [-0.5 -1.0 -1.0 -1.0 0.0  0.0  0.0 0.0  0.0;
              1.0  0.5  1.0  0.0 1.0  0.0  0.0 0.0  0.0;
              1.0  1.0  0.5  0.0 0.0 -1.0  0.0 0.0  0.0;
             -1.0 -1.0  0.0  0.0 0.0  0.0 -1.0 0.0  0.0;
              0.0  0.0  0.0  0.0 0.0  1.0  0.0 1.0 -1.0])
        @test cpxlp.b ≈ Float64[-10,10,5,-10,0]
        @test cpxlp.v_lb ≈ [0.0, 0.0, -0.9999999999999999, 0.0, 0.0, -Inf, 0.0, 0.0, 0.0]
        @test cpxlp.v_ub ≈ [Inf, Inf, 9.999999999999998, Inf, Inf, 4.999999999999999, Inf, Inf, Inf]
        @test cpxlp.c == [zeros(7); 1.; 0.]
        @test newB == [4,5,6,7,8]
    end
end


text:
module FewShotAnomalyDetection

using Flux
using NNlib
using Distributions
using SpecialFunctions
using Adapt
using Random
using LinearAlgebra
using EvalCurves
using Pkg
using Statistics

include("flux_extensions.jl")
include("KNNmemory.jl")
include("bessel.jl")
include("svae.jl")
include("svaebase.jl")
include("svaetwocaps.jl")
include("svae_vamp.jl")
include("svae_vamp_means.jl")
#include("vae.jl")

export SVAEbase, SVAEtwocaps, SVAEvamp, SVAEvampmeans, VAE, loss, wloss, log_pxexpectedz, pz, log_pz, log_px, log_pz_jacobian_encoder, log_pz_jacobian_decoder, samplez, zparams, printing_wloss, mem_wloss, log_det_jacobian_encoder, log_det_jacobian_decoder, hsplit1softp

end


text:
# ---
# title: 859. Buddy Strings
# id: problem859
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/buddy-strings/description/>
# hidden: true
# ---
#
# Given two strings `A` and `B` of lowercase letters, return `true` _if you can
# swap two letters in_`A` _so the result is equal to_`B` _, otherwise,
# return_`false` _._
#
# Swapping letters is defined as taking two indices `i` and `j` (0-indexed) such
# that `i != j` and swapping the characters at `A[i]` and `A[j]`. For example,
# swapping at indices `0` and `2` in `"abcd"` results in `"cbad"`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = "ab", B = "ba"
#     Output: true
#     **Explanation** **:** You can swap A[0] = 'a' and A[1] = 'b' to get "ba", which is equal to B.
#
#
# **Example 2:**
#
#
#
#     Input: A = "ab", B = "ab"
#     Output: false
#     **Explanation** **:** The only letters you can swap are A[0] = 'a' and A[1] = 'b', which results in "ba" != B.
#
#
# **Example 3:**
#
#
#
#     Input: A = "aa", B = "aa"
#     Output: true
#     **Explanation** **:** You can swap A[0] = 'a' and A[1] = 'a' to get "aa", which is equal to B.
#
#
# **Example 4:**
#
#
#
#     Input: A = "aaaaaaabc", B = "aaaaaaacb"
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: A = "", B = "aa"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= A.length <= 20000`
#   * `0 <= B.length <= 20000`
#   * `A` and `B` consist of lowercase letters.
#
#
## @lc code=start
using LeetCode

function buddy_strings(A::String, B::String)::Bool
    if length(A) == length(B)
        different_pairs = Pair{Char,Char}[]
        for (a, b) in zip(A, B)
            if a != b
                push!(different_pairs, a => b)
            end
        end
        if length(different_pairs) == 0
            ## check if A contains at least 2 duplicate letters
            length(Set(A)) < length(A)
        elseif length(different_pairs) == 2
            ((x, y), (z, w)) = different_pairs
            x == w && y == z
        else
            false
        end
    else
        false
    end
end
## @lc code=end


text:
using Documenter
using RandomisationModels

makedocs(
    modules = [RandomisationModels],
    sitename = "RandomisationModels.jl",
    authors = "James Totterdell",
    pages = [
        "Home" => "index.md",
        "Models" => ["models/multi-arm.md"],
        "Library" => ["lib/public.md"],
    ],
    format = Documenter.HTML(prettyurls = false),
)

deploydocs(repo = "github.com/jatotterdell/RandomisationModels.jl.git", devbranch = "main")


text:
module LoremJulia

include("LoremText.jl")
include("LoremFunctions.jl")

end # module

max_stars_repo_name: ColinFay/LoremJulia

max_issues_repo_name: ColinFay/LoremJulia

max_forks_repo_name: ColinFay/LoremJulia


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using HalfEdges
using HalfEdges: halfedges, vertices, vertex, edges, edge, halfedge, face, head, tail, opposite, next
using HalfEdges: orientation, polygons
using Serialization
using StaticArrays

∂(f::Function, y...) = (z...)->f(y...,z...)
second(c) = c[2]
he_ = HalfEdges

@testset "one tri" begin
  topo = HalfEdges.Topology([1,2,3])
  @test 1 == 1  # just make sure it didn't fail out
  @test vertex(HalfEdge(topo, opposite(topo, HalfEdgeHandle(topo, VertexHandle(1))))) == 1
  @test tail(topo, vertices(topo)[1])  == VertexHandle(1)
  @test he_.mapedge( i->i, topo)[1] == edges(topo)[1]
  @test he_.vertex_indices(topo) == 1:nverts(topo)
  @test he_.face_indices(topo) == 1:nfaces(topo)
  @test sort(vertices(topo, HalfEdgeHandle(1))) == sort(vertices(topo, FaceHandle(1)))
  @test length(halfedges(Polygon(topo, FaceHandle(1)))) == 3
  @test eltype(OneRing(topo,VertexHandle(1))) == HalfEdgeHandle
  @test HalfEdges.verts(topo) == vertices(topo)
  @test HalfEdges.edge_indices(topo) == 1:nedges(topo)

  @test isboundary(topo, nothing)
  @test face(topo, nothing) == nothing

  P = [[0.0,0.0,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0]]
  hface = opposite(topo, HalfEdgeHandle(topo, VertexHandle(1)))
  @test isboundary(topo, hface) == false
  @test he_.circumcenter_triangle(topo, P, hface) == [0.5,0.5,0.0]
  @test he_.normals(topo, P) == [[0.0,0.0,1.0]]
  @test facelist(topo) == [[1,2,3]]
  @test edge(topo, EdgeHandle(1)) == (1,2)
  @test sort(he_.angles(topo, P, Polygon(topo, FaceHandle(1)))) ≈ [π/4, π/4, π/2]

  # HalfEdgeHandle for vertex 3 is the most counterclockwise handle with 3 as it's tail
  @test isboundary(topo, vertices(topo)[3])
  @test edge(topo, vertices(topo)[3]) == (VertexHandle(3), VertexHandle(2))
  @test he_.angle(topo, P, vertices(topo)[3]) ≈ π/2
  @test he_.edgehandle(topo, first(UniqueHalfEdges(topo))) ==
            he_.edgehandle(topo, opposite(topo, first(UniqueHalfEdges(topo))))

  @test he_.area(topo, P, FaceHandle(1)) == 0.5
  @test isboundary(topo, VertexHandle(1))
end

@testset "bounded ring" begin
  topo = Topology([[1,2,3],[3,2,4],[4,5,3]],5)

  for i in 1:length(topo.v2he)
    heh = topo.v2he[i]
    @test heh == HalfEdges.opposite(topo, HalfEdges.opposite(topo, heh))
  end

  @test length(collect(HalfEdges.OneRing(topo,VertexHandle(1)))) == 2
  @test sort(OneRingVerts(topo,VertexHandle(1))) == [2, 3]
  @test length(collect(HalfEdges.OneRing(topo,VertexHandle(3)))) == 4
  @test length(collect(HalfEdges.OneRing(topo,VertexHandle(4)))) == 3
  ring3 = [HalfEdges.edge(topo,ih) for ih in HalfEdges.OneRing(topo,VertexHandle(3))]
  @test reduce( (acc,iv)->acc && iv == VertexHandle(3), [HalfEdges.tail(topo,ih) for ih in HalfEdges.OneRing(topo,VertexHandle(3))][1:end], init=true)
  # testing assumption for BoundaryLoop test, not a "real" test ( just assuming hedge 1 is boundary)
  @test HalfEdges.isboundary(topo, HalfEdgeHandle(1)) ||
        HalfEdges.isboundary(topo, HalfEdges.opposite(topo,HalfEdgeHandle(1)))
  @test HalfEdges.BoundaryLoop(topo,HalfEdgeHandle(1)) |> collect |> length == 5
  @test nhalfedges(topo) == length(edges(topo))*2
  alledge = edges(topo); allface = polygons(topo);
  @test orientation(allface[1], alledge[1]) == 1.0
  @test orientation(allface[1], (VertexHandle(2), VertexHandle(3))) == 1.0
  @test orientation(allface[2], (VertexHandle(2), VertexHandle(3))) == -1.0
  @test edges(topo)[IncidentEdges(topo, FaceHandle)[1]] == map(sort, edges(topo, FaceHandle(1)))
  @test edges(topo)[IncidentEdges(topo, FaceHandle)[2]] == map(sort, edges(topo, FaceHandle(2)))
  @test edges(topo)[IncidentEdges(topo, FaceHandle)[3]] == map(sort, edges(topo, FaceHandle(3)))
  @test map(heh->edge(topo,heh), IncidentHalfEdges(topo, FaceHandle)[3]) == map(sort, edges(topo, FaceHandle(3)))
end

@testset "unbroken ring" begin
  topo = Topology([[1,2,3],[1,3,4],[1,4,2]],4)
  P = [[0.0,0.0,0.0],[1.0,-0.25,0.0],[0.0,1.0,0.0],[-1.0,-0.25,0.0]]

  @test length(collect(HalfEdges.OneRing(topo,VertexHandle(1)))) == 3
  @test length(collect(HalfEdges.OneRing(topo,VertexHandle(4)))) == 3
  @test reduce( (acc,iv)->acc && iv == VertexHandle(1), [HalfEdges.tail(topo,ih) for ih in HalfEdges.OneRing(topo,VertexHandle(1))][1:end], init=true)
  @test nhalfedges(topo) == length(edges(topo))*2

  FE = he_.incidence(topo, EdgeHandle, FaceHandle,true)
  EV = he_.incidence(topo, VertexHandle, EdgeHandle,true)
  @test iszero(FE*EV)
  @test length(vertexnormals(topo, P)) == nvertices(topo)
  @test vertexnormals(topo,P)[1] == [0.0, 0.0, 1.0]
  @test length(normals(topo, P)) == nfaces(topo)
  @test normals(topo,P)[1] == [0.0, 0.0, 1.0]
end

@testset "3x3grid" begin
  cell(i,j) = map(v->v+(j-1)*4, [i,i+4,i+1,i+1,i+4,i+5])
  topo = Topology(Iterators.flatten([cell(1,1), cell(2,1), cell(3,1),
                                     cell(1,2), cell(2,2), cell(3,2),
                                     cell(1,3), cell(2,3), cell(3,3)])|>collect)
  @test sort(first(he_.boundary_interior(topo))) == [6,7,10,11]

end

@testset "box48" begin
  f = open(string(Base.@__DIR__, "/box48.spork"), "r")
  mesh = deserialize(f)
  close(f)

  hemesh = Topology(mesh["Poly"], length(mesh["P"]))
  P = mesh["P"]

  @test sum(map(i->he_.area(hemesh,P,i), he_.faces(hemesh))) == 6.0
  @test he_.nedges(hemesh) == he_.edges(hemesh) |> length
  @test he_.nhalfedges(hemesh)/2 == he_.nedges(hemesh) == he_.length(he_.edges(hemesh))
  @test he_.edges(hemesh) |> collect == ∂(edge,hemesh).(he_.UniqueHalfEdges(hemesh))
end

@testset "incidence" begin
  topo = Topology([[1,2,3],[2,4,3]])

  @test edges(topo)[IncidentEdges(topo, FaceHandle)[1]] == map(sort, edges(topo, FaceHandle(1)))
  @test edges(topo)[IncidentEdges(topo, FaceHandle)[2]] == map(sort, edges(topo, FaceHandle(2)))

  P = [[0.0,0.0,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0], [0.5, 0.5,1.0]]
  @test sort(map(h->he_.dihedral_angle(topo, P, h), UniqueHalfEdges(topo)))[[1,end]] ≈ [0.0,π/2]
  # not really working so well, but we will test it anyways
  P = [[0.25,0.25,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0], [0.75, 0.75,0.0]]
  @test he_.improve_mesh(topo, P) > 0
end

@testset "readme" begin
  topo = Topology([[1,2,3],[1,4,2]])

  @test edges(topo) |> collect == ∂(edge, topo).(UniqueHalfEdges(topo))

  # find the only halfedge pointing from vertex 1 to 2
  h = Iterators.filter(h->head(topo, h)==VertexHandle(2), OneRing(topo, VertexHandle(1))) |> first
  @test h == halfedge(topo,(1,2))
  @test tail(topo, h) == VertexHandle(1)
  @test head(topo, h) == VertexHandle(2)
  @test vertex(topo, h) == VertexHandle(2)

  h = opposite(topo, h)
  @test h == halfedge(topo, (2,1))

  topo = Topology([[1,2,3],[1,4,2],[1,3,5,4],[2,4,6]])
  @test length(edges(topo, OneRing(topo, VertexHandle(1)))) == 3
  @test length(edges(topo, Rim(topo, VertexHandle(1)))) == 4
  @test length(edges(topo, BoundaryLoop(topo, first(OneRing(topo, VertexHandle(2)))))) == 5
  @test length(edges(topo, Polygon(topo, FaceHandle(1)))) == 3
  @test length(edges(topo, Polygon(topo, FaceHandle(3)))) == 4

  v6 = VertexHandle(6)
  @test head(topo, next(topo, opposite(topo, next(topo, next(topo, opposite(topo, next(topo, HalfEdgeHandle(1)))))))) == v6
  pivot = opposite(next(next))

  @test head(next)(topo, pivot(topo, opposite(next)(topo, HalfEdgeHandle(1)))) == v6

  @test (halfedge |> next |> opposite |> next |> next |> opposite |> next |> head)(topo, (3,1)) == v6
end

@testset "Degenerate" begin
  topo = Topology([2,3,4])
  P = map(i->rand(SVector{3,Float64}), 1:4)
  @test first(vertexnormals(topo, P)) == zero(SVector{3,Float64})
  @test isempty(floodfill(topo, P)) == false
end

@HandleType Foo
@HandleType Bar
@HandleType Baz

hey(x::Integer) = x+1
hey(x::Foo) = string("foo", x)
hey(x::Bar) = string("bar", x)

@testset "Handles" begin
  @test Foo(1) == 1

  a = sort(map( f->f+1 + Foo(1), (2 .* Foo.([10:-1:-2...]))))
  @test eltype(a) == Foo
  @test a[1] == -2 == Foo(-2)
  @test a[2] < a[3]
  @test a[end] == 22

  @test hey(Foo(2)) == "foo2"
  @test hey(Bar(3)) == "bar3"
  @test hey(Baz(4)) == 5

end

function cube()
  P = [0 0 0
       1 0 0
       0 0 1
       1 0 1
       1 1 1
       0 1 1
       1 1 0
       0 1 0
      ]
  P = SVector{3}.(eachrow(Float64.(P)))
  (Topology([1,2,3, 2,4,3, 3,4,6, 4,5,6, 2,5,4, 2,7,5, 7,8,5, 5,8,6, 1,3,6, 1,6,8, 1,8,2, 2,8,7]), P)
end

@testset "collide" begin
  topo,P = cube()
  col = Collider(topo, P)
  tophalf = query_aabb(col, HalfEdges.BVH.AABB(SVector{3}(0.,0.,0.5), SVector{3}(1.0,1.0,1.0)))
  pre = []
  inorder = []
  post = []
  HalfEdges.BVH.visitdf(col.bvh, x->push!(pre, x), x->push!(post, x), x->push!(inorder, x))
  @test length(pre) == length(post) == length(inorder) > 0

  # should contain all triangles other than bottom two
  @test length(tophalf) == nfaces(topo)-2
  @test isempty(tophalf ∩ [11, 12])

  #pull point into self intersection
  col.mesh.P[1] = SVector{3}(0.5,0.5,2.0)
  update_collider(col)
  #foo(a1,a2,a3, b1,b2,b3) = a2[1] > 1.0
  #@test collide_self(col, foo) == collide_self(topo, col.mesh.P, foo)
  @test collide_self(col) == collide_self(topo, col.mesh.P)

  #==
      2---4
     /|\  |\
    5 | \ | 6    then move 5 and 6 so they intersect
     \|  \|/
      1---3
  ==#
  F = [[1,3,2],[2,3,4],[2,5,1],[3,6,4]]
  topo = Topology(F)
  P = [0 0 0
       0 1 0
       1 0 0
       1 1 0
       1 0.5 1
       0 0.5 1
      ]
  P = SVector{3}.(eachrow(Float64.(P)))
  #foo(a1,a2,a3, b1,b2,b3) = a1[2] > 0.0   # our two triangle with z > 0 are in collision
  #hits = collide_self(Collider(topo, P), foo)
  hits = collide_self(Collider(topo, P))
  @test length(hits) == 1

  # want to get lots of hits to test when hits array is resized
  topo = Topology(repeat(F, 100))
  #hits = collide_self(Collider(topo, P), foo)
  hits = collide_self(Collider(topo, P))
  @test length(hits) > 1
end

@testset "winding numbers" begin
  topo, P = cube()

  Fi = FaceHandle.(1:nfaces(topo))
  centre = sum(P)/length(P)
  @test sum(map(fh->he_.solid_angle(topo, P, he_.halfedge(topo, fh), centre), Fi)) ≈ 4.0*π
  @test sum(map(fh->he_.solid_angle(topo, P, he_.halfedge(topo, fh), P[1]), Fi)) ≈ π/2.0

  # outside, accuracy seems lower for some reason
  @test abs(sum(map(fh->he_.solid_angle(topo, P, he_.halfedge(topo, fh), P[1]-centre), Fi))) < 1e6

  cache = winding_number_cache(topo, P, 1)
  #cache3 = winding_number_cache(topo, P, 1, 3)

  @test winding_number(topo, P, SVector{3}(-1.0,-1.0,-1.0), cache) |> round == 0
  #@test winding_number(topo, P, SVector{3}(-1.0,-1.0,-1.0), cache3) |> round == 0
  @test winding_number(topo, P, SVector{3}(-1.0,-1.0,-1.0), cache.bvh) |> round == 0
  @test winding_number(topo, P, SVector{3}(-1.0,-1.0,-1.0)) |> round == 0
  @test winding_number(topo, P, SVector{3}(0.1,0.2,0.5), cache) |> round == 1
  @test winding_number(topo, P, SVector{3}(0.1,0.2,0.5)) |> round == 1

  @test HalfEdges.isorphan(floodfill(topo, P;verbose=true), 1) == false
  @test floodfill(topo, P) |> length == 1
  topo = Topology(vcat(polygons(topo), [[9, 10, 11]]))
  P = vcat(P, [SVector(0.5,3.0,-0.5), SVector(0.5,-1.0,3.5), SVector(0.5,-1.0,-1.0)])
  @test floodfill(topo, P) |> length == 3
end

@testset "load a mesh" begin
  topo, P = loadmesh("bunnylow.obj")
  @test typeof(topo) == Topology
  @test nfaces(topo) > 0
end

@testset ".Collision" begin
  Point = SVector{3,Float64}
  o,x,y,z = Point(0.,0,0), Point(1.,0,0), Point(0.,1,0), Point(0.,0,1)
  a,b = HalfEdges.Collision.segment_between(((o, x, y),
                                       (0.2x+0.2y-0.5z, 0.2x+0.2y+0.5z, 0.2x - 0.5y+0.5z)),
                                       [(2, (1, 2)), (1, (1, 2))])
  @test a[3] == b[3] == 0.0
  @test HalfEdges.Collision.triangle_edges(o, x, y, Point(0.2, 0.2, 1.0),
                                              Point(0.2, 0.2, -1.0),
                                              Point(0.1, 0.1, -1.0)) |> first == true
  @test HalfEdges.Collision.triangle_edges(o, x, y, Point(0.2, 0.2, -0.1),
                                              Point(0.2, 0.2, -1.0),
                                              Point(0.1, 0.1, -1.0)) |> first == false
end

@testset "DifferentialGeometry" begin
  topo, P = cube()

  P = HalfEdges.poisson(HalfEdges.Δ(topo, P), zero.(P), P, [1,2,3,4])
  # minimal surface on x-z plane
  @test mapreduce(p->p[2], +, P) == 0.0
  @test mapreduce(p->p[1], *, P[5:end]) > 0.0
  @test mapreduce(p->p[1], *, P[5:end]) < 1.0

end

@testset "IslandFinder" begin
  @test length(find_islands([[1,2],[3,4],[5,6],[7,8]])) == 4
  @test length(find_islands([[1,2],[3,4],[5,6],[6,1]])) == 2
  @test length(find_islands([[1,2],[3,4],[5,6],[6,1],[4,5]])) == 1
end

include("BoundingVolumeTrees_runtests.jl")


text:
using JuMP
using Gurobi
using Plasmo
using PlasmoAlgorithms

mp = Model(solver = GurobiSolver())
sp1 = Model(solver = GurobiSolver())
sp2 = Model(solver = GurobiSolver())
sp3 = Model(solver = GurobiSolver())

@variable(mp,x[1:2]>=0)
@constraint(mp,x[1]+x[2]>=50)
@objective(mp,Min,x[1]+2x[2])

@variable(sp1,x[1:2]>=0)
@variable(sp1,y[1:2]>=0)
@constraint(sp1, x[1]+y[1]>=60)
@constraint(sp1, x[2]+y[2]>=70)
@objective(sp1, Min, 3y[1]+4y[2])

@variable(sp2,y[1:2]>=0)
@variable(sp2, z[1:2]>=0)
@constraint(sp2, y[1]+z[1]>=80)
@constraint(sp2, y[2]+z[2]>=90)
@objective(sp2, Min, 5z[1]+6z[2])

@variable(sp3,x[1:2]>=0)
@variable(sp3,t[1:2]>=0)
@constraint(sp3,x[1]+t[2]>=43)
@constraint(sp3,x[2]+t[1]>=75)
@objective(sp3,Min,.5t[1]+5t[2])

g = ModelGraph()
setsolver(g, GurobiSolver())
n1 = add_node(g)
n2 = add_node(g)
n3 = add_node(g)
n4 = add_node(g)

setmodel(n1, mp)
setmodel(n2, sp1)
setmodel(n3, sp2)
setmodel(n4, sp3)

edge1 = Plasmo.add_edge(g,n1,n2)
edge2 = Plasmo.add_edge(g,n2,n3)
edge2 = Plasmo.add_edge(g,n1,n4)

@linkconstraint(g,[i in 1:2], n1[:x][i] == n2[:x][i])
@linkconstraint(g,[i in 1:2], n1[:x][i] == n4[:x][i])
@linkconstraint(g,[i in 1:2], n2[:y][i] == n3[:y][i])

bendersolve(g,max_iterations =10)


text:
module Pointers

    using ..Misc

    export Pointer, shift_start

    """
        Pointer(arr::AbstractArray)

    Returns the pointer for the argument array.
    It also includes auxiliar functions for pointer load and store.

    Usage
    ```julia-repl
    julia> a = [1,2,3,4]
    julia> pt = Pointer(a)
    julia> pt[3]
    3
    ```
    """
    Pointer(arr::AbstractArray) = pointer(arr)

    """
        Pointer(arr::AbstractArray, el_shift::Integer)

    Same as Pointer, but shifts the array in el_shift number of
    elements.
    """
    Pointer(arr::AbstractArray, el_shift::Integer) =
        Pointer(arr) << el_shift

    @inline Base.getindex(ptr::Ptr, i::Integer) = unsafe_load(ptr, i)
    @inline Base.setindex!(ptr::Ptr{T}, val::T, i::Integer) where{T} =
        unsafe_store!(ptr, val, i)

    """
        shift_start(ptr::Ptr, i::Integer)

    Shifts the pointer in the specified i number of elements.
    Another way to call is to use `ptr << nof_elements`
    Usage:
    ```julia-repl
    julia> a = [1,2,3,4]
    julia> pt = Pointer(a)
    julia> pt[1]
    1
    julia> pt = pt << 2
    julia> pt[1]
    3
    ```
    """
    @inline shift_start(ptr::Ptr{T}, i::Integer) where {T} = ptr + sizeof(T) * i
    """
        >>(sl::Ptr{T}, i::Integer)

    Same as shift_start.
    """
    Base.:<<(ptr::Ptr{T}, i::Integer) where {T} = shift_start(ptr, i)

end


text:
using HPAT
using MPI

@acc hpat function calcNaiveBayes(file_name)
    num_classes = 20
    points = DataSource(Matrix{Float64},HDF5,"/points", file_name)
    labels = DataSource(Vector{Float64},HDF5,"/labels", file_name)
    coeffs = HPAT.NaiveBayes(points, labels, num_classes)
    return coeffs
end

c2 = calcNaiveBayes(ENV["SCRATCH"]*"/benchmark_data/naivebayes_train_small.hdf5")
rank = MPI.Comm_rank(MPI.COMM_WORLD)
pes = MPI.Comm_size(MPI.COMM_WORLD)

MPI.Barrier(MPI.COMM_WORLD)
t1 = time_ns()
c2 = calcNaiveBayes(ENV["SCRATCH"]*"/benchmark_data/naivebayes_train_large.hdf5")
t2 = time_ns()

if rank==0
	println("nodes: ", pes, "\nnaive bayes exec time: ", (t2-t1)/1.0e9, "\nCoeffs: ", c2)
end

MPI.Barrier(MPI.COMM_WORLD)


text:
using Images, JLD

function main()
    w = 2000
    h = 2000
    points = load("./sets/15k.jld", "points")
    img = zeros(RGB,w,h)
    for c::ComplexF64 in points
        i::Int = 0
        z::ComplexF64 = c
        while abs(z) < 3 && i < 5000
            x::Int = trunc(Int,(real(z) + 2)*w/3); y::Int = trunc(Int,(imag(z) + 1.5)*h/3)
            0<x<w && 0<y<h && i > 0 ? img[x,y] += RGB(0.0005,0,0) : nothing #log10(i)/3
            z = z*z + c
            i += 1
        end
    end
    return map(clamp01nan, img)
end

@time save("./out/buddhared.png", main())

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
################################################################################
#                             sleepstudy.csv
################################################################################

df        = CSV.File(joinpath(path, "csv", "sleepstudy.csv"); types = [String, Float64, String, String]) |> DataFrame

@testset "  Model 1: sleepstudy.csv SI/SI                            " begin
    #=
    SPSS
    REML 1729.492560
    Random 1375.465318
    Residual 987.588488
    =#
    # Model 1
    lmm = Metida.LMM(@formula(Reaction~Days), df;
    random = Metida.VarEffect(Metida.@covstr(1|Subject), Metida.SI),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 1729.4925602367025 atol=1E-6
end

@testset "  Model 2: sleepstudy.csv CS/SI                            " begin
    # REML SPSS 1904.327
    # 1662.172084
    # 296.693108
    # Model 2
    lmm = Metida.LMM(@formula(Reaction~1), df;
    random = Metida.VarEffect(Metida.@covstr(Days|Subject), Metida.CS),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 1904.3265170722132 atol=1E-6

    #=
    lmm = Metida.LMM(@formula(Reaction~1 + Days), df;
    repeated = Metida.VarEffect(Metida.@covstr(1), Metida.CS, subj = :Subject)
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 1729.4925602367027 atol=1E-6
    =#
end

@testset "  Model 3: sleepstudy.csv CSH/SI                           " begin
    # REML 1772.095
    # Model 3
    lmm = Metida.LMM(@formula(Reaction~1), df;
    random = Metida.VarEffect(Metida.@covstr(Days|Subject), Metida.CSH)
    )
    Metida.fit!(lmm; init = [26.4881, 35.5197, 34.8287, 56.1999, 63.8281, 85.4346, 94.218, 92.8584, 113.679, 129.721, 0.959643, 22.5597])
    @test lmm.result.reml ≈ 1772.0953251997046 atol=1E-6
end

@testset "  Model 4: sleepstudy.csv ARH/SI                           " begin
    # REML 1730.189543
    # Model 4
    lmm = Metida.LMM(@formula(Reaction~1), df;
    random = Metida.VarEffect(Metida.@covstr(Days|Subject), Metida.ARH)
    )
    Metida.fit!(lmm; init = [37.9896, 41.1392, 34.1041, 48.1435, 52.2191, 72.4237, 83.3405, 76.7782, 90.2571, 102.617, 0.900038, 6.83327])
    @test lmm.result.reml ≈ 1730.1895427398322 atol=1E-6
end

################################################################################
#                             Pastes.csv
################################################################################

df        = CSV.File(path*"/csv/Pastes.csv"; types = [String, Float64, String, String, String]) |> DataFrame
@testset "  Model 5: Pastes.csv SI,SI/SI                             " begin
    # REML 246.990746
    # Model 5
    lmm = Metida.LMM(@formula(strength~1), df;
    random = [Metida.VarEffect(Metida.@covstr(1|batch), Metida.SI),
    Metida.VarEffect(Metida.@covstr(1|batch&cask), Metida.SI)]
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 246.99074585348623 atol=1E-6
end

@testset "  Model 6: Pastes.csv ARMA/SI                              " begin
    # SPSS REML 246.818951
    # Model 6
    lmm = Metida.LMM(@formula(strength~1), df;
    random = Metida.VarEffect(Metida.@covstr(cask|batch), Metida.ARMA),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 246.81895071012508 atol=1E-6
end

################################################################################
#                             Penicillin.csv
################################################################################

df        = CSV.File(path*"/csv/Penicillin.csv"; types = [String, Float64, String, String]) |> DataFrame
@testset "  Model 7: Penicillin.csv SI,SI/SI                         " begin
    # SPSS 330.860589
    # Model 7
    lmm = Metida.LMM(@formula(diameter ~ 1), df;
    random = [Metida.VarEffect(Metida.@covstr(1|plate), Metida.SI),
    Metida.VarEffect(Metida.@covstr(1|sample), Metida.SI)],
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 330.86058899109184 atol=1E-6

    #fm = @formula(diameter ~ 1 + (1|plate) + (1|sample))
    #mm  = fit(MixedModel, fm, df, REML=true)

    #=
    #SPSS 432.686
    lmm = Metida.LMM(@formula(diameter~1), df;
    random = Metida.VarEffect(Metida.@covstr(plate|sample), Metida.ARMA)
    )
    Metida.fit!(lmm)
    =#
end

################################################################################
#                             RepeatedPulse.csv
################################################################################

df        = CSV.File(path*"/csv/RepeatedPulse.csv"; types = [String, Float64, String, String]) |> DataFrame
df.Pulse = float.(df.Pulse)
sort!(df, :Day)

@testset "  Model 8: RepeatedPulse.csv SI/AR                         " begin
    #=
    SPSS
    REML 453.339544
    =#
    lmm = Metida.LMM(@formula(Pulse~1), df;
    random = Metida.VarEffect(Metida.@covstr(Time|Time), Metida.SI),
    repeated = Metida.VarEffect(Metida.@covstr(Day|Time), Metida.AR),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 453.3395435627574 atol=1E-6
end
@testset "  Model 9: RepeatedPulse.csv 0/AR                          " begin
    #=
    SPSS
    REML 471.851077
    =#
    lmm = Metida.LMM(@formula(Pulse~1), df;
    repeated = Metida.VarEffect(Metida.@covstr(Day|Time), Metida.AR),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 471.85107712169827 atol=1E-6
end
@testset "  Model 10: RepeatedPulse.csv AR/SI                        " begin
    #=
    SPSS
    REML 453.339555
    =#
    lmm = Metida.LMM(@formula(Pulse~1), df;
    random = Metida.VarEffect(Metida.@covstr(Day|Time), Metida.AR),
    )
    Metida.fit!(lmm)
    #@test lmm.result.reml ≈ 453.3395560121246 atol=1E-6
    @test lmm.result.reml ≈ 453.3395546358513 atol=1E-6
end

################################################################################
#                             ChickWeight.csv
################################################################################
#=
SPSS
REML 5456.638120
=#
df        = CSV.File(path*"/csv/ChickWeight.csv"; types = [String, Float64, Float64, String, String]) |> DataFrame
sort!(df, :Diet)
@testset "  ARH ChickWeight.csv                                      " begin
    lmm = Metida.LMM(@formula(weight~1 + Diet & Time), df;
    random = Metida.VarEffect(Metida.@covstr(1|Chick), Metida.SI),
    repeated = Metida.VarEffect(Metida.@covstr(Diet|Chick), Metida.ARH),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 4439.254893054906 atol=1E-6

    lmm = Metida.LMM(@formula(weight~1 + Diet & Time), df;
    random = Metida.VarEffect(Metida.@covstr(Diet|Chick), Metida.ARH),
    )
    Metida.fit!(lmm)
    @test lmm.result.reml ≈ 5480.751465914058 atol=1E-6
end