text:
import Base.Threads.Atomic

mutable struct Controller
  buttons::Atomic{UInt8}
  index::UInt8
  strobe::UInt8
  Controller() = new(Atomic{UInt8}(0), 0, 0)
end

function read(c::Controller)::UInt8
  val = 0x00
  if c.index < 0x08 && (c.buttons[] & (0x01 << c.index)) != 0x00
    val = 0x01
  end
  c.index += 0x01
  if (c.strobe & 0x01) == 0x01
    c.index = 0x00
  end
  val
end

function write!(c::Controller, val::UInt8)
  c.strobe = val
  if (c.strobe & 0x01) == 0x01
    c.index = 0x00
  end
end

ext: jl


text:
import Base:
    convert,
    done,
    eltype,
    find,
    get,
    getindex,
    isempty,
    keys,
    length,
    linearindexing,
    next,
    setindex!,
    show,
    size,
    start,
    summary,
    values,
    vcat,
    >>,
    .==,
    .!=,
    .>,
    .>=,
    .<,
    .<=,
    .+,
    .-,
    .*,
    ./,
    .÷,
    ~,
    &,
    |

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# level 1 - generic routines

function broutine2x2!(st::AbstractVecOrMat, U::AbstractMatrix, locs::Locations)
    U11 = U[1, 1]; U12 = U[1, 2];
    U21 = U[2, 1]; U22 = U[2, 2];
    for loc in locs
        _broutine2x2!(st, (U11,U12,U21,U22), loc)
    end
    return st
end

# NOTE:
# wrapping the following code seems to have a small overhead
# so we are not wrapping them in a function here
function broutine2x2_kernel_expr(idx_1, idx_2)
    return quote
        ST1 = U11 * st[$idx_1] + U12 * st[$idx_2]
        ST2 = U21 * st[$idx_1] + U22 * st[$idx_2]

        st[$idx_1] = ST1
        st[$idx_2] = ST2
    end
end

function broutine2x2_m_kernel_expr(idx_1, idx_2)
    return quote
        ST1 = U11 * st[b, $idx_1] + U12 * st[b, $idx_2]
        ST2 = U21 * st[b, $idx_1] + U22 * st[b, $idx_2]

        st[b, $idx_1] = ST1
        st[b, $idx_2] = ST2
    end
end

# ST1 = U11 * st[j+i] + U12 * st[j+i+step_1]
# ST2 = U21 * st[j+i] + U22 * st[j+i+step_1]
# st[j+i] = ST1
# st[j+i+step_1] = ST2
body_broutine2x2 = quote
    step_1 = 1 << (plain(locs) - 1)
    step_2 = 1 << plain(locs)

    @inbounds if step_1 == 1
        for j in 0:step_2:size(st, 1)-step_1
            $(broutine2x2_kernel_expr(:(j+1), :(j+2)))
        end
        return st
    end
end

body_broutine2x2_m = quote
    step_1 = 1 << (plain(locs) - 1)
    step_2 = 1 << plain(locs)

    @inbounds if step_1 == 1
        for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
            $(broutine2x2_m_kernel_expr(:(j+1), :(j+2)))
        end
        return st
    end
end

init_broutine2x2_ctrl = quote
    raw = first(plain(locs))
    raw_ctrl = plain(ctrl)
    step_1 = 1 << (raw - 1)
    step_2 = 1 << raw
    ctrl_mask = bmask(raw_ctrl)
    flag_mask = reduce(+, one(Int) << (raw_ctrl[i] - 1) for i in 1:length(ctrl) if ctrl.flags[i])
end

body_broutine2x2_ctrl = quote
    $init_broutine2x2_ctrl
    @inbounds if step_1 == 1
        for j in 0:step_2:size(st, 1)-step_1
            if ismatch(j, ctrl_mask, flag_mask)
                $(broutine2x2_kernel_expr(:(j+1), :(j+2)))
            end
        end
        return st
    end
end

body_broutine2x2_ctrl_m = quote
    $init_broutine2x2_ctrl
    @inbounds if step_1 == 1
        for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
            if ismatch(j, ctrl_mask, flag_mask)
                $(broutine2x2_m_kernel_expr(:(j+1), :(j+2)))
            end
        end
        return st
    end
end

# case 1: we expand the inner loop of step_1 if it's smaller than 16
for k in [2,4,8,16]
    push!(body_broutine2x2.args, quote
        @inbounds if step_1 == $k
            for j in 0:step_2:size(st, 1)-step_1
                Base.Cartesian.@nexprs $k i-> begin
                    $(broutine2x2_kernel_expr(:(j+i), :(j+i+step_1)))
                end
            end
            return st
        end
    end)

    push!(body_broutine2x2_m.args, quote
        @inbounds if step_1 == $k
            for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
                Base.Cartesian.@nexprs $k i-> begin
                    $(broutine2x2_m_kernel_expr(:(j+i), :(j+i+step_1)))
                end
            end
            return st
        end
    end)

    push!(body_broutine2x2_ctrl.args, quote
        @inbounds if step_1 == $k
            for j in 0:step_2:size(st, 1)-step_1
                Base.Cartesian.@nexprs $k i->begin
                    if ismatch(j+i-1, ctrl_mask, flag_mask)
                        $(broutine2x2_kernel_expr(:(j+i), :(j+i+step_1)))
                    end
                end
            end
            return st
        end
    end)

    push!(body_broutine2x2_ctrl_m.args, quote
        @inbounds if step_1 == $k
            for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
                Base.Cartesian.@nexprs $k i->begin
                    if ismatch(j+i-1, ctrl_mask, flag_mask)
                        $(broutine2x2_m_kernel_expr(:(j+i), :(j+i+step_1)))
                    end
                end
            end
            return st
        end
    end)
end

# case 2: we expand the inner loop by 16, if step_1 is larger than 16
push!(body_broutine2x2.args, quote
    for j in 0:step_2:size(st, 1)-step_1
        for i in j:16:j+step_1-1
            Base.Cartesian.@nexprs 16 k->begin
                $(broutine2x2_kernel_expr(:(i+k), :(i+step_1+k)))
            end
        end
    end
    return st
end)

push!(body_broutine2x2_m.args, quote
    for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
        for i in j:16:j+step_1-1
            Base.Cartesian.@nexprs 16 k->begin
                $(broutine2x2_m_kernel_expr(:(i+k), :(i+step_1+k)))
            end
        end
    end
    return st
end)

push!(body_broutine2x2_ctrl.args, quote
    for j in 0:step_2:size(st, 1)-step_1
        for i in j:16:j+step_1-1
            Base.Cartesian.@nexprs 16 k->begin
                if ismatch(i+k-1, ctrl_mask, flag_mask)
                    $(broutine2x2_kernel_expr(:(i+k), :(i+step_1+k)))
                end
            end
        end
    end
    return st
end)

push!(body_broutine2x2_ctrl_m.args, quote
    for j in 0:step_2:size(st, 2)-step_1, b in 1:size(st, 1)
        for i in j:16:j+step_1-1
            Base.Cartesian.@nexprs 16 k->begin
                if ismatch(i+k-1, ctrl_mask, flag_mask)
                    $(broutine2x2_m_kernel_expr(:(i+k), :(i+step_1+k)))
                end
            end
        end
    end
    return st
end)

# create broutine2x2 function expr
def_broutine2x2 = Dict{Symbol, Any}(
    :name => :_broutine2x2!,
    :args => [:(st::AbstractVector), :((U11,U12,U21,U22)), :(locs::Locations)],
    :body => body_broutine2x2
)

def_broutine2x2_m = Dict{Symbol, Any}(
    :name => :_broutine2x2!,
    :args => [:(st::AbstractMatrix), :((U11,U12,U21,U22)), :(locs::Locations)],
    :body => body_broutine2x2_m
)

def_broutine2x2_ctrl = Dict{Symbol, Any}(
    :name => :_broutine2x2!,
    :args => [:(st::AbstractVector), :((U11,U12,U21,U22)), :(locs::Locations), :(ctrl::CtrlLocations)],
    :body => body_broutine2x2_ctrl,
)

def_broutine2x2_ctrl_m = Dict{Symbol, Any}(
    :name => :_broutine2x2!,
    :args => [:(st::AbstractMatrix), :((U11,U12,U21,U22)), :(locs::Locations), :(ctrl::CtrlLocations)],
    :body => body_broutine2x2_ctrl_m,
)

@eval $(combinedef(def_broutine2x2))
@eval $(combinedef(def_broutine2x2_m))
@eval $(combinedef(def_broutine2x2_ctrl))
@eval $(combinedef(def_broutine2x2_ctrl_m))

function broutine2x2!(st::AbstractVecOrMat, U::AbstractMatrix, locs::Locations, ctrl::CtrlLocations)
    U11 = U[1, 1]; U12 = U[1, 2];
    U21 = U[2, 1]; U22 = U[2, 2];
    for loc in locs
        _broutine2x2!(st, (U11,U12,U21,U22), locs, ctrl)
    end
    return st
end

# level 2 - special matrix routines

# level 3 - intrinsic gate routines
@broutine X = [0 1;1 0]
@broutine Y = [0 -im;im 0]
@broutine Z = [1 0;0 -1]
@broutine H = [1/sqrt(2) 1/sqrt(2); 1/sqrt(2) -1/sqrt(2)]
@broutine T = [1 0;0 exp(im * π / 4)]
@broutine S = [1 0;0 im]
@broutine Sdag = [1 0;0 -im]
@broutine Tdag = [1 0;0 exp(-im * π / 4)]
@broutine Rx(θ::Real) = [cos(θ / 2) -im*sin(θ / 2); -im*sin(θ / 2) cos(θ / 2)]
@broutine Ry(θ::Real) = [cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)]
@broutine Rz(θ::Real) = [exp(-im * θ / 2) 0; 0 exp(im * θ / 2)]

@broutine PSWAP(θ) = [-im * sin(θ/2) + cos(θ/2) 0 0 0;
                0 cos(θ/2) -im * sin(θ/2) 0;
                0 -im * sin(θ/2) cos(θ/2) 0;
                0 0 0 -im * sin(θ/2) + cos(θ/2)]

# @broutine SWAP = [1 0 0 0;
#                 0 0 1 0;
#                 0 1 0 0;
#                 0 0 0 1]
#
# NOTE: the compiler seems to be happier without us optimizing the kernel function
#      don't remove T1 and T4, so let's just define this manually :(

const BSWAP = [1 0 0 0;
                0 0 1 0;
                0 1 0 0;
                0 0 0 1]

@inline function unsafe_swap_kernel(st, b3, step_2, step_4)
    @fastmath @inbounds begin
        i1 = b3+1
        i2 = b3+step_2+1
        i3 = b3+step_4+1
        i4 = b3+step_2+step_4+1

        T1 = st[i1]
        T2 = st[i3]
        T3 = st[i2]
        T4 = st[i4]

        st[i1] = T1
        st[i2] = T2
        st[i3] = T3
        st[i4] = T4
    end
end

@inline function unsafe_swap_kernel(st, b, b3, step_2, step_4)
    @fastmath @inbounds begin
        i1 = b3+1
        i2 = b3+step_2+1
        i3 = b3+step_4+1
        i4 = b3+step_2+step_4+1

        T1 = st[b, i1]
        T2 = st[b, i3]
        T3 = st[b, i2]
        T4 = st[b, i4]

        st[b, i1] = T1
        st[b, i2] = T2
        st[b, i3] = T3
        st[b, i4] = T4
    end
end

function _swap_step(locs)
    loc_1 = first(plain(locs))
    loc_2 = last(plain(locs))

    if loc_1 > loc_2
        loc_1 = last(plain(locs))
        loc_2 = first(plain(locs))
    end

    step_1 = 1 << loc_1
    step_2 = 1 << (loc_1 - 1)
    step_3 = 1 << loc_2
    step_4 = 1 << (loc_2 - 1)

    return step_1, step_2, step_3, step_4
end

function broutine!(st::AbstractVector, ::Val{:SWAP}, locs::Locations)
    @assert length(locs) == 2 "number of locations mismatch"
    step = _swap_step(locs)

    @fastmath @inbounds for b1 in 0:step[3]:size(st, 1)-step[3]
        for b2 in b1:step[1]:b1+step[4]-step[1]
            for b3 in b2:b2+step[2]-1
                unsafe_swap_kernel(st, b3, step[2], step[4])
            end
        end
    end
    return st
end

function broutine!(st::AbstractMatrix, ::Val{:SWAP}, locs::Locations)
    @assert length(locs) == 2 "number of locations mismatch"
    step = _swap_step(locs)

    @fastmath @inbounds for b1 in 0:step[3]:size(st, 2)-step[3]
        for b2 in b1:step[1]:b1+step[4]-step[1]
            for b3 in b2:b2+step[2]-1, b in axes(st, 1)
                unsafe_swap_kernel(st, b, b3, step[2], step[4])
            end
        end
    end
    return st
end

function broutine!(st::AbstractVector, ::Val{:SWAP}, locs::Locations, ctrl::CtrlLocations)
    @assert length(locs) == 2 "number of locations mismatch"
    step = _swap_step(locs)

    raw_ctrl = plain(ctrl)
    ctrl_mask = bmask(raw_ctrl)
    flag_mask = reduce(+, one(Int) << (raw_ctrl[i] - 1) for i in 1:length(ctrl) if ctrl.flags[i])

    @fastmath @inbounds for b1 in 0:step[3]:size(st, 1)-step[3]
        for b2 in b1:step[1]:b1+step[4]-step[1]
            for b3 in b2:b2+step[2]-1
                if ismatch(b3, ctrl_mask, flag_mask)
                    unsafe_swap_kernel(st, b3, step[2], step[4])
                end
            end
        end
    end
    return st
end

function broutine!(st::AbstractMatrix, ::Val{:SWAP}, locs::Locations, ctrl::CtrlLocations)
    @assert length(locs) == 2 "number of locations mismatch"
    step = _swap_step(locs)

    raw_ctrl = plain(ctrl)
    ctrl_mask = bmask(raw_ctrl)
    flag_mask = reduce(+, one(Int) << (raw_ctrl[i] - 1) for i in 1:length(ctrl) if ctrl.flags[i])

    @fastmath @inbounds for b1 in 0:step[3]:size(st, 2)-step[3]
        for b2 in b1:step[1]:b1+step[4]-step[1]
            for b3 in b2:b2+step[2]-1, b in axes(st, 1)
                if ismatch(b3, ctrl_mask, flag_mask)
                    unsafe_swap_kernel(st, b, b3, step[2], step[4])
                end
            end
        end
    end
    return st
end


text:
using Random

rng = MersenneTwister(123)
n = 100
m = 50
idx = 3

A = randn(rng, n, m)
P = randn(n, n); @. P = (P + P')/2
H = NullspaceHessianLDL(P, A)
# @show eigvals(H.QR.Q2'*P*H.QR.Q2)
@test norm(H.Z[:,end:-1:1]'*P*H.Z[:,end:-1:1] - H.U'*H.D*H.U) <= 1e-7

remove_constraint!(H, 0)
# A = [A[:, 1:idx-1] A[:, idx+1:m]]
@testset "Update qr - remove column" begin
    @test norm(H.QR.Q'*H.QR.Q - I) <= 1e-7
    @test norm(H.QR.Q1*H.QR.R1 - A) <= 1e-7
end

@testset "Expand hessian matrix" begin
    @test norm(H.Z[:,H.m:-1:1]'*P*H.Z[:,H.m:-1:1] - H.U'*H.D*H.U) <= 1e-7
end

a = randn(n)
add_constraint!(H, a)
A = [A a]
@testset "Update qr - add column" begin
    @test norm(H.QR.Q'*H.QR.Q - I) <= 1e-7
    @test norm(H.QR.Q1*H.QR.R1 - A) <= 1e-7
end

@testset "Shrink hessian matrix" begin
    @test norm(H.Z[:,H.m:-1:1]'*P*H.Z[:,H.m:-1:1] - H.U'*H.D*H.U) <= 1e-7
end

# using BenchmarkTools
# @btime qr_hessenberg!($Q_12, $R_2_T')


text:
# SPD-License-Identifier: MIT

using ArgParse
using BayesianModels
using CUDA
using HDF5
using JLD2
using LinearAlgebra
using Random

CUDA.allowscalar(false)

include("utils.jl")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--compression-level", "-c"
            arg_type = Int
            default = 0
            help = "compression level of the HDF5 archive"
        "model"
            required = true
            help = "input GMM"
        "features"
            required = true
            help = "HDF5 archive containing the features"
        "out"
            required = true
            help = "output posteriors in HDF5 format"
    end
    s.description = """
    Compute the GMM's components posterior. The list of utterances
    to use is read from the standard input.
    """
    parse_args(s)
end

function main(args)
    model, _ = loadmodel(args["model"])

    h5open(args["out"], "w") do fout
        h5open(args["features"], "r") do ffea
            for utt in eachline(stdin)
                X = read(ffea[utt])
                fout[utt, compress = args["compression-level"]] = posterior(model, X)
            end
        end
    end
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
# Autogenerated wrapper script for WaveFD_jll for x86_64-apple-darwin
export libillumination, libprop2DAcoIsoDenQ_DEO2_FDTD, libprop2DAcoTTIDenQ_DEO2_FDTD, libprop2DAcoVTIDenQ_DEO2_FDTD, libprop3DAcoIsoDenQ_DEO2_FDTD, libprop3DAcoTTIDenQ_DEO2_FDTD, libprop3DAcoVTIDenQ_DEO2_FDTD, libspacetime

using CompilerSupportLibraries_jll
using FFTW_jll
JLLWrappers.@generate_wrapper_header("WaveFD")
JLLWrappers.@declare_library_product(libillumination, "@rpath/libillumination.dylib")
JLLWrappers.@declare_library_product(libprop2DAcoIsoDenQ_DEO2_FDTD, "@rpath/libprop2DAcoIsoDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libprop2DAcoTTIDenQ_DEO2_FDTD, "@rpath/libprop2DAcoTTIDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libprop2DAcoVTIDenQ_DEO2_FDTD, "@rpath/libprop2DAcoVTIDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libprop3DAcoIsoDenQ_DEO2_FDTD, "@rpath/libprop3DAcoIsoDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libprop3DAcoTTIDenQ_DEO2_FDTD, "@rpath/libprop3DAcoTTIDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libprop3DAcoVTIDenQ_DEO2_FDTD, "@rpath/libprop3DAcoVTIDenQ_DEO2_FDTD.dylib")
JLLWrappers.@declare_library_product(libspacetime, "@rpath/libspacetime.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll)
    JLLWrappers.@init_library_product(
        libillumination,
        "lib/libillumination.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoIsoDenQ_DEO2_FDTD,
        "lib/libprop2DAcoIsoDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoTTIDenQ_DEO2_FDTD,
        "lib/libprop2DAcoTTIDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoVTIDenQ_DEO2_FDTD,
        "lib/libprop2DAcoVTIDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoIsoDenQ_DEO2_FDTD,
        "lib/libprop3DAcoIsoDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoTTIDenQ_DEO2_FDTD,
        "lib/libprop3DAcoTTIDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoVTIDenQ_DEO2_FDTD,
        "lib/libprop3DAcoVTIDenQ_DEO2_FDTD.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspacetime,
        "lib/libspacetime.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MAE103
using Test
using Literate
using JLD # need this to ensure JLD works inside .jl files for Literate use
using Plots # prevents annoying messages inside the generated notebooks

outputdir = "../notebook"
litdir = "./literate"

for (root, dirs, files) in walkdir(litdir)
    if splitpath(root)[end] == "assets"
        for file in files
            cp(joinpath(root, file),joinpath(outputdir,file),force=true)
        end
    end
end

function replace_includes(str)

    included = ["header.jl"]

    # Here the path loads the files from their proper directory,
    # which may not be the directory of the `examples.jl` file!
    path = litdir

    for ex in included
        content = read(joinpath(litdir,ex), String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    return str
end

for (root, dirs, files) in walkdir(litdir)
    for file in files
        #endswith(file,".jl") && startswith(file,"0.0") && Literate.notebook(joinpath(root, file),outputdir,preprocess = replace_includes)
        endswith(file,".jl") && Literate.notebook(joinpath(root, file),outputdir,preprocess = replace_includes)
    end
end


text:
module TestNA

using FactCheck
using DataCubes
using DataCubes: AbstractArrayWrapper, FloatNAArray, simplify_array, wrap_array, naop_plus

facts("NA tests") do
  @fact AbstractArrayWrapper(nalift([1,2,3])) --> AbstractArrayWrapper([Nullable(1), Nullable(2), Nullable(3)])
  @fact map(AbstractArrayWrapper,Any[nalift([1 2 3]),nalift([4,5])]) --> map(AbstractArrayWrapper,Any[[Nullable(1) Nullable(2) Nullable(3)],[Nullable(4), Nullable(5)]])
  @fact AbstractArrayWrapper(simplify_array(nalift(Any[[1 2 3],[4,5]]))) --> AbstractArrayWrapper(simplify_array(Any[Nullable([1 2 3]), Nullable([4,5])]))
  @fact @nalift([1 2 3;4 NA 5]) --> AbstractArrayWrapper([Nullable(1) Nullable(2) Nullable(3);Nullable(4) Nullable{Int}() Nullable(5)])
  @fact @nalift(Any[[1 2 NA],[3,NA]]) --> nalift(map(AbstractArrayWrapper,Any[[Nullable(1) Nullable(2) Nullable{Int}()], [Nullable(3), Nullable{Int}()]]))
  @fact AbstractArrayWrapper([x for x in FloatNAArray([1.0,2.0,NaN])]) --> AbstractArrayWrapper([Nullable(1.0), Nullable(2.0), Nullable{Float64}()])
  @fact AbstractArrayWrapper(FloatNAArray(reverse([1.0,2.0,3.0]))) --> AbstractArrayWrapper(FloatNAArray([3.0,2.0,1.0]))
  @fact simplify_array(nalift([1,2,3])) --> AbstractArrayWrapper([Nullable(1),Nullable(2),Nullable(3)])
  @fact simplify_array(Any[Nullable(3.5), Nullable(2.0), Nullable{Float64}()]) --> AbstractArrayWrapper(FloatNAArray([3.5, 2.0, NaN]))
  @fact (t=@larr(a=[1 2 3;4 5 6], b=[4 5 6;7 8 9]);dcube.setna!(t.data.data[:a],2,2);dcube.setna!(t,2,1);t) --> @larr(a=[1 2 3;NA NA 6], b=[4 5 6;NA 8 9])
  @fact simplify_array(nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
  @fact DataCubes.simplify_floatarray(nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
  @fact sort(nalift([5.0,3.0,4.0])) --> nalift([3.0,4.0,5.0])
  @fact sort!(@nalift([5.0,3.0,4.0])) --> nalift([3.0,4.0,5.0])
  @fact DataCubes.wrap_array(vcat(FloatNAArray(1.0*[1 2 3]),FloatNAArray(1.0*[4 5 6]))) --> DataCubes.wrap_array(FloatNAArray([1 2 3;4 5 6]*1.0))
  @fact DataCubes.wrap_array(hcat(FloatNAArray(1.0*[1 2 3]),FloatNAArray(1.0*[4 5 6]))) --> DataCubes.wrap_array(FloatNAArray([1 2 3 4 5 6]*1.0))
  @fact DataCubes.simplify_floatarray(Nullable{Float64}[]) --> isempty
  @fact typeof(DataCubes.simplify_floatarray(Nullable{Float64}[])) --> DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}
  @fact nalift(FloatNAArray([1.0 2.0 NaN])) --> @nalift([1.0 2.0 NA])
  @fact (t=@nalift([1.0 NA 3]);DataCubes.setna!(t);t) --> nalift(fill(Nullable{Any}(),1,3))
  @fact (t=@nalift([1.0 NA 3]);DataCubes.setna!(t);typeof(t)) --> DataCubes.AbstractArrayWrapper{Nullable,2,Array{Nullable,2}}
  @fact (t=@nalift([1.0 2.0 3.0]);DataCubes.setna!(t);t) --> nalift(fill(NaN,1,3))
  @fact (t=enumeration([:a,:b,:a]);DataCubes.setna!(t);wrap_array(t)) --> wrap_array(fill(Nullable{Symbol}(),3))
  @fact (t=enumeration([:a,:b,:a]);DataCubes.setna!(t,1);wrap_array(t)) --> wrap_array(@enumeration([NA,:b,:a],[:a,:b]))
  @fact (t=enumeration([:a,:b,:a]);DataCubes.setna!(t);typeof(t)) --> EnumerationArray{Symbol,1,DataCubes.AbstractArrayWrapper{Int64,1,Array{Int64,1}},Int64}
  @fact (t=enumeration([:a,:b,:a]);fill!(t,0);wrap_array(t)) --> wrap_array(fill(Nullable{Symbol}(),3))
  @fact (t=enumeration([:a,:b,:a]);fill!(t,Nullable(:c));wrap_array(t)) --> wrap_array(fill(Nullable{Symbol}(),3))
  @fact (t=enumeration([:a,:b,:a]);fill!(t,Nullable{Symbol}());wrap_array(t)) --> wrap_array(fill(Nullable{Symbol}(),3))
  @fact (t=@nalift([1.0 2.0 3.0]);fill!(t,NaN);t) --> @nalift([NA NA NA])
  @fact (t=@nalift([1.0 2.0 3.0]);fill!(t,Nullable{Float64}());t) --> @nalift([NA NA NA])
  @fact (t=@nalift([1.0 2.0 3.0]);fill!(t,Nullable(5.0));t) --> @nalift([5.0 5.0 5.0])
  @fact igna(Nullable{Int}(), 3) --> 3
  @fact igna(Nullable(1), 3) --> 1
  @fact igna(Nullable(1.0), 3.0) --> 1.0
  @fact igna(Array(Nullable{Int},0)) --> Array(Int,0)
  @fact igna(Array(Nullable{Int},0),1) --> Array(Int,0)
  @fact_throws igna(@nalift([1,2,NA]))
  @fact_throws igna(@nalift([1,2.0,NA]))
  @fact igna(@nalift([1 3 NA]), 1) --> [1 3 1]
  @fact igna(@nalift([1 3.0 NA]), 1) --> [1 3.0 1]
  @fact igna(Array(Nullable,0), 1) --> Array(Nullable,0)
  @fact igna(@nalift([1 3 5]), 1) --> [1 3 5]
  @fact igna(@nalift([1.0 3.0 NA]), 1.0) --> [1.0 3.0 1.0]
  @fact igna(@nalift([1.0 3.0 NA]))[3] --> isnan
  @fact igna(@nalift([1.0 3.0 5.0]), 1.0) --> [1.0 3.0 5.0]
  @fact igna([Nullable(1.0) Nullable(3.0) Nullable{Float64}()])[3] --> isnan
  @fact igna([Nullable(1.0) Nullable(3.0) Nullable(5.0)], 1.0) --> [1.0 3.0 5.0]
  @fact igna([Nullable(1.0) Nullable(3.0) Nullable{Float64}()], 1.0) --> [1.0 3.0 1.0]
  @fact isna(Nullable{Int}()) --> true
  @fact isna(Nullable(1)) --> false
  @fact isna(Nullable(1.0)) --> false
  @fact isna(@nalift([1 3 NA])) --> [false false true]
  @fact isna(@nalift([1 3 NA]), 2:3) --> [false false true][2:3]
  @fact isna(@nalift([1 3 5])) --> [false false false]
  @fact isna(@nalift([1.0 3.0 NA])) --> [false false true]
  @fact isna(@nalift([1.0 3.0 5.0])) --> [false false false]
  @fact isna(@nalift([1.0 3.0 5 NA])) --> [false false false true]
  @fact isna(@enumeration([:a :b NA;:d NA :f])) --> [false false true;false true false]
  @fact isna(@enumeration([:a :b NA;:d NA :f]), 1:2, 3) --> [false false true;false true false][1:2, 3]
  @fact isna([Nullable{Float64}() Nullable(1.0)]) --> [true false]
  @fact ignabool(nalift([true, false])) --> [true, false]
  @fact ignabool(Nullable{Bool}()) --> false
  @fact ignabool(Nullable(true)) --> true
  @fact ignabool(Nullable(false)) --> false
  @fact ignabool([Nullable(true), Nullable{Bool}()]) --> [true, false]
  @fact ignabool([Nullable(true)]) --> [true]
  @fact ignabool([Nullable(false)]) --> [false]
  @fact ignabool(true) --> true
  @fact ignabool(false) --> false
  @fact @nalift([NA NA]) --> wrap_array([Nullable{Any}() Nullable{Any}()])
  @fact @nalift([NA,NA]) --> wrap_array([Nullable{Any}(),Nullable{Any}()])
  @fact @nalift([NA NA;NA NA]) --> wrap_array(reshape(fill(Nullable{Any}(),4),2,2))
  @fact @nalift([1,NA]) --> wrap_array([Nullable(1), Nullable{Int}()])
  @fact @nalift([1 NA]) --> wrap_array([Nullable(1)  Nullable{Int}()])
  @fact @nalift([1 NA;NA NA]) --> wrap_array([Nullable(1)  Nullable{Int}();Nullable{Int}() Nullable{Int}()])
  @fact map(x->x, nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
  @fact map((x,y)->naop_plus(x, y), nalift([1.0,2.0,3.0]), nalift([2.0,3.0,4.0])) --> nalift([3.0,5.0,7.0])
  #@fact nalift(DataCubes.simplify_array(Any[larr(a=[1,2,3]), darr(b=[:a,:b,:c])])) --> DataCubes.simplify_array(Any[larr(a=[1,2,3]), darr(b=[:a,:b,:c])])
  #@fact nalift(DataCubes.simplify_array(Any[darr(a=[1,2,3]), darr(b=[:a,:b,:c])])) --> DataCubes.simplify_array(Any[darr(a=[1,2,3]), darr(b=[:a,:b,:c])])
  #@fact nalift(DataCubes.simplify_array(Any[larr(a=[1,2,3]), larr(b=[:a,:b,:c])])) --> DataCubes.simplify_array(Any[larr(a=[1,2,3]), larr(b=[:a,:b,:c])])
  @fact nalift(nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
  @fact pickaxis(larr(a=[1,2,3]),1) == nalift([1,2,3]) --> false
  @fact nalift([1,2,3]) == pickaxis(larr(a=[1,2,3]),1) --> false

  context("FloatNAArray tests") do
    @fact map(identity, FloatNAArray(Array(Float64,0))) --> Array(Nullable{Any},0)
    @fact nalift([1.0,2.0]) --> DataCubes.wrap_array(FloatNAArray([1.0,2.0]))
    @fact FloatNAArray([1.0,2.0,3.0])[1].value --> 1.0
    @fact FloatNAArray([1.0,2.0,3.0])[1:2].data --> [1.0,2.0]
    @fact wrap_array(DataCubes.simplify_floatarray(FloatNAArray([1.0 2.0]))) --> wrap_array(FloatNAArray([1.0 2.0]))
    @fact typeof(DataCubes.simplify_floatarray(FloatNAArray(Array(Float64,0)))) --> FloatNAArray{Float64,1,Array{Float64,1}}
    @fact nalift(nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
    @fact typeof(igna(nalift(Array(Float64,0)))) --> DataCubes.AbstractArrayWrapper{Float64,1,Array{Float64,1}}
    @fact typeof(igna(nalift(Array(Float64,0)),3.0)) --> DataCubes.AbstractArrayWrapper{Float64,1,Array{Float64,1}}
    @fact (arr=nalift([1.0,2.0,3.0]);arr[1]=0.0;arr) --> nalift([0.0,2.0,3.0])
    @fact (arr=nalift([1.0,2.0,3.0]);arr[1]=Nullable{Float64}();arr) --> @nalift([NA,2.0,3.0])
    @fact (arr=nalift([1.0,2.0,3.0]);arr[1:2]=nalift([0.0,-1.0]) ;arr) --> nalift([0.0,-1.0,3.0])
    @fact (arr=nalift([1.0,2.0,3.0]);arr[1:2]=Nullable(1.0) ;arr) --> nalift([1.0,1.0,3.0])
    @fact typeof(1.0*nalift([1 2])) --> AbstractArrayWrapper{Nullable{Float64},2,DataCubes.FloatNAArray{Float64,2,Array{Float64,2}}}
    @fact sub(1.0*nalift([1 2 3;4 5 6]),1:2,1:1) --> nalift([1.0 4.0]')
    @fact slice(1.0*nalift([1 2 3;4 5 6]),1:2,1) --> nalift([1.0,4.0])
    @fact sub(1.0*nalift([1 2 3;4 5 6]),(1:2,1:1)) --> nalift([1.0 4.0]')
    @fact slice(1.0*nalift([1 2 3;4 5 6]),(1:2,1)) --> nalift([1.0,4.0])
  end
end

end


text:
"""
    weighted_iterator(d)

Return an iterator through pairs of the values and probabilities in distribution `d`.

This is designed to speed up value iteration. Distributions are encouraged to provide a custom optimized implementation if possible.

# Example
```julia-repl
julia> d = BoolDistribution(0.7)
BoolDistribution(0.7)

julia> collect(weighted_iterator(d))
2-element Array{Pair{Bool,Float64},1}:
  true => 0.7
 false => 0.3
```
"""
weighted_iterator(d) = (x=>pdf(d, x) for x in support(d))


text:
using MeshArrays, OceanStateEstimation, NetCDF

p=dirname(pathof(MeshArrays))
include(joinpath(p,"../examples/Demos.jl"))

"""
    random_flow_field(;np=12,nq=18)

Set up a random flow field over a gridded domain of size np,nq

```
ϕ,u,v=random_flow_field()
```
"""
function random_flow_field(;np=12,nq=18)

Γ=simple_periodic_domain(np,nq)
(_,ϕ,_,_)=demo2(Γ)
ϕ .*= 0.5

#For the convergent / scalar potential case, ϕ is interpreted as being
#on center points -- hence the standard gradient function readily gives
#what we need
#(u,v)=gradient(ϕ,Γ)
#return u[1],v[1],ϕ[1]

#For the rotational / streamfunction case, ϕ is interpreted as being
#on S/W corner points -- this is ok since the grid is homegeneous,
#and conveniently yields an adequate substitution u,v <- -v,u; but note
#that doing the same with gradient() would shift indices inconsistenly
u=-(circshift(ϕ[1], (0,-1))-ϕ[1])
v=(circshift(ϕ[1], (-1,0))-ϕ[1])
return u,v,ϕ[1]
end

"""
    solid_body_rotation(np,nz)

Set up an idealized flow field which consists of
[rigid body rotation](https://en.wikipedia.org/wiki/Rigid_body),
plus a convergent term, plus a sinking term.

```
u,v,w=solid_body_rotation(12,4)
```
"""
function solid_body_rotation(np,nz)
    Γ=simple_periodic_domain(np);
    γ=Γ.XC.grid;

    #Solid-body rotation around central location ...
    i=Int(np/2+1)
    u=-(Γ.YG.-Γ.YG[1][i,i])
    v=(Γ.XG.-Γ.XG[1][i,i])

    #... plus a convergent term to / from central location
    d=-0.01
    u=u+d*(Γ.XG.-Γ.XG[1][i,i])
    v=v+d*(Γ.YG.-Γ.YG[1][i,i])

    #Replicate u,v in vertical dimension
    uu=MeshArray(γ,γ.ioPrec,nz)
    [uu[k]=u[1] for k=1:nz]
    vv=MeshArray(γ,γ.ioPrec,nz)
    [vv[k]=v[1] for k=1:nz]

    #Vertical velocity component w
    w=fill(-0.01,MeshArray(γ,γ.ioPrec,nz+1))

    return write(uu),write(vv),write(w)
end

"""
    global_ocean_circulation(;k=10,ny=2)

Set up Global Ocean particle simulation in 2D with seasonally varying flow field.

```
𝑃,𝐷=global_ocean_circulation(k=10,ny=2);
```
"""
function global_ocean_circulation(;k=1,ny=2)

  #k=10 #choice of vertical level
  #ny=2 #number of simulated years (20 for k>20)
  r_reset = 0.01 #fraction of the particles reset per month (0.05 for k<=10)

  #read grid and set up connections between subdomains
  p=dirname(pathof(IndividualDisplacements))
  γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
  Γ=GridLoad(γ)
  Γ=merge(Γ,NeighborTileIndices_cs(Γ))

  func=(u -> update_location_llc!(u,𝐷))
  Γ=merge(Γ,(; update_location! = func))

  #initialize u0,u1 etc
  𝑃,𝐷=set_up_FlowFields(k,Γ,ECCOclim_path);

  #add parameters for use in reset!
  tmp=(frac=r_reset, Γ=Γ)
  𝐷=merge(𝐷,tmp)

  return 𝑃,𝐷

end

"""
    OCCA_FlowFields(;backward_in_time::Bool=false,nmax=Inf)

Define gridded variables and return result as NamedTuple
"""
function OCCA_FlowFields(;backward_in_time::Bool=false,nmax=Inf)

   γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
   Γ=GridLoad(γ)
   n=length(Γ.RC)
   isfinite(nmax) ? n=min(n,Int(nmax)) : nothing

   g=Γ.XC.grid
   func=(u -> IndividualDisplacements.update_location_dpdo!(u,g))

   jj=[:hFacC, :hFacW, :hFacS, :DXG, :DYG, :RAC, :RAZ, :RAS]
   ii=findall([!in(i,jj) for i in keys(Γ)])
   Γ=(; zip(Symbol.(keys(Γ)[ii]), values(Γ)[ii])...)

   backward_in_time ? s=-1.0 : s=1.0
   s=Float32(s)

   function rd(filename, varname,n)
   fil = NetCDF.open(filename, varname)
   siz = size(fil)
   tmp = zeros(siz[1:2]...,n)
   [tmp .+= fil[:,:,1:n,t] for t=1:12]
   tmp ./= 12.0
   tmp[findall(tmp.<-1e22)] .= 0.0
   return tmp
   end

   fileIn=OCCAclim_path*"DDuvel.0406clim.nc"
   u=s*read(rd(fileIn,"u",n),MeshArray(γ,Float32,n))

   fileIn=OCCAclim_path*"DDvvel.0406clim.nc"
   v=s*read(rd(fileIn,"v",n),MeshArray(γ,Float32,n))

   fileIn=OCCAclim_path*"DDwvel.0406clim.nc"
   w=s*rd(fileIn,"w",n)
   w=-cat(w,zeros(360, 160),dims=3)
   w[:,:,1] .=0.0
   w=read(w,MeshArray(γ,Float32,n+1))

   fileIn=OCCAclim_path*"DDtheta.0406clim.nc"
   θ=read(rd(fileIn,"theta",n),MeshArray(γ,Float32,n))

#   fileIn=OCCAclim_path*"DDsalt.0406clim.nc"
#   𝑆=read(rd(fileIn,"salt",n),MeshArray(γ,Float64,n))

   for i in eachindex(u)
      u[i]=u[i]./Γ.DXC[1]
      v[i]=v[i]./Γ.DYC[1]
   end

   for i in eachindex(u)
      u[i]=circshift(u[i],[-180 0])
      v[i]=circshift(v[i],[-180 0])
      θ[i]=circshift(θ[i],[-180 0])
#      𝑆[i]=circshift(𝑆[i],[-180 0])
   end

   for i in eachindex(w)
      w[i]=w[i]./Γ.DRC[min(i[2]+1,n)]
      w[i]=circshift(w[i],[-180 0])
   end

   tmpx=circshift(Γ.XC[1],[-180 0])
   tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
   Γ.XC[1]=tmpx

   tmpx=circshift(Γ.XG[1],[-180 0])
   tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
   Γ.XG[1]=tmpx
   Γ.Depth[1]=circshift(Γ.Depth[1],[-180 0])

   t0=0.0; t1=86400*366*2.0;

   for k=1:n
    (tmpu,tmpv)=exchange(u[:,k],v[:,k],1)
    u[:,k]=tmpu
    v[:,k]=tmpv
   end
   for k=1:n+1
    tmpw=exchange(w[:,k],1)
    w[:,k]=tmpw
   end

   𝑃=FlowFields(u,u,v,v,w,w,[t0,t1],func)

   𝐷 = (θ0=θ, θ1=θ, XC=exchange(Γ.XC), YC=exchange(Γ.YC),
   RF=Γ.RF, RC=Γ.RC,ioSize=(360,160,n))

   return 𝑃,𝐷,Γ

end

"""
    test1_setup()

Call `gcmgrid`, initialize a single point,
rely on `dxdt!`, and just output `sol` at the end.

```
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
𝑃,sol=test1_setup()
```
"""
function test1_setup()

    mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
    XC=MeshArray(mygrid,Float32); XC[1]=vec(2500.:5000.:397500.0)*ones(1,42);
    XG=MeshArray(mygrid,Float32); XG[1]=vec(0.:5000.:395000.0)*ones(1,42);
    YC=MeshArray(mygrid,Float32); YC[1]=ones(80,1)*transpose(vec(2500.:5000.:207500.0));
    YG=MeshArray(mygrid,Float32); YG[1]=ones(80,1)*transpose(vec(0.:5000.:205000.0));

    dx=5000.0
    t0=0.0; t1=18001.0*3600.0
    u=-(YG.-YC[1][40,21])/2000000.
    v=(XG.-XC[1][40,21])/2000000.
    u0=u./dx; u1=u./dx
    v0=v./dx; v1=v./dx

    𝑃=FlowFields(u0[1], u1[1], v0[1], v1[1], [t0,t1])

    u0=[200000.0;0.0]./dx
    du=fill(0.0,2);
    prob = ODEProblem(dxdt!,u0,[0.0,2998*3600.0],𝑃)
    sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

    return 𝑃,sol
end

"""
    test2_periodic_domain(np = 12, nq = 12)

Call `simple_periodic_domain`, initialize 6x6 point cloud,
rely on `dxdt!`, and call `postprocess_xy` at the end.

```
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
df,𝑃=test2_periodic_domain()

using Plots
@gif for t in 𝑃.t0:1.0:𝑃.t1
   scatter_subset(𝑃,df,t)
end
```
"""
function test2_periodic_domain(np = 12, nq = 12)
    #domain and time parameters
    Γ = simple_periodic_domain(np, nq)

    u = 0.1 ./ Γ.DXC
    v = 0.3 ./ Γ.DYC
    (u, v) = exchange(u, v, 1)

    f = (u -> IndividualDisplacements.update_location_dpdo!(u,Γ.XC.grid))
    𝑃=FlowFields(u,u,v,v,[0.0,400.0],f)

    #initial conditions
    x0 = np * (0.4:0.04:0.6)
    y0 = nq * (0.4:0.04:0.6)
    x0 = vec(x0) * ones(1, length(y0))
    y0 = ones(size(x0, 1), 1) * transpose(vec(y0))
    u0 = permutedims([[x0[i];y0[i];1.0] for i in eachindex(x0)])
    du=0*u0

    #solve for trajectories
    prob = ODEProblem(dxdt!, u0, 𝑃.𝑇, 𝑃)
    sol = solve(prob,Euler(),dt=0.1)

    return postprocess_xy(sol, 𝑃),𝑃
end


text:
using HashDicts: HashDict
using Test

import HashDicts

@testset "HashDict construction" begin
    d = HashDict([("A", 1), ("B", 2)])
    d2 = HashDict("A" => 1, "B" => 2)
    @test d["A"] == d2["A"] == 1
    @test d["B"] == d2["B"] == 2
end

@testset "Bucket construction and iteration" begin
    b = HashDicts.Bucket([Pair(1, 2), Pair(2, 8)])
    for (k, v) in b
        @test k == 1 || k == 2
        @test v == 2 || v == 8
    end
    b[8] = 33
    @test b[8] == 33
end

@testset "Number of entries kept up-to-date" begin
    d = HashDict([("A", 1), ("B", 2)])
    @test d.number_entries == 2
    @test length(d) == 2
    d["T"] = 2
    @test d.number_entries == 3
    @test length(d) == 3
    d["T"] = 1
    @test d.number_entries == 3
    d2 = delete!(d, "T")
    @test d2 === d
    @test d.number_entries == 2
    @test_throws KeyError d["T"]
end

@testset "Iteration over HashDict" begin
    d = HashDict{Int,Int}()
    for value in d
    end
    d[1] = 42
    for (k, v) in d
        @test k == 1
        @test v == 42
    end
    d[2] = 42
    for (k, v) in d
        @test k < 3
        @test v == 42
    end
    d2 = Dict{Int, UInt8}()
    idxs = rand(Int, 100)
    vals = rand(UInt8, 100)
    for i in eachindex(idxs)
        d2[idxs[i]] = vals[i]
        @test haskey(d2, idxs[i])
    end
    for (k, v) in d2
        idx = findfirst(==(k), idxs)
        @test idx !== nothing
        d2[idxs[idx]] = vals[idx]
    end
    d3 = Dict{Int, Int}()
    idxs = rand(Int, 100)
    vals = abs.(rand(Int, 100))
    for i in eachindex(idxs)
        d3[idxs[i]] = vals[i]
    end
    for (k, v) in d3
        idx = findfirst(==(k), idxs)
        @test idx !== nothing
        d3[idxs[idx]] = vals[idx]
        @test get(d3, idxs[idx], -1) != -1
    end

    # deleting and default
    for key in idxs[1:10]
        delete!(d3, key)
        @test get(d3, key, -1) == -1
    end
end

@testset "Emptying HashDict" begin
    d = HashDict([("A", 1), ("B", 2)])
    @test length(d) == 2
    empty!(d)
    @test length(d) == 0
    for (k, v) in d
        @test v > 5
    end
end

@testset "Displaying HashDicts" begin
    d = HashDict(2 => "", 3 => "hi")
    io = IOBuffer()
    show(io, d)
    result_str = String(take!(io))
    @test result_str == """
HashDict{Int64,String} with 2 entries:
  2 => ""
  3 => "hi"
"""
end

@testset "Resizing does not affect HashDict" begin
    d = HashDict(2 => "", 3 => "hi")
    initial_size = length(d.buckets)
    @test length(d) == 2
    for (k, v) in d
        @test 2 <= k <= 3
        @test 0 <= length(v) <= 2
    end
    resize!(d)
    @test length(d) == 2
    for (k, v) in d
        @test 2 <= k <= 3
        @test 0 <= length(v) <= 2
    end
    @test length(d.buckets) == 2 * initial_size
end

@testset "Adding elements causes resize" begin
    d = HashDict{Int, String}(base_size = 10)
    @test length(d.buckets) == 10
    for i in 1:10:301
        d[i] = string(i)
    end
    @test length(d.buckets) == 20
    for i in 1:60
        d[i+301] = string(i)
    end
    @test length(d.buckets) == 40
end

@testset "Allocating capacity avoids automatic resize!" begin
    d = HashDict{Int, String}(base_size = 10)
    @test length(d.buckets) == 10
    for i in 1:10:301
        d[i] = string(i)
    end
    @test length(d.buckets) == 20
    resize!(d, 15)
    @test length(d.buckets) == 35
    for i in 1:60
        d[i+301] = string(i)
    end
    @test length(d.buckets) == 35
end

@testset "Collect eltype is Pair{K, V}" begin
    d = HashDict(2 => "", 3 => "hi")
    for p in d
        p isa Pair
    end
    collected = collect(d)
    sort!(collected, by=(p -> p.first))
    @test collected[1] == (2 => "")
    @test collected[2] == (3 => "hi")
end


text:
#######################################
# Iterator structures

struct DTableRowIterator
    d::DTable
end

struct DTableColumnIterator
    d::DTable
end

struct DTablePartitionIterator
    d::DTable
end

#######################################
# DTable functions

Tables.istable(table::DTable) = true
Tables.rowaccess(table::DTable) = true
Tables.rows(table::DTable) = DTableRowIterator(table)
Tables.columnaccess(table::DTable) = true
Tables.columns(table::DTable) = DTableColumnIterator(table)

function determine_schema(table::DTable)
    if table.schema !== nothing
        return table.schema
    end
    # Figure out schema
    chunk_f = chunk -> begin
        r = isnonempty(chunk)
        (r, r ? Tables.schema(Tables.rows(chunk)) : nothing)
    end
    c_idx = 1
    r = (false, nothing)
    while !r[1]
        r = fetch(Dagger.spawn(chunk_f, table.chunks[c_idx]))
        c_idx += 1
    end
    # cache results
    table.schema = r[2]
end

function determine_columnnames(table::DTable)
    s = determine_schema(table)
    s === nothing ? nothing : s.names
end

function _getcolumn(table::DTable, col::Union{Symbol, Int})
    chunk_col = (_chunk, _col) -> Tables.getcolumn(_chunk, _col)
    v = [Dagger.spawn(chunk_col, chunk, col) for chunk in table.chunks]
    SentinelArrays.ChainedVector(fetch.(v))
end

Tables.getcolumn(table::DTable, col::Symbol) = _getcolumn(table, col)
Tables.getcolumn(table::DTable, idx::Int) = _getcolumn(table, idx)

#######################################
# DTableRowIterator functions

Tables.schema(table::DTableRowIterator) = determine_schema(table.d)
length(table::DTableRowIterator) = length(table.d)

function _iterate(iter::DTableRowIterator, chunk_index)
    i = nothing
    row_iterator = nothing
    while i === nothing && chunk_index <= nchunks(iter.d)
        partition = _retrieve(iter.d.chunks[chunk_index])
        row_iterator = Tables.rows(partition)
        i = iterate(row_iterator)
        chunk_index += 1
    end
    if i === nothing
        return nothing
    else
        row, row_state = i
        next_chunk_index = chunk_index
        return (row, (row_iterator, row_state, next_chunk_index))
    end
end

Base.iterate(iter::DTableRowIterator) = _iterate(iter, 1)

function Base.iterate(iter::DTableRowIterator, state)
    (row_iterator, row_state, next_chunk_index) = state
    i = iterate(row_iterator, row_state)
    if i === nothing
        _iterate(iter, next_chunk_index)
    else
        row, row_state = i
        return (row, (row_iterator, row_state, next_chunk_index))
    end
end

#######################################
# DTableColumnIterator functions

Tables.schema(table::DTableColumnIterator) = determine_schema(table.d)
Tables.columnnames(table::DTableColumnIterator) = determine_columnnames(table.d)
Tables.getcolumn(table::DTableColumnIterator, col::Symbol) = Tables.getcolumn(table.d, col)
Tables.getcolumn(table::DTableColumnIterator, idx::Int) = Tables.getcolumn(table.d, idx)
length(table::DTableColumnIterator) = length(Tables.columnnames(table))

function _iterate(table::DTableColumnIterator, column_index)
    columns = Tables.columnnames(table)
    if (columns === nothing || length(columns) < column_index)
        return nothing
    else
        return (Tables.getcolumn(table, column_index), column_index + 1)
    end
end

Base.iterate(table::DTableColumnIterator) = _iterate(table, 1)
Base.iterate(table::DTableColumnIterator, state) = _iterate(table, state)

#######################################
# DTablePartitionIterator functions

Tables.partitions(table::DTable) = DTablePartitionIterator(table)
length(table::DTablePartitionIterator) = nchunks(table.d)

function _iterate(table::DTablePartitionIterator, chunk_index)
    nchunks(table.d) < chunk_index && return nothing
    (_retrieve(table.d.chunks[chunk_index]), chunk_index + 1)
end

Base.iterate(table::DTablePartitionIterator) = _iterate(table, 1)
Base.iterate(table::DTablePartitionIterator, state) = _iterate(table, state)

#######################################
# GDTable
# For normal rows/columns access it should act the same as a DTable

Tables.istable(table::GDTable) = true
Tables.rowaccess(table::GDTable) = true
Tables.rows(table::GDTable) = DTableRowIterator(table.dtable)
Tables.columnaccess(table::GDTable) = true
Tables.columns(table::GDTable) = DTableColumnIterator(table.dtable)
Tables.schema(table::GDTable) = determine_schema(table.dtable)
Tables.getcolumn(table::GDTable, col::Symbol) = Tables.getcolumn(table.dtable, col)
Tables.getcolumn(table::GDTable, idx::Int) = Tables.getcolumn(table.dtable, idx)
Tables.columnnames(table::GDTable) = determine_columnnames(table.dtable)

#######################################
# GDTable partitions
# Here it makes sense to provide partitions as full key groups
# Same as normal iteration over GDTable, but returns partitions only without keys

struct GDTablePartitionIterator
    d::GDTable
end

Tables.partitions(table::GDTable) = GDTablePartitionIterator(table)

function _iterate(table::GDTablePartitionIterator, it)
    if it === nothing
        return nothing
    else
        ((_, partition), index_iter_state) = it
        return (partition, index_iter_state)
    end
end

iterate(table::GDTablePartitionIterator) = _iterate(table, iterate(table.d))
iterate(table::GDTablePartitionIterator, index_iter_state) = _iterate(table, iterate(table.d, index_iter_state))


text:
using UniqueVectors: UniqueVector
using LinearAlgebra: issymmetric, diag, I
using SparseArrays: spzeros

SparseMat = SparseMatrixCSC{Float64,Int}

export Circuit
struct Circuit{T}
    k::SparseMat # inverse inductances
    g::SparseMat # conductances
    c::SparseMat # capacitances
    vertices::UniqueVector{T}
    function Circuit{T}(k::SparseMat, g::SparseMat, c::SparseMat,
        vertices::UniqueVector{T}) where T
        for m in [k, g, c]
            @assert issymmetric(m)
            @assert size(m) == (length(vertices), length(vertices))
            @assert all(diag(m) .== 0)
            @assert all(m .>= 0)
        end
        return new{T}(k, g, c, vertices)
    end
end

Circuit(k::SparseMat, g::SparseMat, c::SparseMat, vertices::AbstractVector{T}) where T =
    Circuit{T}(k, g, c, UniqueVector(vertices))

Circuit(k::AbstractMatrix{<:Real}, g::AbstractMatrix{<:Real}, c::AbstractMatrix{<:Real},
    vertices::AbstractVector) = Circuit(sparse(1.0 * k), sparse(1.0 * g), sparse(1.0 * c),
    vertices)

function Circuit(vertices::AbstractVector{T}) where T
    z() = spzeros(length(vertices), length(vertices))
    return Circuit(z(), z(), z(), vertices)
end

export partial_copy
function partial_copy(circ::Circuit{T};
   k::Union{SparseMat, Nothing}=nothing,
   g::Union{SparseMat, Nothing}=nothing,
   c::Union{SparseMat, Nothing}=nothing,
   vertices::Union{AbstractVector{U}, Nothing}=nothing) where {T, U}
   k = k == nothing ? circ.k : k
   g = g == nothing ? circ.g : g
   c = c == nothing ? circ.c : c
   vertices = vertices == nothing ? circ.vertices : vertices
   return Circuit(k, g, c, vertices)
end

import Base: ==
function ==(circ1::Circuit, circ2::Circuit)
    if typeof(circ1) != typeof(circ2)
        return false
    end
    return all([getfield(circ1, name) == getfield(circ2, name)
        for name in fieldnames(Circuit)])
end

import Base: isapprox
function isapprox(circ1::Circuit, circ2::Circuit)
    if typeof(circ1) != typeof(circ2)
        return false
    end
    return all([getfield(circ1, name) ≈ getfield(circ2, name)
        for name in fieldnames(Circuit)[1:end-1]])
end

export matrices
matrices(c::Circuit) = [c.k, c.g, c.c]

function get_matrix_element(c::Circuit{T}, matrix_name::Symbol, v0::T, v1::T) where T
    i0 = findfirst(isequal(v0), c.vertices)
    i1 = findfirst(isequal(v1), c.vertices)
    return getfield(c, matrix_name)[i0, i1]
end

export get_inv_inductance
get_inv_inductance(c::Circuit{T}, v0::T, v1::T) where T = get_matrix_element(c, :k, v0, v1)

export get_inductance
get_inductance(c::Circuit{T}, v0::T, v1::T) where T = 1/get_inv_inductance(c, v0, v1)

export get_conductance
get_conductance(c::Circuit{T}, v0::T, v1::T) where T = get_matrix_element(c, :g, v0, v1)

export get_resistance
get_resistance(c::Circuit{T}, v0::T, v1::T) where T = 1/get_conductance(c, v0, v1)

export get_capacitance
get_capacitance(c::Circuit{T}, v0::T, v1::T) where T = get_matrix_element(c, :c, v0, v1)

export get_elastance
get_elastance(c::Circuit{T}, v0::T, v1::T) where T = 1/get_capacitance(c, v0, v1)

function set_matrix_element!(c::Circuit{T}, matrix_name::Symbol, v0::T, v1::T,
    value::Real) where T
    @assert value >= 0
    @assert value < Inf
    i0 = findfirst(isequal(v0), c.vertices)
    i1 = findfirst(isequal(v1), c.vertices)
    @assert i0 != i1
    matrix = getfield(c, matrix_name)
    matrix[i0, i1] = value
    matrix[i1, i0] = value
    return c
end

export set_inv_inductance!
set_inv_inductance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_matrix_element!(c, :k, v0, v1, value)

export set_inductance!
set_inductance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_inv_inductance!(c, v0, v1, 1/value)

export set_conductance!
set_conductance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_matrix_element!(c, :g, v0, v1, value)

export set_resistance!
set_resistance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_conductance!(c, v0, v1, 1/value)

export set_capacitance!
set_capacitance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_matrix_element!(c, :c, v0, v1, value)

export set_elastance!
set_elastance!(c::Circuit{T}, v0::T, v1::T, value::Real) where T =
    set_capacitance!(c, v0, v1, 1/value)

#######################################################
# spanning trees
#######################################################

export SpanningTree
struct SpanningTree{T}
    root::T
    edges::Vector{Tuple{T, T}}
end

# depth-1 spanning tree
function SpanningTree(vertices::AbstractVector{T}) where T
    @assert length(vertices) > 0
    root = vertices[1]
    edges = [(v, root) for v in vertices[2:end]]
    return SpanningTree{T}(root, edges)
end

export coordinate_matrix
# T matrix in section IIC with a given spanning tree.
# The direction of the edges passed in is ignored and enforced as being outwards
# from the root
function coordinate_matrix(c::Circuit{T}, tree::SpanningTree{T}) where T
    V = length(c.vertices)
    @assert tree.root in c.vertices
    @assert length(tree.edges) == V-1 # the correct number of edges for a spanning tree
    mapping = Dict(tree.root => spzeros(V-1)) # the function l in equation 2
    edge_set = Set(enumerate(tree.edges))
    while !isempty(edge_set)
        for e in edge_set
            (index, (v0, v1)) = e
            @assert !(v0 in keys(mapping) && v1 in keys(mapping)) # edges contains a loop
            if v0 in keys(mapping) || v1 in keys(mapping)
                delete!(edge_set, e)
                if v0 in keys(mapping) # enforces all edges are oriented away from the root
                    mapping[v1] = 1.0 * mapping[v0]
                    mapping[v1][index] += 1
                else
                    mapping[v0] = 1.0 * mapping[v1]
                    mapping[v0][index] += 1
                end
            end
        end
    end
    @assert length(mapping) == V # all vertices have been assigned a coordinate
    return hcat([mapping[v] for v in c.vertices]...)
end

# T matrix in section IIC for depth-1 spanning tree
coordinate_matrix(num_vertices::Int) = transpose([spzeros(1, num_vertices-1); I])
coordinate_matrix(c::Circuit) = coordinate_matrix(length(c.vertices))

#######################################################
# cascade and unite
#######################################################

export cascade
function cascade(circs::AbstractVector{Circuit{T}}) where T
    circuit_matrices = [cat(m..., dims=(1,2))
        for m in zip([matrices(circ) for circ in circs]...)]
    vertices = vcat([circ.vertices for circ in circs]...)
    return Circuit(circuit_matrices..., vertices)
end

cascade(circs::Vararg{Circuit{T}}) where T = cascade(collect(circs))

export unite_vertices
function unite_vertices(circ::Circuit{T}, vertices::AbstractVector{T}) where T
    if length(vertices) <= 1
        return circ
    end
    n = length(circ.vertices)
    vertex_inds = [findfirst(isequal(v), circ.vertices) for v in vertices]
    keep_inds = filter(x -> !(x in vertex_inds[2:end]), 1:n)
    function unite_indices(mat::AbstractMatrix)
        m = copy(mat)
        ind1 = vertex_inds[1]
        for ind2 in vertex_inds[2:end]
            m[ind1, :] += mat[ind2, :]
            m[:, ind1] += mat[:, ind2]
        end
        m = m[keep_inds, keep_inds]
        for i in 1:size(m, 1)
            m[i, i] = 0
        end
        return m
    end
    circuit_matrices = [unite_indices(mat) for mat in matrices(circ)]
    return Circuit(circuit_matrices..., circ.vertices[keep_inds])
end

unite_vertices(circ::Circuit{T}, vertices::Vararg{T}) where T =
    unite_vertices(circ, collect(vertices))

export cascade_and_unite
# cascade all circuits and unite vertices with the same name
function cascade_and_unite(circs::AbstractVector{Circuit{T}}) where T
    @assert length(circs) >= 1
    if length(circs) == 1
        return circs[1]
    end
    # number the vertices so that the names are all distinct and then cascade
    vertex_number = 1
    function rename(circ)
        vertices = [(vertex_number + i - 1, vertex)
            for (i, vertex) in enumerate(circ.vertices)]
        vertex_number += length(vertices)
        return Circuit(matrices(circ)..., vertices)
    end
    circ = cascade(map(rename, circs))
    # merge all vertices with the same name
    original_vertices = vcat([c.vertices for c in circs]...)
    for vertex in unique(original_vertices)
        inds = findall([v[2] == vertex for v in circ.vertices])
        circ = unite_vertices(circ, circ.vertices[inds])
    end
    # remove numbering
    return Circuit(matrices(circ)..., [v[2] for v in circ.vertices])
end

cascade_and_unite(circs::Vararg{Circuit{T}}) where T = cascade_and_unite(collect(circs))


text:
function _trajectory(e::AbstractEnvironment, T::Integer)
    (
        states = Array(undef, statespace(e), T),
        obses = Array(undef, observationspace(e), T),
        acts = Array(undef, actionspace(e), T),
        rews = Array(undef, rewardspace(e), T),
    )
end

function _rollout(e::AbstractEnvironment, actions::AbstractMatrix)
    T = size(actions, 2)
    traj = _trajectory(e, T)
    traj.acts .= actions
    for t=1:T
        st = view(traj.states, :, t)
        at = view(traj.acts, :, t)
        ot = view(traj.obses, :, t)

        getstate!(st, e)
        getobservation!(ot, e)
        setaction!(e, at)
        traj.rews[t] = getreward(st, at, ot, e)
        step!(e)
    end
    traj
end

function testenv_correctness(etype::Type{<:AbstractEnvironment}, args...; kwargs...)
    makeenv() = etype(args...; kwargs...)
    randstate() = (e = makeenv(); randreset!(e); getstate(e))
    randobs() = (e = makeenv(); randreset!(e); getobservation(e))
    function randaction()
       e = makeenv()
       a = rand(actionspace(e)) .- 0.5 # random controllers in range [-0.5, 0.5]
       setaction!(e, a)
       # TODO hack for case where env has action limits until BoundedShape exists
       getaction!(a, e)
       return a
    end

    @testset "Interface" begin

        @test isconcretetype(typeof(makeenv()))

        @testset "time consistency" begin
            e = makeenv()
            @test time(e) isa Float64
            @test time(e) == 0 # time should always start from 0
            t1 = time(e)
            @test time(e) === t1
            @test step!(e) === e
            t2 = time(e)
            @test (t2 - t1) == timestep(e)
        end

        # For each of state/action/observation/reward, test:
        # 1. return type of getfoo(env) is compatible with foospace(env)
        # 2. getfoo!(x, env) returns x
        # 3. getfoo!(x, env) == getfoo(env)
        # 4. if settable, that setfoo!(env, x) returns env
        # 5. if settable, that getfoo(env) == x after setfoo!(env, x)
        # 6. all functions are inferrable

        @testset "state" begin
            let e = makeenv()
                @test statespace(e) isa Shapes.AbstractVectorShape
                @test eltype(getstate(e)) == eltype(statespace(e))
                @test eltype(getstate!(rand(statespace(e)), e)) == eltype(statespace(e))
                @test axes(getstate(e)) == axes(statespace(e))
                @test axes(getstate!(rand(statespace(e)), e)) == axes(statespace(e))
            end

            @test let e = makeenv(), x = rand(statespace(e))
                x === getstate!(x, e)
            end
            @test let e = makeenv()
                getstate!(rand(statespace(e)), e) == getstate(e)
            end

            @test let e = makeenv()
                e === setstate!(e, rand(statespace(e)))
            end
            @test let e = makeenv(), x = randstate()
                setstate!(e, x)
                x == getstate!(rand(statespace(e)), e) == getstate(e)
            end
        end

        @testset "action" begin
            let e = makeenv()
                @test actionspace(e) isa Shapes.AbstractVectorShape
                @test eltype(getaction(e)) == eltype(actionspace(e))
                @test eltype(getaction!(rand(actionspace(e)), e)) == eltype(actionspace(e))
                @test axes(getaction(e)) == axes(actionspace(e))
                @test axes(getaction!(rand(actionspace(e)), e)) == axes(actionspace(e))
            end

            @test let e = makeenv(), x = rand(actionspace(e))
                x === getaction!(x, e)
            end
            @test let e = makeenv()
                getaction!(rand(actionspace(e)), e) == getaction(e)
            end

            @test let e = makeenv()
                e === setaction!(e, rand(actionspace(e)))
            end
            @test let e = makeenv(), x = randaction()
                setaction!(e, x)
                x == getaction!(rand(actionspace(e)), e) == getaction(e)
            end
        end

        @testset "obs" begin
            let e = makeenv()
                @test observationspace(e) isa Shapes.AbstractVectorShape
                @test eltype(getobservation(e)) == eltype(observationspace(e))
                @test eltype(getobservation!(rand(observationspace(e)), e)) == eltype(observationspace(e))
                @test axes(getobservation(e)) == axes(observationspace(e))
                @test axes(getobservation!(rand(observationspace(e)), e)) == axes(observationspace(e))
            end

            @test let e = makeenv(), x = rand(observationspace(e))
                x === getobservation!(x, e)
            end
            @test let e = makeenv()
                getobservation!(rand(observationspace(e)), e) == getobservation(e)
            end
        end

        # For getreward and isdone, test that they are functions of the
        # passed in (state, action, observation) and not any internal data in the env.
        # NOTE requires that `randreset!` generating a new state with different
        # reward/isdone

        @testset "reward" begin
            let e = makeenv()
                @test rewardspace(e) isa ScalarShape
                @test typeof(getreward(e)) == eltype(rewardspace(e))
            end

            @test let e = makeenv(), s = getstate(e), a = getaction(e), o = getobservation(e)
                r1 = getreward(s, a, o, e)
                r2 = getreward(e)
                randreset!(e)
                r3 = getreward(s, a, o, e)
                r1 == r2 == r3
            end
        end

        @testset "isdone" begin
            @test let e = makeenv()
                isdone(e) isa Bool
            end

            @test let e = makeenv(), s = getstate(e), a = getaction(e), o = getobservation(e)
                d1 = isdone(s, o, e)
                d2 = isdone(e)
                randreset!(e)
                d3 = isdone(s, o, e)
                d1 == d2 == d3
            end
        end

        @testset "constructor consistency" begin
            e1, e2 = makeenv(), makeenv()
            @test getstate(e1) == getstate(e2)
            @test getaction(e1) == getaction(e2)
            @test getobservation(e1) == getobservation(e2)
            @test getreward(e1) == getreward(e2)
            @test isdone(e1) == isdone(e2)
        end

        @testset "reset" begin
            let e = makeenv()
                @test e === reset!(e)
            end

            let e1 = makeenv(), e2 = makeenv()
                setaction!(e1, rand(actionspace(e1)))
                setaction!(e2, rand(actionspace(e2)))
                for _=1:100
                    step!(e1)
                    step!(e2)
                end

                reset!(e1)
                reset!(e2)

                @test getstate(e1) == getstate(e2)
                @test getaction(e1) == getaction(e2)
                @test getobservation(e1) == getobservation(e2)
                @test getreward(e1) == getreward(e2)
                @test isdone(e1) == isdone(e2)
            end
        end

        @testset "randreset" begin
            let e = makeenv()
                @test e === randreset!(e)
                @test e === randreset!(Random.default_rng(), e)
            end

            let e1 = makeenv(), e2 = makeenv(), rng = Random.MersenneTwister()
                Random.seed!(rng, 1)
                randreset!(rng, e1)
                randreset!(rng, e2)

                @test getstate(e1) != getstate(e2)
                @test getaction(e1) == getaction(e2)
                @test getobservation(e1) != getobservation(e2)
                # TODO For some environments, resetting state may not yield
                # different reward
                # @test getreward(e1) != getreward(e2)
            end

            let e1 = makeenv(), e2 = makeenv(), rng = Random.MersenneTwister()
                Random.seed!(rng, 1)
                randreset!(rng, e1)
                Random.seed!(rng, 1)
                randreset!(rng, e2)

                @test getstate(e1) == getstate(e2)
                @test getaction(e1) == getaction(e2)
                @test getobservation(e1) == getobservation(e2)
                # TODO For some environments, resetting state may not yield
                # different reward
                # @test getreward(e1) == getreward(e2)
            end
        end

    end

    @testset "Determinism" begin
        # execute a random control _trajectory and check for repeatability
        let e1 = makeenv(), e2 = makeenv()
            actions = rand(actionspace(e1), 1000)
            reset!(e1)
            reset!(e2)
            t1 = _rollout(e1, actions)
            t2 = _rollout(e2, actions)
            reset!(e1)
            t3 = _rollout(e1, actions)
            @test t1 == t2 == t3
        end
    end
end

function testenv_allocations(etype::Type{<:AbstractEnvironment}, args...; kwargs...)
    @testset "Allocations" begin
        e = makeenv()
        s, a, o = getstate(e), getaction(e), getobservation(e)

        @test_noalloc statespace($e)
        @test_noalloc getstate!($s, $e)
        @test_noalloc setstate!($e, $s)

        @test_noalloc actionspace($e)
        @test_noalloc getaction!($a, $e)
        @test_noalloc setaction!($e, $a)

        @test_noalloc observationspace($e)
        @test_noalloc getobservation!($o, $e)

        @test_noalloc rewardspace($e)
        @test_noalloc getreward($s, $a, $o, $e)

        @test_noalloc reset!($e)
        @test_noalloc randreset!($e)
        @test_noalloc step!($e)
        @test_noalloc isdone($s, $o, $e)
        @test_noalloc time($e)
        @test_noalloc timestep($e)
    end
end

function testenv_inferred(etype::Type{<:AbstractEnvironment}, args...; kwargs...)
    @testset "Type Stability" begin
        e = makeenv()
        s, a, o = getstate(e), getaction(e), getobservation(e)

        @test_inferred statespace(e)
        @test_inferred getstate!(s, e)
        @test_inferred setstate!(e, s)

        @test_inferred actionspace(e)
        @test_inferred getaction!(a, e)
        @test_inferred setaction!(e, a)

        @test_inferred observationspace(e)
        @test_inferred getobservation!(o, e)

        @test_inferred rewardspace(e)
        @test_inferred getreward(s, a, o, e)

        @test_inferred reset!(e)
        @test_inferred randreset!(e)
        @test_inferred step!(e)
        @test_inferred isdone(s, o, e)
        @test_inferred time(e)
        @test_inferred timestep(e)
    end
end


text:
## export data for JFinM_Charts
##-----------------------------

## call data with default js variable name
function writeData(data::Array{Int, 2}, chrt::VineGraph)
    ## create JavaScript data with default name
    return writeData(data, "vineArrayData", chrt)
end

## TODO: allow additional information about data:
## factors / other grouping

## call data with given js variable name
function writeData(data::Array{Int, 2}, dataName::String,
                   chrt::VineGraph)
    ## transform vine tree to graph data Javascript variable
    layers = Copulas.linkLayers(Vine(data))

    nVars = size(data, 1)
    nodesCode = "{\n \"nodes\": [\n"
    for ii=1:nVars
        if ii < nVars
            currNode = """
{"name": $ii, "group": 1},
"""
        else
            currNode = """
{"name": $ii, "group": 1}
"""
        end
        nodesCode = [nodesCode, currNode]
    end
    nodesCode = [nodesCode, "],\n"]
    linksCode = "\"links\": [\n"
    nVars = size(data, 1)
    for ii=1:nVars
        for jj=(ii+1):nVars
            if ii == (nVars-1)
                currLink = """
{"source": $(ii-1), "target": $(jj-1), "value": $(layers[ii, jj])}
"""
            else
                currLink = """
{"source": $(ii-1), "target": $(jj-1), "value": $(layers[ii, jj])},
"""
            end
            linksCode = [linksCode, currLink]
        end
    end
    linksCode = [linksCode, "]\n}"]
    dataCode = [nodesCode, linksCode]

    dataScript = """
<script>
var $dataName = $(dataCode...)
</script>
"""

    return dataScript, dataName
end


text:
#sim_res = run_sim(FT(22),FT(1.5),FT(5),FT(10),FT(0.025*0.4),FT(32));

#sim_res1 = run_sim(FT(1.8e1), FT(8.2e-1), FT(1.5e1), FT(9.5e-4), FT(2.5), FT(2.0), FT(1.1e3), FT(4.2e-1), 48*365*1 - 52*48, 48*365, 0.35);

using StatsBase
using DataFrames
using CSV
using Parameters
using Land
using Land.Photosynthesis
using Land.CanopyLayers
using Land.PlantHydraulics
using Land.SoilPlantAirContinuum
using Land.StomataModels
using Thermodynamics

using CLIMAParameters
using CLIMAParameters:AbstractEarthParameterSet
using CLIMAParameters.Planet: grav

const FT = Float32

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

KG_H_2_MOL_S = FT(55.55 / 3600);
mpa2mm = FT(10^6/9.8);

#deltaT = FT(60*60)

K_STEFAN = FT(Stefan());

#this function creates a SPAC node and modifies it by editing the
#storage volume and PV curves

function create_moflux_node(vcmax_par, k_plant, z_soil, smc0, storage_mult,
	nslice,deltaT,alpha,n,full_can,full_angles,root_dist_par,
	canopy_pvslope, trunk_pvslope)
#parlist = convert(Array{FT,1}, [22, 1.5, 5, 1e-5, 2.5,2,700,0.45,0.38,1])
#N = 10
#istart = 1
#vcmax_par, p_crit, k_plant, k_soil, b_soil, p20, z_soil, n_soil, smc0, storage_mult = parlist

node = create_spac(OSMWang{FT}(),vcmax_par,k_plant, z_soil, FT(40),
 nslice,alpha,n,full_can,full_angles,root_dist_par);

@unpack angles, can_opt, can_rad, canopy_rt, envirs, f_SL, ga, in_rad,
		latitude, leaves_rt, n_canopy, photo_set, plant_hs, plant_ps,
		rt_con, soil_opt, stomata_model, wl_set, n_root, soil_bounds = node;

#don't need to do the following two lines now that air layers are fixed
#node.plant_hs.branch[1].v_maximum = node.plant_hs.branch[2].v_maximum
#node.plant_hs.branch[1].v_storage = node.plant_hs.branch[2].v_storage

#old
canopyPV = PVCurveLinear{Float32}(canopy_pvslope, FT(1e-4));
trunkPV = PVCurveLinear{Float32}(trunk_pvslope, FT(1e-4));

for i_can in 1:n_canopy
	plant_hs.leaves[i_can].pv = canopyPV
	plant_hs.branch[i_can].pv = canopyPV
end

plant_hs.trunk.pv = trunkPV;

for i_root in 1:n_root
	plant_hs.roots[i_root].pv = trunkPV
end

total_water_vol = FT(storage_mult * 12 * node.ga * 1000/18.02); #in mol over the tree's area

#actually the branch volumes should be normalized by each branch's share of total branch length
for i in 1:n_canopy
	plant_hs.leaves[i].v_maximum = FT(0.1)*total_water_vol / n_canopy /  plant_hs.leaves[i].area
	plant_hs.branch[i].v_maximum[1] = FT(0.1)*total_water_vol / n_canopy
end

plant_hs.trunk.v_maximum[1] = FT(0.8)*total_water_vol

for i in 1:n_root
	#plant_hs.roots[i].v_maximum[1] = FT(0.05)*total_water_vol / n_root
	plant_hs.roots[i].v_maximum[1] = FT(0.05)*total_water_vol*FT(1/8) #* plant_hs.roots[i].area
end

w_soil = FT(smc0);
node.swc = ones(FT, length(node.swc))*w_soil;
ψ_soil = soil_p_25_swc(plant_hs.roots[1].sh, w_soil);
for root in plant_hs.roots
	root.p_ups = ψ_soil;
end;

#set the plant water storage so the corresponding xylem potential is slightly below the soil water potential
#pot_init = 1+0.05*ψ_soil
#pot_init = FT(0.5)
pot_init = ψ_soil - 0.025

for i in 1:n_canopy
#pot = volumes[1][i] ./ leafI.v_maximum .- FT(1)) / leafI.pv.slope
	plant_hs.leaves[i].v_storage = (FT(pot_init)*plant_hs.leaves[i].pv.slope + FT(1)) * plant_hs.leaves[i].v_maximum
	plant_hs.branch[i].v_storage = (FT(pot_init)*plant_hs.branch[i].pv.slope + FT(1)) * plant_hs.branch[i].v_maximum
end

plant_hs.trunk.v_storage = (FT(pot_init)*plant_hs.trunk.pv.slope + FT(1)) * plant_hs.trunk.v_maximum

for i in 1:n_root
	plant_hs.roots[i].v_storage = (FT(pot_init)*plant_hs.roots[i].pv.slope + FT(1)) * plant_hs.roots[i].v_maximum
end

#then calculate the steady state plant water storage
#and set plant state to that
update_pk_tree!(plant_hs)

volumes0 = get_v_prof2(plant_hs);
dmat0, dvec0 = create_deriv_mat(plant_hs)
set_vol!(plant_hs, -inv(dmat0)*dvec0)

update_pk_tree!(plant_hs)

return node

end


text:
abstract AbstractOptimizer

immutable SGD <: AbstractOptimizer
end

immutable Adam <: AbstractOptimizer
end


text:
using PointSampler
using Test

@testset "PointSampler.jl" begin
    include("monte_carlo_th.jl");
    include("monte_carlo_th_iterator.jl");
end


text:
using TensorKits
using TensorKits: spin

# Create random iMPS
const BOND = 50

const mat = begin
    SS = spin([(1,"xx"), (1,"yy"), (1,"zz")], D=3)
    hamiltonian = SS + 1/3 * SS^2
    exp(-hamiltonian * 0.1)
end

# Create AKLT Hamiltonian and iTEBD engine
a = rand(ComplexF64, 1, 3, 1)
b = rand(ComplexF64, 1, 3, 1)
l = [1.0]

# Setup TEBD
for i=1:1000
    global BOND, mat, a, b, l
    for i=1:1000
        a,b,l = itebd2!(mat, a, b, l, bound=BOND)
    end
end

# Exact AKLT ground state
aklt = begin
    aklt_tensor = zeros(2,3,2)
    aklt_tensor[1,1,2] = +sqrt(2/3)
    aklt_tensor[1,2,1] = -sqrt(1/3)
    aklt_tensor[2,2,2] = +sqrt(1/3)
    aklt_tensor[2,3,1] = -sqrt(2/3)
    aklt_tensor
end


text:
"""
    request(method, cred, path; query_kwargs...)

Send a request to Trello and get back the reponse.
 - `method` is the HTTP method to use, e.g `"GET"`, `"POST"`, `"PUT"` or `"DELETE"`
    the `get_request`, `post_request`, `put_request` and `delete_request` are convenience wrapppers of `request` with the relevant method.
 - `cred`: a `TrelloCred` object for the account making the request.
 - `path` the endpoint being targetted, e.g. `"/1/boards/5d721712fa52696736da60cf/lists"`.
    Often you will need to interpolate a board or list id or similar into this. (like `5d721712fa52696736da60cf`)
 - `query_kwargs`: the query to be used with the request,
    e.g `idList="5d72175351f2b73b0f872e03", name="My new Card", desc="The card description", ...`

When invoking `request` it is important to look at the Trello API docs to workout
what arguements are allowed, and how they are spelled etc.
https://developers.trello.com/reference

Returns the parsed response.
Generally as a `NamedTuple`, or a `Vector` of `NamedTuples`

Throws errors if there are HTTP errors.
"""
function request(method, cred, path; query_kwargs...)
    query = ["key"=>cred.api_key, "token"=>cred.api_token]
    for (key, value) in query_kwargs
        isempty(value) && continue
        push!(query, string(key) => string(value))
    end
    uri = HTTP.URI(host="api.trello.com", scheme="https", path=path, query=query)
    resp = HTTP.request(method, uri)
    return JSON2.read(String(resp.body))
end

get_request(cred, path; query_kwargs...) = request("GET", cred, path; query_kwargs...)
post_request(cred, path; query_kwargs...) = request("POST", cred, path; query_kwargs...)
put_request(cred, path; query_kwargs...) = request("PUT", cred, path; query_kwargs...)
delete_request(cred, path; query_kwargs...) = request("DELETE", cred, path; query_kwargs...)

function indexed_collection(collection, index=:name)
    keys = getproperty.(collection, index)
    return LittleDict(keys, collection)
end

sortbypos(collection) = sort(collection; by=x->x.pos)
filteroutclosed(collection) = filter(x->!x.closed, collection)


text:
module boundaries

#
# Perfect conductor
#

function perfect_conductor!(field, pos=1)
    field[pos] = 0
end

#
# ABC
#

# Trivial
function trivial_abc!(field, pos=200, right=true)
  if right
    field[pos] = field[pos-1]
  else
    field[pos] = field[pos+1]
  end
end

# Mur ABC, First order

type MurABCFirst
    right::Bool
    coef::Real
    prev::Real
    pos::Int
end

function setup_first_order_abc(eps, mu, pos=200, right=true)
    abc = MurABCFirst( right, 0.0, 0.0, pos);
    if right
        temp = 1. / sqrt(mu[pos - 1] * eps[pos])
        abc.coef = (temp - 1.0) / (temp + 1.0)
    else
        temp = 1. / sqrt(mu[pos] * eps[pos])
        abc.coef = (temp - 1.0) / (temp + 1.0)
    end
    return abc
end

function first_order_diff_abc!(field, abc)
    if abc.right
        field[abc.pos] = abc.prev + abc.coef * (field[abc.pos - 1] - field[abc.pos]);
        abc.prev = field[abc.pos - 1];
    else
        field[abc.pos] = abc.prev + abc.coef * (field[abc.pos + 1] - field[abc.pos]);
        abc.prev = field[abc.pos + 1];
    end
end

#
# CPML
#
### TODO
type CPML
    right::Bool;
end

function setup_CPML()
end

function CPML(ez, hy, cpml)
end

#
# TFSF
#

function totalfield_scatteredfield!(ez, hy, eps, mu, incident, pos=1)
    ez_inc, hy_inc = incident;
    hy[pos] -= hy_inc / globals.imp0;
    ez[pos+1] += ez_inc / sqrt(eps[i+1] * mu[i+1]);
end

end


text:
export projector, print_blocktree
using InteractiveUtils

"""
    projector(x)

Return projector on `0` or projector on `1`.
"""
projector(x) = code == 0 ? mat(P0) : mat(P1)

"""
    print_subtypetree(::Type[, level=1, indent=4])

Print subtype tree, `level` specify the depth of the tree.
"""
function print_subtypetree(t::Type, level = 1, indent = 4)
    level == 1 && println(t)
    for s in subtypes(t)
        println(join(fill(" ", level * indent)) * string(s))
        print_subtypetree(s, level + 1, indent)
    end
end

"""
    rmlines(ex)

Remove `LineNumberNode` from an `Expr`.
"""
rmlines(ex::Expr) = begin
    hd = ex.head
    hd == :macrocall && return ex
    tl = map(rmlines, filter(!islinenumbernode, ex.args))
    Expr(hd, tl...)
end
rmlines(@nospecialize(a)) = a
islinenumbernode(@nospecialize(x)) = x isa LineNumberNode


text:
const MAX_RETX_ATTEMPTS = 8

remaining_outbound_capacity(conn::TCPConnection) = remain_cap(conn.sender.stream)

bytes_in_flight(conn::TCPConnection) = bytes_in_flight(conn.sender)

unassembled_bytes(conn::TCPConnection) = unassembled_bytes(conn.receiver)

time_since_last_segment_received(conn::TCPConnection) = conn.time_since_last_segment_received

isactive(conn::TCPConnection) = conn.active

inbound_stream(conn::TCPConnection) = conn.receiver.reassembler.output

function segment_received!(conn::TCPConnection, seg::TCPSegment)
    !isactive(conn) && return
    conn.time_since_last_segment_received = 0
    # STATE: LISTEN => SYN_RECV
    # receive SYN / send SYN + ACK
    if ackno(conn.receiver) === nothing && next_seqno_absolute(conn.sender) == 0
        !seg.header.syn && return
        segment_received!(conn.receiver, seg)
        connect!(conn)
        return
    end
    # STATE: SYN_SENT => ESTABLISHED
    # receive SYN + ACK / send ACK
    if ackno(conn.receiver) === nothing && next_seqno_absolute(conn.sender) != 0 &&
        bytes_in_flight(conn.sender) == next_seqno_absolute(conn.sender)
        length(seg.payload) != 0 && return
        if seg.header.syn && !seg.header.ack
            segment_received!(conn.receiver, seg)
            send_empty_segment!(conn.sender)
            send_segs!(conn)
            return
        end
        if seg.header.rst
            set_error(conn.sender.stream)
            set_error(conn.receiver |> stream_out)
            conn.active = false
            return
        end
    end
    # give the seg to receiver
    segment_received!(conn.receiver, seg)
    # tell the sender ackno and window_size
    seg.header.ack && ack_received!(conn.sender, seg.header.ackno, seg.header.win)
    if seg.header.rst
        send_empty_segment!(conn.sender)
        unclean_shutdown!(conn)
        return
    end

    # send at least one segment to reply
    # make sure there is at least one segment
    isempty(conn.sender.segments_out) && (length_in_sequence_space(seg) != 0) &&
        send_empty_segment!(conn.sender)
    # send segments
    send_segs!(conn)
end

function write!(conn::TCPConnection, data::AbstractString)
    write_size = write!(conn.sender.stream, data)
    fill_window!(conn.sender)
    send_segs!(conn)
    write_size
end

function tick!(conn::TCPConnection, ms_since_last_tick::Int)
    !conn.active && return
    conn.time_since_last_segment_received += ms_since_last_tick
    tick!(conn.sender, ms_since_last_tick)
    (consecutive_retransmissions(conn.sender) > MAX_RETX_ATTEMPTS) && unclean_shutdown!(conn)
    send_segs!(conn)
end

function send_segs!(conn::TCPConnection)
    while conn.sender |> segments_out |> !isempty
        seg = dequeue!(conn.sender.segments_out)
        if ackno(conn.receiver) !== nothing
            seg.header.ack = true
            seg.header.ackno = ackno(conn.receiver)
            seg.header.win = window_size(conn.receiver)
        end
        enqueue!(conn.segments_out, seg)
    end
    clean_shutdown!(conn)
end

function end_input_stream!(conn::TCPConnection)
    end_input!(conn.sender.stream)
    fill_window!(conn.sender)
    send_segs!(conn)
end

function connect!(conn::TCPConnection)
    fill_window!(conn.sender)
    send_segs!(conn)
end

function unclean_shutdown!(conn::TCPConnection)
    set_error(conn.sender.stream)
    set_error(conn.receiver |> stream_out)
    conn.active = false
    seg = dequeue!(conn.sender.segments_out)
    if ackno(conn.receiver) !== nothing
        seg.header.ack = true
        seg.header.ackno = ackno(conn.receiver)
        seg.header.win = window_size(conn.receiver)
    end
    seg.header.rst = true
    enqueue!(conn.segments_out, seg)
end

function clean_shutdown!(conn::TCPConnection)
    !input_ended(stream_out(conn.receiver)) && return
    if !eof(conn.sender.stream)
        conn.linger_after_streams_finish = false
    elseif bytes_in_flight(conn.sender) == 0
        if !conn.linger_after_streams_finish ||
            conn.time_since_last_segment_received >= 10 * conn.rt_timeout
            conn.active = false
        end
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


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

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

# ╔═╡ 1dfd0480-6b0d-11eb-1f62-dbdc1a677571
begin
#	import Pkg;
#	Pkg.add("Plots")
#	Pkg.add("PlutoUI")
	using PlutoUI
	using Plots
	using LinearAlgebra
	using SparseArrays
end

# ╔═╡ e6375832-6b29-11eb-38b2-7582cac61e64
md"""

# **The strange dice** - _The central limit theorem_

**Bernoulli and Laplace** invented this interesting dice game with colored faces (red for odd numbers and green for even numbers of pips).
The side facing the gambler works as a bonus / malus factor.

In the case of a green side facing the gambler, the points get doubled, in case of a red face the points get subtracted by minus 1.

At the trinautic tournament the whole crew is rolling the strange dice and the points are summed up.

Look at the distribution below that shows the probabilities for the different number of points.

**Use the slider to add crew members to play and look how the average distribution of points (total points divided by the number of gamblers) changes**

"""

# ╔═╡ 86b1d820-6b2b-11eb-19cf-6390ab5fccc8
HTML("

")

# ╔═╡ 9f155e80-6be6-11eb-22c1-1b4f6e30ec72
md"Choose number of 👉 $(@bind iterations Slider(1:60, show_value = true)) crew members that roll the strange dice."

# ╔═╡ 7ecb15e0-6b26-11eb-1d2b-37dd4e07ef05
HTML("

")

# ╔═╡ 185ef260-6b0b-11eb-1fd8-b5c70ec616d7
p = [1/6 1/6 1/6 1/6 1/6 1/6]

# ╔═╡ 943dea30-6b0b-11eb-19f0-75227e40a492
A = [ 3/4 1/4; 1/2 1/2; 3/4 1/4; 1/4 3/4; 1/2 1/2; 1/4 3/4]

# ╔═╡ 036a9bb0-6b0c-11eb-0ece-d54174c4d7a0
probability_table = A.*p'

# ╔═╡ 3eb465e0-6b14-11eb-19a4-df75db84f12d
sum(probability_table)

# ╔═╡ 427d9fb0-6b10-11eb-043c-653dbc050130
begin
	B = zeros(13,2)
	B[(1:6)',2] = probability_table[:,2]
	B[:,1] = (0:12)'
	index = collect(1:6)*2
	B[index.+1,2] = B[index.+1,2] + probability_table[:,1]
end

# ╔═╡ 29969590-6b12-11eb-11d7-9d831132fe4f
begin
	plot(
    B[:,1], B[:,2],
    line = (1.0, 0.0, :bar),
    normalize = false,
    bins = 10,
	bar_width = 0.1,
    marker = (6, 0.5, :x),
    markerstrokewidth = 5.,
    color = [:steelblue],
    fill = 0.9,
    orientation = :v,
    title = "The strange dice distribution",
	ylabel = "probability mass function",
	xlabel = "Points",
	label = :none,
)
plot!( B[:,1], B[:,2],
    line = (1.0, 0, :path),
    normalize = false,
    bins = 10,
	bar_width = 0.2,
    marker = (4, 0.5, :o),
    markerstrokewidth = 5.,
    color = [:steelblue],
    fill = 0.9,
    orientation = :v,
    title = "The strange dice distribution",
	ylabel = "probability mass function",
	xlabel = "Points",
	label = :none,
)
end

# ╔═╡ 6df0d2c0-6b10-11eb-04d1-1b71b887b8bb
B,sum(B[:,2])

# ╔═╡ 462e27c0-6b14-11eb-2b2f-3bd578fb8712
dice_mean = sum(B[:,2].*B[:,1])

# ╔═╡ 691d5aa0-6b71-11eb-2f53-db4a2217c8e0
dice_std = sqrt(sum(B[:,2].* (B[:,1].-dice_mean).^2))

# ╔═╡ 18e26640-6be7-11eb-20b0-5542c5620d8d
iter = iterations-1

# ╔═╡ 7dbbfe90-6b20-11eb-39c8-c99fd55ccf05
sparse_vector = repeat(B[:,2]',iter*12+1)

# ╔═╡ 80f0ce50-6b21-11eb-2b31-6982188b8a41
sparse_vector[:,1]

# ╔═╡ e107cae0-6b1d-11eb-3ff5-b7861639d659
P = Array(spdiagm(0 => sparse_vector[:,1], 1 => sparse_vector[:,2], 2 => sparse_vector[:,3], 3 => sparse_vector[:,4], 4 => sparse_vector[:,5], 5 => sparse_vector[:,6], 6 => sparse_vector[:,7], 7 => sparse_vector[:,8], 8 => sparse_vector[:,9], 9 => sparse_vector[:,10], 10=> sparse_vector[:,11], 11 => sparse_vector[:,12], 12 => sparse_vector[:,13]))

# ╔═╡ 58fe1a00-6b22-11eb-207c-095c2b30a0e6
pp = append!(B[:,2], zeros((iter)*12))

# ╔═╡ b3279d30-6b22-11eb-0658-dbe663f09558
p_final = Array(pp'*P^(iter))

# ╔═╡ 6f072a20-6b23-11eb-0add-8b62821b708a
r = Array(LinRange(0,12,(iter-1)*12+13*2-1))

# ╔═╡ b7a9ad70-6b23-11eb-295b-4324ef257047
begin

	plot(
    [B[:,1], r, [dice_mean, dice_mean]], [B[:,2]./0.9,p_final[1,:]./((r[2]-r[1])*0.9), [0,1]],
    line = ([1 0.3 1], [0 0.0 3], [:bar :bar :line]),
	linestyle = [:auto :auto :dash],
    normalize = true,
    marker = (6, 0.8, [:none :none]),
    markerstrokewidth = 5.,
    color = [:steelblue :orangered :black],
    fill = [1. 1.],
	bar_width = [1*0.9 (r[2]-r[1])*0.9],
    orientation = :v,
    title = string("Average distribution of \nstrange dice with ", string(iter+1), " gamblers"),
	ylabel = "probability / width of bar",
	xlabel = "Points",
	usetex = true,
	label = [:none :none "Mean"],
	ylim = [0,1],
)
	plot!([dice_mean - dice_std, dice_mean + dice_std, 3, dice_mean+dice_std, dice_mean-dice_std], 0.5.*[1, 1,NaN, 1, 1],
		line = (2., 2., :path),
		color = :magenta,
		arrow = true,
		label = "Std",
		)
	plot!([dice_mean - dice_std/sqrt(iter), dice_mean + dice_std/sqrt(iter), 3, dice_mean+dice_std/sqrt(iter), dice_mean-dice_std/sqrt(iter)], 0.08*[1, 1,NaN, 1, 1]/sqrt(((r[2]-r[1]))),
		line = (2., 2., :path),
		color = :green,
		arrow = true,
		label = string("Std/sqrt(",string(iter+1), ")"),
		)

end

# ╔═╡ 89f8f100-6b24-11eb-3292-c15016c60a5b
p_final[1,:]

# ╔═╡ cec01d8e-6be2-11eb-1465-ad0f3ffcc2f5
function pretty(M::Matrix{T} where T<:String)
	max_length = maximum(length.(M))
	dv="
"
	HTML(dv*join([join("
".*M[i,:].*"
", " ") for i in 1:size(M,1)]
			, "
$dv")*"
")
end

# ╔═╡ e7e029f0-6b9c-11eb-179d-1d2fca9eb8af
[string.(rand("🧔👩👨👧🧑", 1, iterations));string.(rand("🎲", 1, iterations))] |> pretty

# ╔═╡ Cell order:
# ╟─1dfd0480-6b0d-11eb-1f62-dbdc1a677571
# ╟─e6375832-6b29-11eb-38b2-7582cac61e64
# ╟─86b1d820-6b2b-11eb-19cf-6390ab5fccc8
# ╟─29969590-6b12-11eb-11d7-9d831132fe4f
# ╟─e7e029f0-6b9c-11eb-179d-1d2fca9eb8af
# ╟─9f155e80-6be6-11eb-22c1-1b4f6e30ec72
# ╟─b7a9ad70-6b23-11eb-295b-4324ef257047
# ╟─7ecb15e0-6b26-11eb-1d2b-37dd4e07ef05
# ╟─185ef260-6b0b-11eb-1fd8-b5c70ec616d7
# ╟─943dea30-6b0b-11eb-19f0-75227e40a492
# ╟─036a9bb0-6b0c-11eb-0ece-d54174c4d7a0
# ╟─3eb465e0-6b14-11eb-19a4-df75db84f12d
# ╟─427d9fb0-6b10-11eb-043c-653dbc050130
# ╟─6df0d2c0-6b10-11eb-04d1-1b71b887b8bb
# ╟─691d5aa0-6b71-11eb-2f53-db4a2217c8e0
# ╟─462e27c0-6b14-11eb-2b2f-3bd578fb8712
# ╟─18e26640-6be7-11eb-20b0-5542c5620d8d
# ╟─7dbbfe90-6b20-11eb-39c8-c99fd55ccf05
# ╟─80f0ce50-6b21-11eb-2b31-6982188b8a41
# ╟─e107cae0-6b1d-11eb-3ff5-b7861639d659
# ╟─58fe1a00-6b22-11eb-207c-095c2b30a0e6
# ╟─b3279d30-6b22-11eb-0658-dbe663f09558
# ╟─6f072a20-6b23-11eb-0add-8b62821b708a
# ╟─89f8f100-6b24-11eb-3292-c15016c60a5b
# ╟─cec01d8e-6be2-11eb-1465-ad0f3ffcc2f5

max_stars_repo_name: fonsp/bayes-imoox

max_issues_repo_name: fonsp/bayes-imoox

max_forks_repo_name: fonsp/bayes-imoox


text:
using DeepUncertainty
include("layers.jl")

function CGCNN(
    input_feature_length;
    num_conv = 2,
    conv_activation = softplus,
    atom_conv_feature_length = 80,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    num_hidden_layers = 1,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
)
    @assert atom_conv_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(
        AGNConv(
            input_feature_length => atom_conv_feature_length,
            conv_activation,
            initW = initW,
        ),
        [
            AGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
        [
            Dense(
                pooled_feature_length,
                pooled_feature_length,
                hidden_layer_activation,
                init = initW,
            ) for i = 1:num_hidden_layers-1
        ]...,
        Dense(pooled_feature_length, output_length, output_layer_activation, init = initW),
    )
    return model
end

function MC_CGCNN(
    input_feature_length;
    num_conv = 2,
    conv_activation = softplus,
    atom_conv_feature_length = 80,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    num_hidden_layers = 1,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
    dropout_rate = 0.2,
)
    @assert atom_conv_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(
        MCAGNConv(
            input_feature_length => atom_conv_feature_length,
            dropout_rate,
            conv_activation,
            initW = initW,
        ),
        [
            MCAGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                dropout_rate,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
        [
            MCLayer(
                MCDense(
                    pooled_feature_length,
                    pooled_feature_length,
                    dropout_rate,
                    hidden_layer_activation,
                ),
                dropout_rate,
            ) for i = 1:num_hidden_layers-1
        ]...,
        MCDense(
            pooled_feature_length,
            output_length,
            dropout_rate,
            output_layer_activation,
        ),
    )
    return model
end

function BayesianCGCNN(
    input_feature_length;
    num_conv = 2,
    conv_activation = softplus,
    atom_conv_feature_length = 80,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    num_hidden_layers = 1,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
)
    @assert atom_conv_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(
        VariationalAGNConv(
            input_feature_length => atom_conv_feature_length,
            conv_activation,
            initW = initW,
        ),
        [
            VariationalAGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
        [
            VariationalDense(
                pooled_feature_length,
                pooled_feature_length,
                hidden_layer_activation,
                init = initW,
            ) for i = 1:num_hidden_layers-1
        ]...,
        VariationalDense(
            pooled_feature_length,
            output_length,
            output_layer_activation,
            init = initW,
        ),
    )
    return model
end


text:
# Defines structus and functions relating to display.
# Defines `render` for all game objects.

function renderProgressBar(percent, cam::Camera, renderer, center::AbstractPos{C}, dims::D, color, bgColor, boxColor) where D<:AbstractDims{C} where C
    # bg
    renderRectCentered(cam, renderer, center, dims, bgColor)
    # progress
    origin = topLeftPos(center, dims)
    renderRectFromOrigin(cam, renderer, origin, D(round(percent * dims.w), dims.h), color)
    # outline
    renderRectCentered(cam, renderer, center, dims, nothing; outlineColor=boxColor)
end
function renderUnit(o::UnitTypes, playerColor, cam::Camera, renderer, dims::WorldDims, color)
    # First render the player color, then the unit color (for transparency)
    renderRectCentered(cam, renderer, o.pos, dims, playerColor)
    renderRectCentered(cam, renderer, o.pos, dims, color)

    # render health bar
    healthBarPos = WorldPos(o.pos.x, o.pos.y + dims.h/2 + healthBarRenderOffset)
    healthBarDims = WorldDims(healthBarRenderWidth, healthBarRenderHeight)
    renderProgressBar(health_percent(o), cam, renderer, healthBarPos,
          healthBarDims, healthBarColor, kBuildOpsBgColor,
          if toScreenPixelDims(healthBarDims, cam).h <= 2 # pixels
              nothing
          else
              healthBarOutlineColor
          end)
end
function render(o::Collector, playerColor, cam::Camera, renderer)
    unitColor = blendAlphaColors(playerColor, kCollectorColor)
    dims = WorldDims(collectorRenderWidth, collectorRenderWidth)
    tl = topLeftPos(o.pos, dims)
    for y in [tl.y, tl.y - dims.h]
        for x in tl.x+1 : 5 : tl.x-1 + dims.w
            nudge = Vector2D(rand(0.0:0.01:2.0,2)...)
            renderRectCentered(cam, renderer, WorldPos(x,y)+nudge, WorldDims(2,10), unitColor)
        end
    end
    for x in [tl.x, tl.x + dims.w]
        for y in tl.y-1 : -5 : tl.y+1 - dims.h
            nudge = Vector2D(rand(0.0:0.01:2.0,2)...)
            renderRectCentered(cam, renderer, WorldPos(x,y)+nudge, WorldDims(10,2), unitColor)
        end
    end
    renderUnit(o, playerColor, cam, renderer, dims, kCollectorColor)
end
function render(o::Fighter, playerColor, cam::Camera, renderer)
    dims = WorldDims(unitRenderWidth, unitRenderWidth)
    renderUnit(o, playerColor, cam, renderer, dims, kFighterColor)
end

abstract type AbstractButton end
mutable struct MenuButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    dims::UIPixelDims
    text::String
    callBack
end
mutable struct KeyButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    dims::UIPixelDims
    text::String
    callBack
end

mutable struct CheckboxButton
    toggled::Bool
    button::MenuButton
end

import Base.run
run(b::AbstractButton) = b.callBack()
function run(b::CheckboxButton)
    b.toggled = !b.toggled
    b.button.callBack(b.toggled)
end

function render(b::AbstractButton, cam::Camera, renderer, color, fontSize)
    if (!b.enabled)
         return
    end
    topLeft = topLeftPos(b.pos, b.dims)
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos..., toScreenPixelDims(b.dims, cam)...)
    x,y = Int[0], Int[0]
    SDL2.GetMouseState(pointer(x), pointer(y))
    if clickedButton == b
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 50
        else
            color = color - 30
        end
    else
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 10
        end
    end
    SetRenderDrawColor(renderer, color)
    SDL2.RenderFillRect(renderer, Ref(rect) )
    renderText(renderer, cam, b.text, b.pos; fontSize = fontSize)
end

function render(b::MenuButton, cam::Camera, renderer)
    render(b, cam, renderer, kMenuButtonColor, kMenuButtonFontSize)
end
function render(b::KeyButton, cam::Camera, renderer)
    render(b, cam, renderer, kKeySettingButtonColor, kKeyButtonFontSize)
end

function render(b::CheckboxButton, cam::Camera, renderer)
    # Hack: move button text offcenter before rendering to accomodate checkbox
    offsetText = " "
    text_backup = b.button.text
    b.button.text = offsetText * b.button.text
    render(b.button, cam, renderer)
    b.button.text = text_backup

    # Render checkbox
    render_checkbox_square(b.button, 6, SDL2.Color(200,200,200, 255), cam, renderer)

    if b.toggled
        # Inside checkbox "fill"
        render_checkbox_square(b.button, 8, SDL2.Color(100,100,100, 255), cam, renderer)
    end
end

function render_checkbox_square(b::AbstractButton, border, color, cam, renderer)
    checkbox_radius = b.dims.h/2.0 - border  # (checkbox is a square)
    topLeft = topLeftPos(b.pos, b.dims)
    topLeft = topLeft + border
    screenPos = toScreenPos(topLeft, cam)
    screenDims = toScreenPixelDims(UIPixelDims(checkbox_radius*2, checkbox_radius*2), cam)
    rect = SDL2.Rect(screenPos..., screenDims...)
    SetRenderDrawColor(renderer, color)
    SDL2.RenderFillRect(renderer, Ref(rect) )
end


text:
# Beg et al. (2007): https://doi.org/10.1073/pnas.0609845104.
module BegData

    using DataFrames
    import CSV

    using ProjAssistant
    @gen_sub_proj

    ## ------------------------------------------------------------------------
    # FILES AND DIRS
    # const BEG_RAW_DATA_DIR = joinpath(ChF.RAW_DATA_DIR, "beg2007___data")
    # const BEG_ENZIMATIC_DATA_FILE =
    #     joinpath(BEG_RAW_DATA_DIR, "beg2007___enzymatic_data.tsv")

    # The average crowding coefficient a was fitted to obtain the minimum square
    # deviation between the measured and model predicted growth rates,
    # resulting in a = 0.0040 ± 0.0005 h•g/mmol, in which g is grams dry weight.
    # However, the maximum growth rates on glucose and glycerol are more consistent
    # with a = 0.0031 ± 0.0001 h•g/mmol and a = 0.0053 ± 0.0001 h•g/mmol, respectively.
    const AVE_COST  = 0.0031 # I will use this to set the unknown costs

    ## ------------------------------------------------------------------------
    # enzymatic costs from Beg et al. (2007):
    # https://doi.org/10.1073/pnas.0609845104.
    load_enz_data() = CSV.read(rawdir(BegData, "beg2007___enzymatic_data.tsv"), DataFrame)

    function __init__()
        @create_proj_dirs
    end
end