text:
using POMDPs, Crux, Flux, POMDPGym, Random, Distributions, POMDPPolicies, Plots
Crux.set_function("isfailure", POMDPGym.isfailure)

# Define the disturbance distribution based on a normal distribution
xnom = Normal(0f0, 0.5f0)
xs = Float32[-2., -0.5, 0, 0.5, 2.]
ps = exp.([logpdf(xnom, x) for x in xs])
ps ./= sum(ps)
px = DiscreteNonParametric(xs, ps)
xlogprobs = Base.log.(ps)

# Action space of the protagonist
as = Float32[-2., -0.5, 0, 0.5, 2.]

# Construct the MDP
mdp = AdditiveAdversarialMDP(InvertedPendulumMDP(actions=as, λcost=0, include_time_in_state=true), px)
S = state_space(mdp)

# construct the model
QS(outputs) = DiscreteNetwork(Chain(Dense(3, 64, relu), Dense(64, 64, relu), Dense(64, length(outputs))), outputs)
Pf(outputs) = DiscreteNetwork(Chain(Dense(3, 64, relu), Dense(64, 64, relu), Dense(64, length(outputs), (x)->-softplus(-(x-2)))), outputs)
# Pf(outputs) = DiscreteNetwork(Chain(Dense(3, 256, tanh), Dense(256, 256, tanh), Dense(256, length(outputs), sigmoid)), outputs, (x) -> x ./ sum(x, dims=1))
AdvPol(protag = QS(as)) = AdversarialPolicy(protag, Pf(xs))

# Solve with DQN
𝒮_dqn = DQN(π=QS(as), S=S, N=50000, buffer_size=Int(1e4), buffer_init=1000, required_columns=[:fail])
π_dqn = solve(𝒮_dqn, mdp)

# show the nominal distriubtion of paths
D = steps!(Sampler(mdp, 𝒮_isarl.π.P, S=S, required_columns=[:fail]), Nsteps=10000)
D2 = steps!(Sampler(mdp, π_dqn, S=S, required_columns=[:fail]), Nsteps=10000)
scatter(D[:s][1, :], D[:s][3, :], marker_z = D[:done][1,:], xlabel="θ", ylabel="t", alpha=0.5)
scatter(D2[:s][1, :], D2[:s][3, :], marker_z = D2[:done][1,:], xlabel="θ", ylabel="t", alpha=0.5)

# solve with IS
𝒮_isarl = ISARL_Discrete(π=AdvPol(),
                         S=S,
                         ϵ_init = 1e-5,
                         N=100000,
                         xlogprobs=xlogprobs,
                         px=px,
                         buffer_size=100_000,
                         buffer_init=1000,
                         c_opt = (;batch_size=128),
                         x_c_opt=(;batch_size=1024),
                         π_explore=ϵGreedyPolicy(Crux.LinearDecaySchedule(1f0, 0.1f0, 20000), as),
                         )
π_isarl = solve(𝒮_isarl, mdp)

y = Crux.IS_DQN_target(antagonist(𝒮_isarl.π),𝒮_isarl.𝒫, 𝒮_isarl.buffer, 1f0)

indices = 𝒮_isarl.buffer[:done][:]
y[indices]

value(antagonist(𝒮_isarl.π), 𝒮_isarl.buffer[:s][:,indices], 𝒮_isarl.buffer[:x][:,indices])

# check the number of failures in the buffers (compared to successs)
sum(𝒮_isarl.buffer[:fail])
sum(𝒮_isarl.buffer[:done])

# Plot the distribution of disturbances
histogram(Flux.onecold(𝒮_isarl.buffer[:x]))

# Plot the value function in 2D
v(θ, t) = sum(ps .* value(antagonist(𝒮_isarl.π), [θ, -1, t]))
heatmap(deg2rad(-20):0.01:deg2rad(20), 0:0.1:5, v)

# Show the distribution of data in the replay buffer
scatter(𝒮_isarl.buffer[:s][1, :], 𝒮_isarl.buffer[:s][3, :], marker_z = 𝒮_isarl.buffer[:done][1,:], xlabel="θ", ylabel="t", alpha=0.5)
vline!([deg2rad(20), deg2rad(-20)])

# pfail_dqn = Crux.failure(Sampler(mdp, π_dqn, S=S, max_steps=100), Neps=Int(1e5), threshold=100)
# println("DQN Failure rate: ", pfail_dqn)

# solve with RARL
# 𝒮_rarl = RARL(π=AdvPol(), S=S, N=N)
# π_rarl = solve(𝒮_rarl, mdp)

# pfail_rarl = Crux.failure(Sampler(mdp, protagonist(π_rarl), S=S, max_steps=100), Neps=Int(1e5), threshold=100)
# println("RARL Failure rate: ", pfail_rarl)

# pfail_isarl = Crux.failure(Sampler(mdp, protagonist(π_isarl), S=S, max_steps=100), Neps=Int(1e5), threshold=100)
# println("IS Failure rate: ", pfail_isarl)


text:
using LibGit2
using LibGit2.Consts: DELTA_ADDED, DELTA_MODIFIED, DELTA_DELETED
using JuliaFormatter

const DEFAULT_DIFF_FILTER = Set((
    LibGit2.Consts.DELTA_ADDED,
    LibGit2.Consts.DELTA_MODIFIED,
    # LibGit2.Consts.DELTA_DELETED # only need to format existing files
))

# modified based on LibGit2.diff_files
function diff_files(repo, tree, pathspecs=""; filter=DEFAULT_DIFF_FILTER)
    tree = LibGit2.GitTree(repo, tree * "^{tree}")
    files = String[]
    try
        diff = LibGit2.diff_tree(repo, tree, pathspecs)
        for i in 1:LibGit2.count(diff)
            delta = diff[i]
            delta === nothing && break
            if LibGit2.Consts.DELTA_STATUS(delta.status) in filter
                push!(files, unsafe_string(delta.new_file.path))
            end
        end
        close(diff)
    catch
        close(tree)
    end
    return files
end

"""
    format_diff_file(root, tree="HEAD", pathspecs=""; filter)

Apply JuliaFormatter to git diff (.jl) files.

# Inputs

- `root::String`: the root path to the git repo
- `tree::String`: git tree spec
- `pathspecs::String=""`: same to `git diff -- <pathspecs>`
- `filter=identity`: a function that used to filter out files of interested.
"""
function format_diff_file(root::String, tree="HEAD", pathspecs::String=""; filter=identity)
    repo = LibGit2.GitRepo(root)

    jl_diff = Base.filter(
        x -> endswith(x, ".jl"), joinpath.(root, diff_files(repo, tree, pathspecs))
    )
    for file in filter(jl_diff)
        format(file, BlueStyle(); verbose=true)
    end
end

ext: jl


text:
# Miscellaneous helper macros and functions
#

# Helper macro to iterate over the DG grid. Generates the needed loops
# and indices: `eh`, `ev`, `e`, `k,`, `j`, `i`, `ijk`.
macro traverse_dg_grid(grid_info, topl_info, expr)
    return esc(
        quote
            for eh in 1:($(topl_info).nhorzrealelem)
                for ev in 1:($(topl_info).nvertelem)
                    e = ev + (eh - 1) * $(topl_info).nvertelem
                    for k in 1:($(grid_info).Nqk)
                        evk = $(grid_info).Nqk * (ev - 1) + k
                        for j in 1:$(grid_info).Nq[2]
                            for i in 1:$(grid_info).Nq[1]
                                ijk =
                                    i +
                                    $(grid_info).Nq[1] *
                                    ((j - 1) + $(grid_info).Nq[2] * (k - 1))
                                $expr
                            end
                        end
                    end
                end
            end
        end,
    )
end

# Helper macro to iterate over a 3D array. Used for walking the
# interpolated (GCM) grid.
macro traverse_interpolated_grid(nlong, nlat, nlevel, expr)
    return esc(quote
        for lo in 1:($nlong)
            for la in 1:($nlat)
                for le in 1:($nlevel)
                    $expr
                end
            end
        end
    end)
end

# Helpers to extract data from the various state arrays
function extract_state(dg, state, ijk, e, st::AbstractStateType)
    bl = dg.balance_law
    FT = eltype(state)
    num_state = number_states(bl, st)
    local_state = MArray{Tuple{num_state}, FT}(undef)
    for s in 1:num_state
        local_state[s] = state[ijk, s, e]
    end
    return Vars{vars_state(bl, st, FT)}(local_state)
end


text:
using PMIx_jll

function __init__()
    # Required for PMIx relocateable binaries
    # TODO: this should be done in PMIx_jll package
    # https://github.com/JuliaPackaging/Yggdrasil/issues/390
    ENV["PMIX_INSTALL_PREFIX"] = PMIx_jll.artifact_dir
    nothing
end

const pid_t = Cint

const __time_t = Clong
const __suseconds_t = Clong
const time_t = __time_t

const UINT32_MAX = typemax(UInt32)
const UINT8_MAX = typemax(UInt8)

# FIXME:
# correct for 4.1.0

const PMIX_INTERNAL_ERR_BASE = -1330

max_issues_repo_name: JuliaParallel/PMIx

max_forks_repo_name: JuliaParallel/PMIx


text:
__precompile__()

module Espresso

export
    # utils
    ExH,
    to_block,
    parse_call_args,
    parse_call_expr,
    make_call_expr,
    with_keywords,
    without_keywords,
    # rewrite
    matchex,
    matchingex,
    findex,
    subs,
    rewrite,
    tryrewrite,
    rewrite_all,
    without,
    set_default_placeholders,
    # simplification
    simplify,
    @simple_rule,
    # funexpr
    funexpr,
    func_expr,
    func_name,
    get_or_generate_argnames,
    # indexing
    split_indexed,
    make_indexed,
    with_indices,
    get_vars,
    get_var_names,
    get_indices,
    find_vars,
    find_var_names,
    find_indices,
    forall_sum_indices,
    forall_indices,
    sum_indices,
    # ExNode
    ExNode,
    getcategory,
    getvar,
    setvar!,
    varname,
    varidxs,
    getexpr,
    getexpr_kw,
    setexpr!,
    getguards,
    setguards!,
    getvalue,
    setvalue!,
    indexof,
    dependencies,
    to_expr,
    to_expr_kw,
    isindexed,
    # ExGraph core
    AbstractExGraph,
    ExGraph,
    parse!,
    reparse,
    evaluate!,
    cat,
    fuse_assigned,
    # ExGraph utils
    dependents,
    external_vars,
    topsort,
    collect_deps,
    expand_deps,
    expand_const,
    remove_unused,
    eliminate_common,
    inline_nodes,
    graph_hash,
    # expr utils
    mergeex,
    sanitize,
    # tracking
    TrackedArray,
    TrackedReal,
    tracked_val,
    tracked_exgraph,
    swap_default_graph!,
    reset_default_graph!,
    # conversions
    to_buffered,
    to_inplace,
    @inplacerule,
    # destruct
    make_func_expr,
    isstruct,
    field_values,
    named_field_values,
    destruct,
    destruct_inputs,
    # codegens
    generate_code,
    VectorCodeGen,
    BufCodeGen,
    CuCodeGen,
    GPUCodeGen,
    autoselect_codegen,
    # helpers
    sum_1,
    sum_2,
    squeeze_sum,
    squeeze_sum_1,
    squeeze_sum_2,
    __construct,
    # re-export
    mul!

include("core.jl")

end


text:
const DS = DataStructure101

export PriorityQueue
export update!, pop_min!, pop_max!

mutable struct PriorityQueue{T}
    priority::DS.MinMaxHeap
    data::Vector{T}

    PriorityQueue{T}() where T = new(DS.MinMaxHeap(), T[])
end

function Base.show(io::IO, pq::PriorityQueue{T}) where T
    print(io, "PriorityQueue{$T}($(pq.data)")
end

function Base.push!(pq::PriorityQueue{T}, v::T) where T
    push!(pq.data, v)
    priority = length(pq.data)
    DS.push_bubble!(pq.priority, priority)
end

Base.minimum(pq::PriorityQueue) = minimum(pq.priority) => pq.data[minimum(pq.priority)]

Base.maximum(pq::PriorityQueue) = maximum(pq.priority) => pq.data[maximum(pq.priority)]

function pop_min!(pq::PriorityQueue)
    pop_min!(pq.priority)
    pq.priority.data .-= 1

    return popfirst!(pq.data)
end

function pop_max!(pq::PriorityQueue)
    pop_max!(pq.priority)

    return pop!(pq.data)
end

function update!(pq::PriorityQueue, p::Int64, offset::Int64)
    space = offset > 0 ? 1 : -1
    for i = p:space:p+(offset-space)
        pq.data[i], pq.data[i+space] = pq.data[i+space], pq.data[i]
    end
end


text:
@testset verbose=true "𝔇" begin

    @testset "Compare H/D indexing ($T)" for T in [Float64, Float32]
        # Here, we check that we can pass in either an "H wedge" array to be used with
        # WignerHindex, or a full 𝔇 array used with WignerDindex, and obtain the same
        # H recurrence results
        ℓₘₐₓ = 8
        for m′ₘₐₓ in 0:ℓₘₐₓ
            expiβ = cis(rand(0:eps(T):π))
            expiβNaNCheck = complex(NaNCheck{T}(expiβ.re), NaNCheck{T}(expiβ.im))
            NCTN = NaNCheck{T}(NaN)
            Hw = fill(NCTN, WignerHsize(ℓₘₐₓ, m′ₘₐₓ))
            H!(Hw, expiβNaNCheck, ℓₘₐₓ, m′ₘₐₓ, H_recursion_coefficients(ℓₘₐₓ, T))
            𝔇 = fill(NCTN, WignerDsize(ℓₘₐₓ, m′ₘₐₓ))
            H!(
                𝔇, expiβNaNCheck, ℓₘₐₓ, m′ₘₐₓ,
                H_recursion_coefficients(ℓₘₐₓ, T), WignerDindex
            )
            for n in 0:ℓₘₐₓ
                for m′ in -min(n, m′ₘₐₓ):min(n, m′ₘₐₓ)
                    for m in abs(m′):n
                        Hnm′m = Hw[WignerHindex(n, m′, m, m′ₘₐₓ)]
                        𝔇nm′m = 𝔇[WignerDindex(n, m′, m, m′ₘₐₓ)]
                        @test Hnm′m == 𝔇nm′m
                    end
                end
            end
        end
    end

    @testset "Compare 𝔇 to formulaic d ($T)" for T in [BigFloat, Float64, Float32]
        # Now, we're ready to check that d_{n}^{m′,m}(β) matches the expected values
        # for a range of β values
        for ℓₘₐₓ in 0:4
            H_rec_coeffs = H_recursion_coefficients(ℓₘₐₓ, T)
            𝔇 = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ, ℓₘₐₓ))
            expimα = Array{Complex{T}}(undef, ℓₘₐₓ+1)
            expimγ = Array{Complex{T}}(undef, ℓₘₐₓ+1)
            expiα = complex(one(T))
            expiγ = complex(one(T))
            for β in βrange(T)
                expiβ = cis(β)
                R = from_euler_angles(zero(T), β, zero(T))
                D!(𝔇, R, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                for n in 0:ℓₘₐₓ
                    for m′ in -n:n
                        for m in -n:n
                            𝔇_formula = ExplicitWignerMatrices.D_formula(
                                n, m′, m, expiα, expiβ, expiγ
                            )
                            𝔇_recurrence = 𝔇[WignerDindex(n, m′, m)]
                            @test 𝔇_formula ≈ 𝔇_recurrence atol=200eps(T) rtol=200eps(T)
                        end
                    end
                end
            end
        end
    end

    @testset "Compare 𝔇 to formulaic 𝔇 ($T)" for T in [BigFloat, Float64, Float32]
        # Now, we're ready to check that 𝔇_{n}^{m′,m}(β) matches the expected values
        # for a range of α, β, γ values
        Random.seed!(123)
        ℓₘₐₓ = T===BigFloat ? 4 : 8
        H_rec_coeffs = H_recursion_coefficients(ℓₘₐₓ, T)
        𝔇 = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ, ℓₘₐₓ))
        expimα = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        expimγ = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        @showprogress "Compare 𝔇 to formulaic 𝔇 ($T)" for α in αrange(T, 5)
            for β in βrange(T, 5)
                for γ in γrange(T, 5)
                    R = from_euler_angles(α, β, γ)
                    expiα, expiβ, expiγ = to_euler_phases(R)
                    D!(𝔇, R, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                    for n in 0:ℓₘₐₓ
                        for m′ in -n:n
                            for m in -n:n
                                𝔇_formula = ExplicitWignerMatrices.D_formula(
                                    n, m′, m, expiα, expiβ, expiγ
                                )
                                𝔇_recurrence = 𝔇[WignerDindex(n, m′, m)]
                                @test 𝔇_formula ≈ 𝔇_recurrence atol=30eps(T) rtol=30eps(T)
                            end
                        end
                    end
                end
            end
        end
    end

    @testset "Group characters $T" for T in [BigFloat, Float64, Float32]
        # χʲ(β) ≔ Σₘ dʲₘₘ(β) ≡ Σₘ 𝔇ʲₘₘ(exp(v̂ β/2)) = sin((2j+1)β/2) / sin(β/2)
        # Here, v̂ is any unit vector; group characters are constant on conjugacy classes and
        # conjugacy classes of SO(3) are rotations through the same angle about any axis.
        ℓₘₐₓ = T===BigFloat ? 10 : 20
        m′ₘₐₓ = ℓₘₐₓ
        H_rec_coeffs = H_recursion_coefficients(ℓₘₐₓ, T)
        d = Array{T}(undef, WignerDsize(ℓₘₐₓ, m′ₘₐₓ))
        𝔇 = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ, m′ₘₐₓ))
        expimα = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        expimγ = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        @showprogress "Group characters $T" for β in βrange(T)
            expiβ = cis(β)
            d!(d, expiβ, ℓₘₐₓ, H_rec_coeffs)
            for j in 0:ℓₘₐₓ
                sin_ratio = sin((2j+1)*β/2) / sin(β/2)
                if abs(β) < 10eps(T)
                    sin_ratio = T(2j+1)
                elseif abs(β-π) < 10eps(T)
                    sin_ratio = T(-1)^j
                end
                χʲ = sum(d[WignerDindex(j, m, m)] for m in -j:j)
                @test χʲ ≈ sin_ratio atol=500eps(T) rtol=500eps(T)
                for v̂ in v̂range(T)
                    R = exp(β/2 * v̂)
                    D!(𝔇, R, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                    χʲ = sum(𝔇[WignerDindex(j, m, m)] for m in -j:j)
                    @test χʲ ≈ sin_ratio atol=500eps(T) rtol=500eps(T)
                end
            end
        end
    end

    @testset "Representation property ($T)" for T in [Float64, Float32, BigFloat]
        # For each l, 𝔇ˡₙ,ₘ(R₁ R₂) = Σₚ 𝔇ˡₙ,ₚ(R₁) 𝔇ˡₚ,ₘ(R₂)
        tol = 3eps(T)
        ℓₘₐₓ = 10
        𝔇₁ = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ))
        𝔇₂ = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ))
        𝔇₁₂ = Array{Complex{T}}(undef, WignerDsize(ℓₘₐₓ))
        H_rec_coeffs = H_recursion_coefficients(ℓₘₐₓ, T)
        expimα = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        expimγ = Array{Complex{T}}(undef, ℓₘₐₓ+1)
        @showprogress "Representation property ($T)" for R₁ in Rrange(T)
            for R₂ in Rrange(T)
                D!(𝔇₁, R₁, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                D!(𝔇₂, R₂, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                D!(𝔇₁₂, R₁*R₂, ℓₘₐₓ, H_rec_coeffs, expimα, expimγ)
                for ℓ in 0:ℓₘₐₓ
                    i = WignerDindex(ℓ, -ℓ, -ℓ)
                    j = WignerDindex(ℓ, ℓ, ℓ)
                    𝔇₁ˡ = transpose(reshape(𝔇₁[i:j], 2ℓ+1, 2ℓ+1))
                    𝔇₂ˡ = transpose(reshape(𝔇₂[i:j], 2ℓ+1, 2ℓ+1))
                    𝔇₁₂ˡ = transpose(reshape(𝔇₁₂[i:j], 2ℓ+1, 2ℓ+1))
                    @test 𝔇₁ˡ * 𝔇₂ˡ ≈ 𝔇₁₂ˡ atol=(2ℓ+1)^2*tol rtol=(2ℓ+1)^2*tol
                end
            end
        end
    end

end


text:
module MouseButtons
export MouseButton
@enum MouseButton::UInt8 begin
    LEFT = 1
    MIDDLE = 2
    RIGHT = 3
    WHEEL_UP = 4
    WHEEL_DOWN = 5
end
end

module Keys
import SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer
export Key
@enum Key::UInt32 begin
    UNKNOWN = SDL2.SDLK_UNKNOWN;
    BACKSPACE = SDL2.SDLK_BACKSPACE;
    TAB = SDL2.SDLK_TAB;
    CLEAR = SDL2.SDLK_CLEAR;
    RETURN = SDL2.SDLK_RETURN;
    PAUSE = SDL2.SDLK_PAUSE;
    ESCAPE = SDL2.SDLK_ESCAPE;
    SPACE = SDL2.SDLK_SPACE;
    EXCLAIM = SDL2.SDLK_EXCLAIM;
    QUOTEDBL = SDL2.SDLK_QUOTEDBL;
    HASH = SDL2.SDLK_HASH;
    DOLLAR = SDL2.SDLK_DOLLAR;
    AMPERSAND = SDL2.SDLK_AMPERSAND;
    QUOTE = SDL2.SDLK_QUOTE;
    LEFTPAREN = SDL2.SDLK_LEFTPAREN;
    RIGHTPAREN = SDL2.SDLK_RIGHTPAREN;
    ASTERISK = SDL2.SDLK_ASTERISK;
    PLUS = SDL2.SDLK_PLUS;
    COMMA = SDL2.SDLK_COMMA;
    MINUS = SDL2.SDLK_MINUS;
    PERIOD = SDL2.SDLK_PERIOD;
    SLASH = SDL2.SDLK_SLASH;
    K_0 = SDL2.SDLK_0;
    K_1 = SDL2.SDLK_1;
    K_2 = SDL2.SDLK_2;
    K_3 = SDL2.SDLK_3;
    K_4 = SDL2.SDLK_4;
    K_5 = SDL2.SDLK_5;
    K_6 = SDL2.SDLK_6;
    K_7 = SDL2.SDLK_7;
    K_8 = SDL2.SDLK_8;
    K_9 = SDL2.SDLK_9;
    COLON = SDL2.SDLK_COLON;
    SEMICOLON = SDL2.SDLK_SEMICOLON;
    LESS = SDL2.SDLK_LESS;
    EQUALS = SDL2.SDLK_EQUALS;
    GREATER = SDL2.SDLK_GREATER;
    QUESTION = SDL2.SDLK_QUESTION;
    AT = SDL2.SDLK_AT;
    LEFTBRACKET = SDL2.SDLK_LEFTBRACKET;
    BACKSLASH = SDL2.SDLK_BACKSLASH;
    RIGHTBRACKET = SDL2.SDLK_RIGHTBRACKET;
    CARET = SDL2.SDLK_CARET;
    UNDERSCORE = SDL2.SDLK_UNDERSCORE;
    BACKQUOTE = SDL2.SDLK_BACKQUOTE;
    A = SDL2.SDLK_a;
    B = SDL2.SDLK_b;
    C = SDL2.SDLK_c;
    D = SDL2.SDLK_d;
    E = SDL2.SDLK_e;
    F = SDL2.SDLK_f;
    G = SDL2.SDLK_g;
    H = SDL2.SDLK_h;
    I = SDL2.SDLK_i;
    J = SDL2.SDLK_j;
    K = SDL2.SDLK_k;
    L = SDL2.SDLK_l;
    M = SDL2.SDLK_m;
    N = SDL2.SDLK_n;
    O = SDL2.SDLK_o;
    P = SDL2.SDLK_p;
    Q = SDL2.SDLK_q;
    R = SDL2.SDLK_r;
    S = SDL2.SDLK_s;
    T = SDL2.SDLK_t;
    U = SDL2.SDLK_u;
    V = SDL2.SDLK_v;
    W = SDL2.SDLK_w;
    X = SDL2.SDLK_x;
    Y = SDL2.SDLK_y;
    Z = SDL2.SDLK_z;
    DELETE = SDL2.SDLK_DELETE;
    KP0 = SDL2.SDLK_KP_0;
    KP1 = SDL2.SDLK_KP_1;
    KP2 = SDL2.SDLK_KP_2;
    KP3 = SDL2.SDLK_KP_3;
    KP4 = SDL2.SDLK_KP_4;
    KP5 = SDL2.SDLK_KP_5;
    KP6 = SDL2.SDLK_KP_6;
    KP7 = SDL2.SDLK_KP_7;
    KP8 = SDL2.SDLK_KP_8;
    KP9 = SDL2.SDLK_KP_9;
    KP_PERIOD = SDL2.SDLK_KP_PERIOD;
    KP_DIVIDE = SDL2.SDLK_KP_DIVIDE;
    KP_MULTIPLY = SDL2.SDLK_KP_MULTIPLY;
    KP_MINUS = SDL2.SDLK_KP_MINUS;
    KP_PLUS = SDL2.SDLK_KP_PLUS;
    KP_ENTER = SDL2.SDLK_KP_ENTER;
    KP_EQUALS = SDL2.SDLK_KP_EQUALS;
    UP = SDL2.SDLK_UP;
    DOWN = SDL2.SDLK_DOWN;
    RIGHT = SDL2.SDLK_RIGHT;
    LEFT = SDL2.SDLK_LEFT;
    INSERT = SDL2.SDLK_INSERT;
    HOME = SDL2.SDLK_HOME;
    END = SDL2.SDLK_END;
    PAGEUP = SDL2.SDLK_PAGEUP;
    PAGEDOWN = SDL2.SDLK_PAGEDOWN;
    F1 = SDL2.SDLK_F1;
    F2 = SDL2.SDLK_F2;
    F3 = SDL2.SDLK_F3;
    F4 = SDL2.SDLK_F4;
    F5 = SDL2.SDLK_F5;
    F6 = SDL2.SDLK_F6;
    F7 = SDL2.SDLK_F7;
    F8 = SDL2.SDLK_F8;
    F9 = SDL2.SDLK_F9;
    F10 = SDL2.SDLK_F10;
    F11 = SDL2.SDLK_F11;
    F12 = SDL2.SDLK_F12;
    F13 = SDL2.SDLK_F13;
    F14 = SDL2.SDLK_F14;
    F15 = SDL2.SDLK_F15;
    NUMLOCK = SDL2.SDLK_NUMLOCKCLEAR;
    CAPSLOCK = SDL2.SDLK_CAPSLOCK;
    SCROLLOCK = SDL2.SDLK_SCROLLLOCK;
    RSHIFT = SDL2.SDLK_RSHIFT;
    LSHIFT = SDL2.SDLK_LSHIFT;
    RCTRL = SDL2.SDLK_RCTRL;
    LCTRL = SDL2.SDLK_LCTRL;
    RALT = SDL2.SDLK_RALT;
    LALT = SDL2.SDLK_LALT;
    # RMETA = SDL2.SDLK_RGUI;
    # LMETA = SDL2.SDLK_LGUI;
    LSUPER = SDL2.SDLK_LGUI;
    RSUPER = SDL2.SDLK_RGUI;
    MODE = SDL2.SDLK_MODE;
    HELP = SDL2.SDLK_HELP;
    SYSREQ = SDL2.SDLK_SYSREQ;
    MENU = SDL2.SDLK_MENU;
    POWER = SDL2.SDLK_POWER;
    EURO = SDL2.SDLK_CURRENCYUNIT;
end
end

module Keymods
export Keymod
@enum Keymod::UInt16 begin
    NONE = 0
    LSHIFT = 1
    RSHIFT = 2
    LCTRL = 64
    RCTRL = 128
    LALT = 256
    RALT = 512
    LMETA = 1024
    RMETA = 2048
    NUM = 4096
    CAPS = 8192
    MODE = 16384
    CTRL = 192
    SHIFT = 3
    ALT = 768
    META = 3072
end

end

using .MouseButtons
using .Keys
using .Keymods

struct Keyboard
    pressed::Array{UInt32, 1}
end

Keyboard() = Keyboard(Array{Any, 1}())

function Base.getproperty(k::Keyboard, s::Symbol)
    s = Symbol(uppercase(string(s)))
    if isdefined(GameZero.Keys, s)
        return UInt32(getfield(GameZero.Keys, s)) in getfield(k, :pressed)
    end
    if isdefined(GameZero.Keymods, s)
        return UInt32(getfield(GameZero.Keymods, s)) in getfield(k, :pressed)
    end
    return false
end

Base.push!(k::Keyboard, item) = push!(getfield(k, :pressed), item)
function Base.delete!(k::Keyboard, item)
    a = getfield(k, :pressed)
    deleteat!(a, a.==item)
end


text:
@testset "Noise Extraction Test" begin
  using Random, DiffEqNoiseProcess, DiffEqBase, Test
  seed = 100
  Random.seed!(seed)

  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = (tstart, tend)

  Random.seed!(seed)
  W = WienerProcess(0.0,0.0,0.0)
  prob = NoiseProblem(W,trange)
  sol = solve(prob, dt=dt)

  _sol = deepcopy(sol)
  _sol.save_everystep = false
  sol.save_everystep = false
  tarray = [tstart]
  for i = 1:1000
    t = tarray[i]
    sol(t)
    push!(tarray,t+dt)
  end
  @test length(_sol) == length(sol.W) == length(tarray)
  @test _sol == sol

  Random.seed!(seed)
  W2 = WienerProcess!(0.0,[0.0],[0.0])
  prob2 = NoiseProblem(W2,trange)
  sol2 = solve(prob2, dt=dt)

  sol2.save_everystep = false
  for i in 1:1:length(tarray)
    t = tarray[i]
    sol2(t)
  end

  @test length(sol2.W) == length(tarray)
  @test minimum(_sol .== sol2)

  W3 = NoiseGrid(reverse(_sol.t),reverse(_sol.W))
  prob3 = NoiseProblem(W3,reverse(trange))
  sol3 = solve(prob3, dt=-dt)

  @test minimum(reverse(sol3.W) .== _sol)

  W4 = NoiseWrapper(_sol, reverse=true)
  for i in 1:1:length(tarray)
    t = tarray[end-i+1]
    tmp = DiffEqNoiseProcess.interpolate!(W4,nothing,nothing,t)
    @test _sol.W[end-i+1] ≈ tmp[1]
    @test _sol.Z[end-i+1] ≈ tmp[2]
  end

end


text:
"""
Copyright Marcus Greiff 2020

Tests for the matrix_double_symmetric_real_eigenvalues() function
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture(m)
    A = rand(m,m)
    A = A + A'
    Ac = deepcopy(A)
    Am = matrix_double_t(A)
    B  = zeros(m,1);
    Bc = deepcopy(B)
    Bm = matrix_double_t(B)
    return A, Ac, Am, B, Bc, Bm
end

# Compact call for the matrix_double_symmetric_real_eigenvalues function
function test_call_matrix_double_symmetric_real_eigenvalues(Am, Bm)
    status = ccall((:matrix_double_symmetric_real_eigenvalues, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{matrix_double_t},Ref{matrix_double_t}),
        Am,
        Bm)
    return status
end

@verbose(1, "Testing matrix PSD solver...")
@testset "Matrix solve PD linear systems" begin
    @testset "Good inputs" begin
        for m = 1:10
            # Fixture
            A, Ac, Am, B, Bc, Bm = fixture(m)
            # Solve system in Julia
            X = eigvals(A);
            # Call
            status = test_call_matrix_double_symmetric_real_eigenvalues(Am, Bm)
            # Test
            @test isapprox(status, 1)
            @test isapprox(A, Ac)
            @test !isapprox(B, Bc)
            @test isapprox(X, B)
        end
    end
end


text:
using TensorOperations

function tensordot(A,B,n)
   @assert size(A,n) == size(B,2)
   T = eltype(A)
   d = length(size(A))
   szC = ntuple(i -> (i==n) ? size(B,1) : size(A,i), d)
   C = Array{T}(undef, szC)
   tensordot!(C,A,B,n)
end

function tensordot!(C,A,B,n)
   d = length(size(A))
   T = eltype(C)
   @assert 1<=n<=d
   idA = ntuple(identity, d)
   idB = (0, n)
   idC = ntuple(i -> (i==n) ? 0 : i, d)
   tensorcontract!(one(T), A, idA, 'N', B, idB, 'N', zero(T), C, idC)
end


text:
using BenchmarkTools
using LinearAlgebra
using Test
using TrajectoryOptimization: to_array
using Logging

println("\nRUNNING BENCHMARKS\n")
system = "Simple Pendulum"

# Set up models and objective
u_bound = 2.
model,obj = TrajectoryOptimization.Dynamics.pendulum!
opts = TrajectoryOptimization.SolverOptions()
integration = :rk3

obj.Q .= Diagonal(I,2)*1e-3
obj.R .= Diagonal(I,1)*1e-3
obj.Qf .= Diagonal(I,2)*30
obj.tf = 5.
obj_c = TrajectoryOptimization.ConstrainedObjective(obj, u_min=-u_bound, u_max=u_bound) # constrained objective
obj = TrajectoryOptimization.to_static(obj)
obj_c =TrajectoryOptimization.to_static(obj_c)
dt = 0.1

# Unconstrained
solver = TrajectoryOptimization.Solver(model,obj,dt=dt,opts=opts,integration=integration)
solver.opts.verbose = false
U = ones(model.m,solver.N)
solver.opts.use_static = true
results, = TrajectoryOptimization.solve(solver,U)

err = norm(results.X[end]-obj.xf)
@test err < 1e-3
println("$system - Unconstrained")
disable_logging(Logging.Info)
@btime TrajectoryOptimization.solve(solver,U)

# Constrained
solver = TrajectoryOptimization.Solver(model,obj_c,dt=dt,opts=opts,integration=integration)
solver.opts.verbose = false
solver.opts.cost_tolerance_intermediate = 1e-2
solver.opts.use_static = true
results_c,stats = TrajectoryOptimization.solve(solver,U)

max_c = TrajectoryOptimization.max_violation(results_c)
err = norm(results_c.X[end]-obj.xf)

@test norm(results_c.X[end]-obj.xf) < 1e-3
@test max_c < 1e-2
println("$system - Constrained")
@btime TrajectoryOptimization.solve(solver,U)

### Infeasible Start
obj_c2 = TrajectoryOptimization.update_objective(obj_c)
solver = TrajectoryOptimization.Solver(model, obj_c2, dt=dt, opts=opts, integration=integration)
X_interp = TrajectoryOptimization.line_trajectory(obj.x0, obj.xf,solver.N)
solver.opts.verbose = false
solver.opts.solve_feasible = false
solver.opts.use_static = true
results_inf, stats = TrajectoryOptimization.solve(solver,X_interp,U)

max_c = TrajectoryOptimization.max_violation(results_inf)
err = norm(results_inf.X[end]-obj.xf)

@test norm(results_inf.X[end]-obj.xf) < 1e-3
@test max_c < 1e-2
println("$system - Infeasible Start")
@btime TrajectoryOptimization.solve(solver,U)


text:
function _taskregistry(; name = "Learning tasks")
    registry = Registry(
        (;
            id = Field(String, name = "ID", formatfn = x -> sprint(show, x)),
            name = Field(
                String,
                name = "Name",
                description = "The name of the learning task",
                computefn = (row, key) -> get(row, key, row.id)),
            blocks = Field(
                Any,
                name = "Block types",
                description = "Types of the blocks that are compatible with this task",
                filterfn = blocktypesmatch,
                formatfn = b -> FeatureRegistries.code_format(_formatblock(b))),
            category = Field(
                String,
                name = "Category",
                description = "Kind of task, e.g. \"supervised\""),
            description = Field(
                String,
                name = "Description",
                optional = true,
                description = "More information about the learning task",
                formatfn=FeatureRegistries.md_format),
            constructor = Field(
                Any,
                name = "Learning task",
                description = "Function instance to create a corresponding learning task.",
                formatfn = FeatureRegistries.code_format),
            package = Field(
                Module,
                name = "Package",
                formatfn = FeatureRegistries.code_format),
        );
        name,
        loadfn = row -> row.constructor,
        description = """
        A registry for learning tasks. `load`ing an entry will return a function
        that can be used to construct a `LearningTask` given `blocks`.

        ```julia
        taskfn = load(learningtasks(id))
        task = taskfn(blocks; kwargs...)
        ```

        Inspect `?taskfn` for documentation on the arguments that the function accepts.

        See `datarecipes` to load these datasets in a format compatible with learning
        tasks.
        """
    )
end

const TASKS = _taskregistry()

"""
    learningtasks(; filters...)

Show a registry of available learning tasks. Pass in filters as keyword
arguments to look at a subset.

See also [finding functionality](/documents/docs/discovery.md), [`datasets`](#),
and [`datarecipes`](#). For more information on registries, see
[FeatureRegistries.jl](https://github.com/lorenzoh/FeatureRegistries.jl).

## Examples

Show all available learning tasks:

{cell}
```julia
using FastAI
learningtasks()
```

Show all computer vision tasks:

{cell}
```julia
learningtasks(package=FastAI.Vision)
```

Show all classification tasks, i.e. where the target block is a [`Label`](#):

{cell}
```julia
learningtasks(blocks=(Any, Label))
```

Get an explanation of fields in the learning task registry:

{cell}
```julia
info(learningtasks())
```
"""
learningtasks(; kwargs...) = isempty(kwargs) ? TASKS : filter(TASKS; kwargs...)


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")
include("CESMReader.jl")

using .CESMReader

function mreplace(x)
    return replace(x, missing=>NaN)
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Ocn domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global area = mreplace(ds["area"][:])
    global mask = mreplace(ds["mask"][:])
    global lat  = mreplace(ds["yc"][:])

    area .*= 4π * Re^2 / sum(area)

    global idx_GLB = ( mask .== 1.0 )
    global idx_NH  = ( mask .== 1.0 ) .& ( lat .>= 0.0 )
    global idx_SH  = ( mask .== 1.0 ) .& ( lat .<  0.0 )

end

function calAreaWeightedSum(v, idx)
    return sum( area[idx] .* v[idx] )
end

function calAreaWeightedSumTimeseries(v, idx)
    return [ calAreaWeightedSum( view(v, :, :, t), idx ) for t=1:size(v)[3] ]
end

let
    global ice_volume_GLB, ice_volume_NH, ice_volume_SH
    global ice_area_GLB, ice_area_NH, ice_area_SH

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    vice, aice = getData(fh, ["vice", "aice"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    ice_volume_GLB = calAreaWeightedSumTimeseries(vice, idx_GLB)
    ice_volume_NH  = calAreaWeightedSumTimeseries(vice, idx_NH)
    ice_volume_SH  = calAreaWeightedSumTimeseries(vice, idx_SH)

    ice_area_GLB = calAreaWeightedSumTimeseries(aice, idx_GLB)
    ice_area_NH  = calAreaWeightedSumTimeseries(aice, idx_NH)
    ice_area_SH  = calAreaWeightedSumTimeseries(aice, idx_SH)

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("ice_volume_GLB", ice_volume_GLB, ("time",), Dict()),
        ("ice_volume_NH",  ice_volume_NH,  ("time",), Dict()),
        ("ice_volume_SH",  ice_volume_SH,  ("time",), Dict()),
        ("ice_area_GLB",   ice_area_GLB,   ("time",), Dict()),
        ("ice_area_NH",    ice_area_NH,    ("time",), Dict()),
        ("ice_area_SH",    ice_area_SH,    ("time",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# Calculate the LX76 complexity as described in:
#  F. Kaspar and H.G. Schuster, "Easily Calculable Measure for the
#  Complexity of Spatiotemporal Patterns", Physical Review A, vol 36,
#  num 2, pg 842, 1987.
# This implementation is based on the Matlab code by Stephen Faul:
#  http://www.mathworks.co.uk/matlabcentral/fileexchange/6886-kolmogorov-complexity
#
# This is nice and pure Julia but Zlib.compress is much faster so not much use for this...
function lz76_complexity_kaspar_shuster(s, normalize = true)
  n = length(s)
  c = 1
  l = 1

  i = 0
  k = 1
  k_max = 1
  stop = 0

  while stop==0
    if s[i+k] != s[l+k]
      if k > k_max
        k_max=k
      end
      i=i+1

      if i==l
        c=c+1
        l=l+k_max
        if l+1>n
          stop=1
        else
          i=0
          k=1
          k_max=1
        end
      else
        k=1
      end
    else
      k=k+1
      if l+k > n
        c=c+1
        stop=1
      end
    end
  end

  if normalize
    # a la Lempel and Ziv (IEEE trans inf theory it-22, 75 (1976),
    # h(n)=c(n)/b(n) where c(n) is the kolmogorov complexity
    # and h(n) is a normalised measure of complexity.
    b=n/log2(n)
    return c/b
  else
    return c
  end
end


text:
"""
This tries to add questions to the markdown files

The rendering of the questions depends on the output:

* html: make self-grading questions using JavaScript
* latex: hide the questions
* ipynb: hide the questions
* questions: use idiosyncratic dcf-style markup for questions

"""
nothing

using Mustache, LaTeXStrings
import Base: show

abstract type Question end

MaybeString = Union{String, AbstractString, Nothing}

mutable struct Numericq <: Question
    val::Real
    tol::Real
    reminder
    answer_text::MaybeString
    m::Real
    M::Real
    units
    hint
end

mutable struct Radioq <: Question
    choices::Vector
    answer::Int
    reminder
    answer_text::MaybeString
    values
    labels
    hint
    inline::Bool

end

mutable struct Multiq <: Question
    choices::Vector
    answer::Vector{Int}
    reminder::AbstractString
    answer_text::MaybeString
    values
    labels
    answers
    hint
    inline::Bool
end

mutable struct Shortq <: Question
    answer
    reminder::AbstractString
    answer_text::MaybeString
    hint
end

mutable struct Longq <: Question
    reminder::AbstractString
    answer_text::MaybeString
    hint
    rows::Int
    cols::Int
end

"""
A numeric question graded with a tolerance

Arguments:

* `val::Real` answer

* `tol::Real` tolerance. Answer is right if `|ans - val| <= tol`

* `reminder` a reminder as to what question is, student may see

* `answer_text`: reminder of what answer is, student does not see

* `hint`: a possible hint for the student

* `units`: a string holding the units, if specified.

Returns an object of type `Question`.

Example

```
numericq(10, 1e-3, "what is 5 + 5?", units="An integer")
```
"""
function numericq(val, tol=1e-3, reminder="", args...; hint::AbstractString="", units::AbstractString="")
    answer_text= "[$(round(val-tol,digits=5)), $(round(val+tol,digits=5))]"
    Numericq(val, tol, reminder, answer_text, val-tol, val+tol, units, hint)
end

numericq(val::Int; kwargs...) = numericq(val, 0; kwargs...)

"""
Multiple choice question

Arguments:

* `choices`: vector of choices.

* `answer`: index of correct choice

* `inline::Bool`: hint to render inline (or not) if supported

Example
```
radioq(["beta", L"\beta", "`beta`"], 2, "a reminder", hint="which is the Greek symbol")
```
"""
function radioq(choices, answer, reminder="", answer_text=nothing;  hint::AbstractString="", inline::Bool=(hint!=""),
                keep_order::Bool=false)
    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices # map(markdown_to_latex,choices) |> x -> map(chomp, x) ##|> x -> join(x, " | ")
    !keep_order && shuffle!(inds)

    Radioq(choices[inds], findfirst(isequal(answer), inds), reminder, answer_text, values, labels[inds], hint, inline)
end

"""
True of false questions:

Example:

```
booleanq(true, reminder="Does it hurt...")
```

"""
function booleanq(ans::Bool, reminder="", answer_text=nothing;labels::Vector=["true", "false"], hint::AbstractString="", inline::Bool=true)
    choices = labels[1:2]
    ans = 2 - ans
    radioq(choices, ans, reminder, answer_text; hint=hint, inline=inline, keep_order=true)
end

"""

Boolean question with `yes` or `no` labels.

Examples: `yesnoq("yes")` or `yesnoq(true)`

"""
yesnoq(ans::AbstractString, args...; kwargs...) = radioq(["Yes", "No"], ans == "yes" ? 1 : 2, args...; keep_order=true, kwargs...)
yesnoq(ans::Bool, args...; kwargs...) = yesnoq(ans ? "yes" : "no", args...;kwargs...)

function multiq(choices, answer, reminder="", answer_text=nothing; hint::AbstractString="", inline::Bool=false)
    values = join(1:length, " | ")
    labels =  map(markdown_to_latex, choices) |> x -> map(chomp, x) |> x -> join(x, " | ")
    answers = join(answer, " | ")
    Multiq(choices, answer, reminder, answer_text,
           values, labels, answers,
           hint,
           inline
           )
end

"""

Short question, has regular expression grading:

Example:

```
shortq("x^2", L"a expression using powers to compute x\\cdot x")
```

Do not use format in answer part. The reminder defaults to an empty string.
"""
function shortq(answer, reminder="", answer_text=nothing;hint::AbstractString="")
    Shortq(answer, reminder, answer_text, hint)
end

"""

Long questions, not graded

Example

```
longq("a reminder", "an answer for the grader to see")
```
"""
function longq(reminder="", answer_text=nothing;hint::AbstractString="", rows=3, cols=60)
    Longq(reminder, answer_text, hint, rows, cols)
end

export numericq, radioq, booleanq, yesnoq, shortq, longq, multiq

## we have different display mechanisms based on the output type

## application/x-latexq
latexq_templates=Dict()

latexq_templates["Numericq"] = mt"""
\\begin{answer}
    type: numeric
    reminder: {{{reminder}}}
    answer: [{{m}}, {{M}}]
    {{#answer_text}}answer_text: {{{answer_text}}} {{/answer_text}}
\\end{answer}
"""

latexq_templates["Radioq"] = mt"""
\\begin{answer}
type: radio
reminder: {{{:reminder}}}
values: {{{:values}}}
labels: {{{:labels}}}
answer: {{{:answer}}}
{{#:answer_text}}answer_text: {{{:answer_text}}} {{/:answer_text}}
\\end{answer}
"""

latexq_templates["Multiq"] = mt"""
\\begin{answer}
type: checkbox
reminder: {{{reminder}}}
values: {{{values}}}
labels: {{{labels}}}
answer: {{{answer}}}
{{#answer_text}}answer_text: {{{answer_text}}} {{/answer_text}}
\\end{answer}
"""

latexq_templates["Shortq"] =  mt"""
\\begin{answer}
type: shorttext
reminder: {{{reminder}}}
answer: {{{answer}}}
{{#answer_text}}answer_text: {{{answer_text}}} {{/answer_text}}
\\end{answer}
"""

latexq_templates["Longq"] = mt"""
\\begin{answer}
type: longtext
reminder: {{{reminder}}}
{{#answer_text}}answer_text: {{{answer_text}}} {{/answer_text}}
rows: {{{rows}}}
cols: {{{cols}}}
\\end{answer}
"""

function show(io::IO, m::MIME"application/x-latexq", x::Radioq)
    Mustache.render(io, latexq_templates["Radioq"],
                    Dict(:reminder=>x.reminder,
                         :values=> join(x.values, " | "),
                         :labels=> join(map(WeavePynb.markdown_to_latex, x.labels), " | "),
                         :answer=> x.answer,
                         :answer_text=> x.answer_text
                         )
                    )
end
function show(io::IO, m::MIME"application/x-latexq", x::Question)
    tof = split(string(typeof(x)), ".")[end]
    Mustache.render(io, latexq_templates[tof],x)
end

##################################################
## application/x-latex
## questions for printed output of projects
latex_templates=Dict()

latex_templates["Numericq"] = mt"""
\vspace{12pt}
\noindent\textit{Enter a number:} \underline{\hspace{5cm}}
\vspace{12pt}
"""

latex_templates["Radioq"] = mt"""
\vspace{12pt}
\noindent\textit{Make a selection:}
\begin{enumerate}
{{#labels}}\item {{{.}}} {{/labels}}
\end{enumerate}
\verb+  -------------  +
\vspace{12pt}
"""

latex_templates["Multiq"] = mt"""
\vspace{12pt}
\noindent\textit{Make one or more selections:}
\begin{enumerate}
{{#labels}}\item {{{.}}} {{/labels}}
\end{enumerate}
\verb+  -------------  +
\vspace{12pt}
"""

latex_templates["Shortq"] =  mt"""
\vspace{12pt}
\noindent\textit{Your answer:}
\vspace{0.5in}\\
"""

latex_templates["Longq"] = mt"""
\vspace{12pt}
\textit{Your answer:}
\vspace{2in}\\
"""

function show(io::IO, m::MIME"application/x-latex", x::Question)
    tof = split(string(typeof(x)), ".")[end]
    Mustache.render(io, latex_templates[tof],x)
end

## text/html
html_templates=Dict()

html_templates["Numericq"] = mt"""
<form name='WeaveQuestion' data-id='{{ID}}' data-controltype='{{TYPE}}'>

{{{form}}}
{{#hint}}

{{/hint}}

<input id="{{ID}}" type="number" class="form-control">
{{#units}}{{{units}}}{{/units}}

</form>

"""

function show(io::IO, m::MIME"text/html", x::Numericq)
    d = Dict()
    d["ID"] = randstring()
    d["TYPE"] = "numeric"
    d["selector"] = "#" * d["ID"]
    d["status"] = ""
    d["hint"] = ""# x.hint
    d["units"] = x.units
    d["correct"] = "Math.abs(this.value - $(x.val)) <= $(x.tol)"
    out =  Mustache.render(html_templates["Numericq"], d)
    Mustache.render(io, html_templates["Numericq"], d)
end

html_templates["Radioq"] = mt"""
{{#items}}

<input type='radio' name='radio_{{ID}}' value='{{value}}'>
{{{label}}}

{{/items}}

"""

html_templates["question_tpl"] = mt"""
<form name="WeaveQuestion" data-id="{{ID}}" data-controltype="{{TYPE}}">

{{{form}}}
{{#hint}}

{{/hint}}

</form>

"""

html_templates["script_tpl"] = mt"""
$("{{{selector}}}").on("change", function() {
  correct = {{{correct}}};

  if(correct) {
     $("#{{ID}}_message").html("
 Correct
");
  } else {
     $("#{{ID}}_message").html("
 Incorrect
");
  }
});
"""

function markdown(x)
    length(x) == 0 && return("")

    x = Markdown.parse(x)
    x = sprint(io -> WeavePynb.tohtml(io, x))
#    x[4:end-4]                  # strip out

    x
end

function show(io::IO, m::MIME"text/html", x::Radioq)
    ID = randstring()

    tpl = mt"""
    {{#items}}

      <input type="radio" name="radio_{{ID}}" value="{{value}}">{{{label}}}

    {{/items}}
"""

choices = map(string, x.choices)

    items = Dict[]
    ## make items
    for i in 1:length(choices)
        item = Dict("no"=>i,
                "label"=>markdown(choices[i]),
                "value"=>i
                )
        push!(items, item)
    end

    script = Mustache.render(html_templates["script_tpl"],
                             Dict("ID"=>ID,
                              "selector"=>"input:radio[name='radio_$ID']",
                              "correct"=>"this.value == $(x.answer)"))

    form = Mustache.render(tpl, Dict("ID"=>ID, "items"=>items,
                                 "inline" => x.inline ? " inline" : ""
                                 ))

    Mustache.render(io, html_templates["question_tpl"],
    Dict("form"=>form, "script"=>script,
         "TYPE"=>"radio",
         "ID"=>ID, "hint"=>markdown(x.hint)))

end

function show(io::IO, m::MIME"text/html", x::Question)
     println(io, "
Question type $(typeof(x)) is not supported in this format
")
end

##
## text/markdown
md_templates=Dict()

md_templates["Radioq"] = mt"""
{{#:items}}
{{{.}}}
{{/:items}}

"""


text:
"""
	duration(o)

Return the sum of the lengthes of the presence
intervals of the given object.
"""
duration(o::StreamObject)=length(o.presence)

"""
	duration(s)

Return the duration of the stream, that is tend - tstart.
"""
duration(s::AbstractStream)=length(s.T)

"""
	duration(W)

Return the sum of the nodes' durations in W.
"""
function duration(W::Dict{AbstractString,Node})
	length(W) != 0 ? sum([duration(v) for (k,v) in W]) : 0.0
end

"""
	node_duration(ls)

Return the node duration in the given Link Stream.
By definition, this is just the duration of the link stream itself.
"""
node_duration(ls::Union{LinkStream,DirectedLinkStream})=duration(ls)

"""
	node_duration(s)

Return the node duration in the given Stream Graph.
The node duration is computed as:
```math
k = \\frac{\\left| W \\right|}{\\left| V \\right|}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_duration(s::Union{StreamGraph,DirectedStreamGraph})
	length(s.V)>0 ? duration(s.W)/length(s.V) : 0
end

"""
	node_duration(s,tc)

Return the node duration in the given Stream Graph using a TimeCursor.
In this case, the node duration is computed as:
```math
k = \\int_{k \\in T} k_t dt
```
Where
```math
k_t=\\frac{\\left| V_t \\right|}{\\left| V \\right|}
```
is the node contributions of a time instant t.
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_duration(s::Union{StreamGraph,DirectedStreamGraph},tc::TimeCursor)
	k::Float64=0.0
	length(s.V)==0 && return 0.0
	start!(tc)
	k+=duration(tc.S)*number_of_nodes(tc)
	while tc.S.t1 < s.T.list[1][2]
		next!(tc)
		k+=duration(tc.S)*number_of_nodes(tc)
	end
	k/length(s.V)
end

"""
	link_duration(s)

Return the link duration in the given Stream.
The link duration is computed as:
```math
l = \\frac{\\left| E \\right|}{\\left| V \\otimes V \\right|}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function link_duration(s::AbstractStream)
	length(s.V)>1 ? 2*sum([duration(l) for (k,v) in s.E for (kk,l) in v])/(length(s.V)*(length(s.V)-1)) : 0
end

"""
	link_duration(s,tc)

Return the link duration in the given Stream using a TimeCursor.
In this case, the link duration is computed as:
```math
l = \\int_{t \\in T} l_t dt
```
Where
```math
l_t = \\frac{ \\left| E_t \\right|}{\\left| V \\otimes V \\right|}
```
is the link contributions of a time instant t.
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function link_duration(s::AbstractStream,tc::TimeCursor)
	k::Float64=0.0
	length(s.V)<=1 && return 0.0
	start!(tc)
	k+=duration(tc.S)*number_of_links(tc)
	while tc.S.t1 < s.T.list[1][2]
		next!(tc)
		k+=duration(tc.S)*number_of_links(tc)
	end
	k/(length(s.V)*(length(s.V)-1))
end

"""
	duration(s)

Return the duration of a given state.
"""
duration(s::State)=s.t1-s.t0

"""
	duration(j)

Return the duration of a given DurationJump.
"""
duration(j::DurationJump)=j.δ

"""
	duration(p)

Return the duration of a given Path.
"""
duration(p::AbstractPath)=finish(p)-start(p)


text:
"""
The structure used for user requests.
"""
mutable struct UserRequest
  id::String
  name::String
  email::String
  address::String
  organization::String
  licenseType::String
  billingId::String
end

"""
The response structure for user calls.
"""
mutable struct UserResponse
  userId::String
  email::String
  address::String
  organization::String
  licenseType::String
  billingId::String
  id::String
  name::String
  createdTimestamp::String
  lastUpdatedTimestamp::String
  links::Dict{String, Any}
end


text:
using Flux; Flux.CuArrays.allowscalar(false)
using Flux: params, throttle, Optimise.@epochs
using SparseGaussianProcesses
import Plots

gp = SparseGaussianProcess(SquaredExponentialKernel(1), inducing_points = PseudoDataInducingPoints)

f(x) = 2 * sin.(x |> vec) + randn(length(x))/10
x = reshape(-5:0.1:5, (1,:)) |> collect
y = f(x)' |> collect
gp.inducing_points(reshape(-5:1:5,1,:), gp.kernel)

(gp,x,y) = gpu.((gp,x,y))

rand!(gp.prior_basis, gp.kernel; num_features=67)
rand!(gp; num_samples=17)
loss(gp,x,y)|>sum

opt = ADAM(0.01)
dataset = Iterators.repeated((x,y), 200)

Flux.train!((x,y) -> loss(gp,x,y)|>sum, params(gp), dataset, opt; cb = throttle(() -> @show(loss(gp,x,y)|>sum), 0.1))

(gp,x,y) = cpu.((gp,x,y))

plot_gp_intervals(gp, x, y)


text:
__precompile__()

module BlockArrays

# AbstractBlockArray interface exports
export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
export Block, getblock, getblock!, setblock!, nblocks, blocksize, blockcheckbounds, BlockBoundsError

export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat
export PseudoBlockArray, PseudoBlockMatrix, PseudoBlockVector, PseudoBlockVecOrMat

import Base: @propagate_inbounds, full
using Base.Cartesian

include("abstractblockarray.jl")
include("blocksizes.jl")
include("blockindices.jl")
include("blockarray.jl")
include("pseudo_blockarray.jl")
include("operations.jl")
include("utilities.jl")
include("show.jl")

end # module


text:
ismacro(x :: Expr) = Meta.isexpr(x, :macrocall)
ismacro(_) = false

function flatten_macros(node :: Expr)
    @match node begin
    Expr(:macrocall, op :: Symbol, ::LineNumberNode, arg) ||
    Expr(:macrocall, op :: Symbol, arg) =>

    @match arg begin
    Expr(:tuple, args...) || a && Do(args = [a]) =>

    @match args begin
    [args..., function ismacro end && tl] => [(op |> get_op, args), flatten_macros(tl)...]
    _ => [(op |> get_op, args)]

    end
    end
    end
end

function codegen(node)
    ops = flatten_macros(node)
    let rec(vec) =
        @match vec begin
            [] => []
            [(&generate_groupby, args1), (&generate_having, args2), tl...] =>
                [generate_groupby(args1, args2), rec(tl)...]
            [(hd, args), tl...] =>
                [hd(args), rec(tl)...]
        end
        init = quote
            let iter = $get_records($ARG),
                fields = $get_fields($ARG),
                types =$type_unpack($length(fields), $eltype(iter))
                (fields, types, iter)
            end
        end
        fn_body = foldl(rec(ops), init = init) do last, mk
            mk(last)
        end
        quote
            @inline function ($ARG :: $TYPE_ROOT, ) where {$TYPE_ROOT}
                let ($IN_FIELDS, $IN_TYPES, $IN_SOURCE) = $fn_body
                    $build_result(
                        $TYPE_ROOT,
                        $IN_FIELDS,
                        $IN_TYPES,
                        $IN_SOURCE
                    )
                end
            end
        end
    end
end

function generate_select
end

function generate_where
end

function generate_groupby
end

function generate_orderby
end

function generate_having
end

function generate_limit
end

const registered_ops = Dict{Symbol, Any}(
    Symbol("@select") => generate_select,
    Symbol("@where") => generate_where,
    Symbol("@groupby") => generate_groupby,
    Symbol("@having") => generate_having,
    Symbol("@limit") => generate_limit,
    Symbol("@orderby") => generate_orderby
)

function get_op(op_name)
    registered_ops[op_name]
end

macro select(node)
    codegen(Expr(:macrocall, Symbol("@select"), __source__, node)) |> esc
end

macro where(node)
    codegen(Expr(:macrocall, Symbol("@where"), __source__, node)) |> esc
end

macro groupby(node)
    codegen(Expr(:macrocall, Symbol("@groupby"), __source__, node)) |> esc
end

macro having(node)
    codegen(Expr(:macrocall, Symbol("@having"), __source__, node)) |> esc
end

macro orderby(node)
    codegen(Expr(:macrocall, Symbol("@orderby"), __source__, node)) |> esc
end

macro limit(node)
    codegen(Expr(:macrocall, Symbol("@limit"), __source__, node)) |> esc
end


text:
function Base.iterate(layer::AbstractFeatureLayer, state::Int=0)
    layer.ptr == C_NULL && return nothing
    state == 0 && resetreading!(layer)
    ptr = GDAL.ogr_l_getnextfeature(layer.ptr)
    if ptr == C_NULL
        resetreading!(layer)
        return nothing
    else
        return (Feature(ptr), state+1)
    end
end

Base.eltype(layer::AbstractFeatureLayer) = Feature

Base.length(layer::AbstractFeatureLayer) = nfeature(layer, true)

struct BlockIterator
    rows::Cint
    cols::Cint
    ni::Cint
    nj::Cint
    n::Cint
    xbsize::Cint
    ybsize::Cint
end

function blocks(raster::AbstractRasterBand)
    (xbsize, ybsize) = blocksize(raster)
    rows = height(raster)
    cols = width(raster)
    ni = ceil(Cint, rows / ybsize)
    nj = ceil(Cint, cols / xbsize)
    BlockIterator(rows, cols, ni, nj, ni * nj, xbsize, ybsize)
end

function Base.iterate(obj::BlockIterator, iter::Int=0)
    iter == obj.n && return nothing
    j = floor(Int, iter / obj.ni)
    i = iter % obj.ni
    nrows = if (i + 1) * obj.ybsize < obj.rows
        obj.ybsize
    else
        obj.rows - i * obj.ybsize
    end
    ncols = if (j + 1) * obj.xbsize < obj.cols
        obj.xbsize
    else
        obj.cols - j * obj.xbsize
    end
    (((i, j), (nrows, ncols)), iter+1)
end

struct WindowIterator
    blockiter::BlockIterator
end
Base.size(i::WindowIterator) = (i.blockiter.ni, i.blockiter.nj)
Base.length(i::WindowIterator) = i.blockiter.n
Base.IteratorSize(::Type{WindowIterator}) = Base.HasShape{2}()
Base.IteratorEltype(::Type{WindowIterator}) = Base.HasEltype()
Base.eltype(::WindowIterator) = Tuple{UnitRange{Int}, UnitRange{Int}}

windows(raster::AbstractRasterBand) = WindowIterator(blocks(raster))

function Base.iterate(obj::WindowIterator, iter::Int=0)
    handle = obj.blockiter
    next = Base.iterate(handle, iter)
    next == nothing && return nothing
    (((i, j), (nrows, ncols)), iter) = next
    (((1:ncols) .+ j * handle.xbsize, (1:nrows) .+ i * handle.ybsize), iter)
end

mutable struct BufferIterator{T <: Real}
    raster::AbstractRasterBand
    w::WindowIterator
    buffer::Array{T, 2}
end

function bufferwindows(raster::AbstractRasterBand)
    BufferIterator(
        raster,
        windows(raster),
        Array{pixeltype(raster)}(undef, blocksize(raster)...)
    )
end

function Base.iterate(obj::BufferIterator, iter::Int=0)
    next = Base.iterate(obj.w, iter)
    next == nothing && return nothing
    ((cols, rows), iter) = next
    rasterio!(obj.raster, obj.buffer, rows, cols)
    return (obj.buffer[1:length(cols), 1:length(rows)], iter)
end


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

mutable struct EncryptionService <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    lastEnabledTime::Any # spec type: Union{ Nothing, DateTime } # spec name: lastEnabledTime

    function EncryptionService(;enabled=nothing, lastEnabledTime=nothing)
        o = new()
        validate_property(EncryptionService, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(EncryptionService, Symbol("lastEnabledTime"), lastEnabledTime)
        setfield!(o, Symbol("lastEnabledTime"), lastEnabledTime)
        o
    end
end # type EncryptionService

const _property_map_EncryptionService = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("lastEnabledTime")=>Symbol("lastEnabledTime"))
const _property_types_EncryptionService = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("lastEnabledTime")=>"DateTime")
Base.propertynames(::Type{ EncryptionService }) = collect(keys(_property_map_EncryptionService))
Swagger.property_type(::Type{ EncryptionService }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_EncryptionService[name]))}
Swagger.field_name(::Type{ EncryptionService }, property_name::Symbol) =  _property_map_EncryptionService[property_name]

function check_required(o::EncryptionService)
    true
end

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


text:
using ArrayFire
using AutoGrad
importall Base
using AutoGrad: matmul2arg,broadcast2arg,unbroadcast,id2,math2arg,float2arg,broadcast2cmp
import AutoGrad: tofloat,_dbg

tofloat{T<:AbstractFloat}(a::AFArray{T})=a
_dbg{T}(x::AFArray{T})=Symbol("AF$(join([id2(x),size(x)...,8*sizeof(T)],'_'))")
typealias AForN Union{AFArray,Number}
for (f,d) in matmul2arg
    @eval @primitive $f(x1::AFArray,x2::AFArray)::y $(d[1]) $(d[2])
end
for (f,g) in broadcast2arg
    @eval @primitive $f(x1::AFArray,x2::AFArray)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
    @eval @primitive $f(x1::Number,x2::AFArray)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
    @eval @primitive $f(x1::AFArray,x2::Number)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
end
for (f,g) in math2arg
    @eval @primitive $f(x1::AFArray,x2::AFArray)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
    @eval @primitive $f(x1::Number,x2::AFArray)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
    @eval @primitive $f(x1::AFArray,x2::Number)::y  unbroadcast(y,x1,dy->dy.*$(g[1]))  unbroadcast(y,x2,dy->dy.*$(g[2]))
end
for (f,g) in float2arg
    @eval @primitive $f(x1::AFArray,x2::AFArray)::y unbroadcast(y,x1,$(g[1])) unbroadcast(y,x2,$(g[2]))
    @eval @primitive $f(x1::Number,x2::AFArray)::y unbroadcast(y,x1,$(g[1])) unbroadcast(y,x2,$(g[2]))
    @eval @primitive $f(x1::AFArray,x2::Number)::y unbroadcast(y,x1,$(g[1])) unbroadcast(y,x2,$(g[2]))
end
for f in broadcast2cmp
    @eval begin
        @zerograd $f(x1::AFArray,x2::AFArray)
        @zerograd $f(x1::AFArray,x2::Number)
        @zerograd $f(x1::Number,x2::AFArray)
    end
end

Base.eltype{T}(x::AFArray{T})=T
@primitive  sum(x::AFArray)  (dy->convert(eltype(x),dy).+zeros(x))
@primitive  sum(x::AFArray,i...)  (dy->dy.+zeros(x))
Base.size(x::AFArray,i)=(if i>ndims(x); 1; else; size(x)[i]; end)

include("profile.jl")
d0af = map(i->map(AFArray,i),d0)
# w1af = map(AFArray,w1)
w2af = map(AFArray,w2)
# timeall(w2af,d0af,10)
# loop(fun[9],w2af,d0af,10)
# (x1af,y1af)=first(d0af)
# fun[10](w2af,x1af,y1af)

function timeall_af(w=w2af,d=d0af,t=10)
    for i=1:length(fun)
        print(i); printfun(fun[i])
        for j=1:2
            sleep(2)
            gc_enable(false)
            @time (loop_af(fun[i],w,d,t); sync())
            gc_enable(true)
            @show deviceMemInfo()
        end
    end
end

function loop_af(f,w,d,t)
    for i in 1:t
        for (x,y) in d
            # ArrayFire.af_device_gc()
            f(w,x,y)
        end
    end
end

gc_enable(false); @time loop_af(fun[9],w2af,d0af,10); gc_enable(true)
gc_enable(false); @time loop_af(fun[9],w2af,d0af,10); gc_enable(true)
gc_enable(false); @time loop_af(fun[9],w2af,d0af,10); gc_enable(true)
gc_enable(false); @time loop_af(fun[9],w2af,d0af,10); gc_enable(true)
gc_enable(false); @time loop_af(fun[9],w2af,d0af,10); gc_enable(true)
gc_enable(false); @profile loop_af(fun[9],w2af,d0af,10); gc_enable(true)
Profile.print()
Profile.print(format=:flat)


text:
#!/usr/bin/env julia
using Bio.Seq
import JSON
using DataStructures

function readfile_sampleids(filename::String)
    ctr = counter(String)
    open(FASTQReader, filename, quality_encoding=:sanger) do reader
        rec = FASTQSeqRecord{DNASequence}()
        while !eof(reader)
            read!(reader, rec)
            if endswith(rec.name, "/2")
                # Skip R2, otherwise counts are inflated
                continue
            end
            rdata = JSON.parse(rec.metadata.description)
            sample = rdata["id"]
            push!(ctr, sample)
        end
    end
    return ctr
end

type Assessment
    total::Int
    correct::Int
    incorrect::Int
    unassigned::Int
    Assessment() = new(0,0,0,0)
    Assessment(t, c, i, u) = new(t, c, i, u)
end

function assess_samples(truthfile, sampledir)
    # Read true counts
    samplecounts = JSON.parsefile(truthfile)

    # Find unknown file and read it
    unknowncounts = readfile_sampleids(joinpath(sampledir, "unknown.fastq"))

    sample_stats = Dict{String, Assessment}()
    for sample in keys(samplecounts)
        fname = joinpath(sampledir, "$(sample).fastq")
        got = readfile_sampleids(fname)
        nread_expect = samplecounts[sample]
        nread_got_sample = got[sample]
        nread_unassigned = unknowncounts[sample]
        # Wrong is those not unassigned, and not in correct file. Must be
        # elsewhere
        nread_got_wrong = nread_expect - nread_unassigned -  nread_got_sample
        a = Assessment(nread_expect, nread_got_sample, nread_got_wrong,
                       nread_unassigned)
        sample_stats[sample] = a
    end
    return sample_stats
end

function append!(a::Assessment, b::Assessment)
    a.total += b.total
    a.correct += b.correct
    a.incorrect += b.incorrect
    a.unassigned += b.unassigned
end

function summarise_assessment(assessment)
    total = Assessment()
    for (s, a) in assessment
        append!(total, a)
    end
    return total
end

function dump_assessment(file, assess)
    open(file, "w") do f
        println(f, "Sample\tTotal\tCorrect\tIncorrect\tUnassigned")
        for (s, a) in assess
            println(f, "$s\t$(a.total)\t$(a.correct)\t$(a.incorrect)\t$(a.unassigned)")
        end
    end
end

##################
#  END OF FUNCS  #
##################

const afile = ARGS[1]
const seed = ARGS[2]
const barcode_set = ARGS[3]
const demuxer = ARGS[4]
const truthfile = ARGS[5]
const sampledir = ARGS[6]

const assessment = assess_samples(truthfile, sampledir)
dump_assessment(afile, assessment)

const tot = summarise_assessment(assessment)

println("$seed\t$demuxer\t$barcode_set\t$(tot.correct)\t$(tot.incorrect)\t$(tot.unassigned)")


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    jdcnv(date) -> julian_days

### Purpose ###

Convert proleptic Gregorian Calendar date in UTC standard to number of Julian
days.

### Explanation ###

Takes the given proleptic Gregorian date in UTC standard and returns the number
of Julian calendar days since epoch `-4713-11-24T12:00:00`.

### Argument ###

* `date`: date in proleptic Gregorian Calendar.  Can be either a single date or
  an array of dates.  Each element can be either a `DateTime` type or anything
  that can be converted directly to `DateTime`.  In the case of vectorial input,
  each element is considered as a date, so you cannot provide a date by parts.

### Output ###

Number of Julian days, as a floating point.

### Example ###

Find the Julian days number at 2016 August 23, 03:39:06.

``` julia
jdcnv(DateTime(2016, 08, 23, 03, 39, 06))
# => 2.4576236521527776e6
jdcnv(2016, 08, 23, 03, 39, 06)
# => 2.4576236521527776e6
jdcnv("2016-08-23T03:39:06")
# => 2.4576236521527776e6
```

### Notes ###

This is the inverse of `daycnv`.

`get_juldate` returns the number of Julian days for current time.  It is
equivalent to `jdcnv(now(Dates.UTC))`.

For the conversion of Julian date to number of Julian days, use `juldate`.
"""
const jdcnv = Dates.datetime2julian
jdcnv(date...) = jdcnv(DateTime(date...))

function jdcnv{T<:Any}(dt::AbstractArray{T})
    jd = similar(dt, Float64)
    for i in eachindex(dt)
        jd[i] = jdcnv(dt[i])
    end
    return jd
end


text:
module MultipleDispatchKoans

struct Dog
    name::String
end

struct Cat
    name::String
end

function hello()
end

function hello(name)
end

function pet_talk(pet::Dog)
end

function pet_talk(pet::Cat)
end

end


text:
# Determine the local system's time zone
# Based upon Python's tzlocal https://pypi.python.org/pypi/tzlocal
using Mocking: Mocking, @mock

if Sys.iswindows()
    import TimeZones.WindowsTimeZoneIDs: WINDOWS_TRANSLATION
end

"""
    localzone() -> TimeZone

Returns a `TimeZone` object that is equivalent to the system's current time zone.
"""
function localzone()
    # Only allow creating a TimeZone using standard and legacy IANA time zone database
    # names. We allow the use of legacy names here as most operating systems still use the
    # legacy names.
    mask = Class(:STANDARD) | Class(:LEGACY)

    @static if Sys.isunix()
        name = ""

        # Try getting the time zone from the "TZ" environment variable
        # https://linux.die.net/man/3/tzset
        #
        # Note: The macOS man page tzset(3) states some additional information about the
        # colon being optional:
        #
        # > If TZ appears in the environment and its value begins with a colon (`:'), the
        # > rest of its value is used as a pathname of a tzfile(5)-format file from which to
        # > read the time conversion information. If the first character of the pathname is
        # > a slash (`/'), it is used as an absolute pathname; otherwise, it is used as a
        # > pathname relative to the system time conversion information directory.
        # >
        # > If its value does not begin with a colon, it is first used as the pathname of a
        # > file (as described above) from which to read the time conversion information. If
        # > that file cannot be read, the value is then interpreted as a direct
        # > specification (the format is described below) of the time conversion
        # > information.
        if haskey(ENV, "TZ")
            name = ENV["TZ"]
            force_tzfile = false

            # If the TZ format starts with a colon we'll prefer using the time zone
            # information from the specified file.
            if startswith(name, ':')
                name = name[2:end]  # Either a relative or absolute path
                force_tzfile = true
            else
                # Check if name matches pre-compiled time zone
                istimezone(name, mask) && return TimeZone(name, mask)
            end

            if startswith(name, '/')
                return @mock open(name) do f
                    read_tzfile(f, something(_path_tz_name(name, mask), "local"))
                end
            else
                # The system time zone directory used depends on the (g)libc version
                tzdirs = ["/usr/lib/zoneinfo", "/usr/share/zoneinfo"]
                haskey(ENV, "TZDIR") && pushfirst!(tzdirs, ENV["TZDIR"])

                for dir in tzdirs
                    filepath = joinpath(dir, name)
                    (@mock isfile(filepath)) || continue
                    return @mock open(filepath) do f
                        read_tzfile(f, name)
                    end
                end

                if !force_tzfile
                    return parse_tz_format(name)
                else
                    error("Unable to locate tzfile: $name")
                end
            end
        end

        # Look for Linux distribution configuration files that contain the time zone name.
        # Since we don't expect these files on macOS we'll avoid doing these checks to
        # improve performance.
        @static if Sys.islinux()
            filename = "/etc/timezone"
            if @mock isfile(filename)
                @mock open(filename) do file
                    name = read(file, String)

                    # Get rid of host definitions and comments:
                    name = strip(replace(name, r"#.*" => ""))
                    name = replace(name, ' ' => '_')
                end

                istimezone(name, mask) && return TimeZone(name, mask)
            end

            # CentOS has a ZONE setting in /etc/sysconfig/clock,
            # OpenSUSE has a TIMEZONE setting in /etc/sysconfig/clock and
            # Gentoo has a TIMEZONE setting in /etc/conf.d/clock

            zone_re = r"(TIME)?ZONE\s*=\s*\"(?<name>.*?)\""
            for filepath in ("/etc/sysconfig/clock", "/etc/conf.d/clock")
                (@mock isfile(filepath)) || continue
                @mock open(filepath) do file
                    for line in readlines(file)
                        matched = match(zone_re, line)
                        if matched != nothing
                            name = matched["name"]
                            name = replace(name, ' ' => '_')
                            break
                        end
                    end
                end

                istimezone(name, mask) && return TimeZone(name, mask)
            end
        end

        # systemd distributions use symlinks that include the zone name,
        # see man page of localtime(5) and timedatectl(1)
        link = "/etc/localtime"
        if @mock islink(link)
            target = @mock readlink(link)
            name = _path_tz_name(target, mask)
            name !== nothing && return TimeZone(name, mask)
        end

        # No explicit setting existed. Use localtime
        for filepath in ("/etc/localtime", "/usr/local/etc/localtime")
            (@mock isfile(filepath)) || continue
            return @mock open(filepath) do f
                read_tzfile(f, "local")
            end
        end
    elseif Sys.iswindows()
        # The Windows Time Zone Utility (tzutil) is pre-installed from Windows XP and later.
        # This approach was used as it is the fastest without without adding additional
        # package dependencies.
        #
        # https://docs.microsoft.com/en-us/windows-server/administration/windows-commands/tzutil
        # https://docs.microsoft.com/en-us/previous-versions/windows/it-pro/windows-server-2012-r2-and-2012/hh875624(v=ws.11)
        #
        # Alternative approaches:
        # - Call .NET via Powershell: `powershell -Command "[TimeZoneInfo]::Local.Id"`
        # - Read the Windows registry
        win_name = @mock read(`tzutil /g`, String)

        if haskey(WINDOWS_TRANSLATION, win_name)
            return TimeZone(WINDOWS_TRANSLATION[win_name], mask)
        else
            error("unable to translate to POSIX time zone name from: \"$win_name\"")
        end
    end

    error("Failed to find local time zone")
end

# Extract a time zone name from a path.
# e.g. "/usr/share/zoneinfo/Europe/Warsaw" becomes "Europe/Warsaw"
function _path_tz_name(path::AbstractString, mask::Class=Class(:ALL))
    i = 0
    while i !== nothing
        name = SubString(path, i + 1)
        istimezone(name, mask) && return name
        i = findnext(isequal('/'), path, i + 1)
    end

    return nothing
end

"""
    parse_tz_format(str) -> TimeZone

Parse the time zone format typically provided via the "TZ" environment variable. Details on
the format can be found under the man page for
[tzset](https://man7.org/linux/man-pages/man3/tzset.3.html).
"""
function parse_tz_format(str::AbstractString)
    x = _parsesub_tz(str)
    if x isa Tuple
        tz, i = x
        return tz
    else
        throw(ParseNextError("Unhandled TZ environment variable. $(x.msg)", x.str, x.s, x.e))
    end
end

"""
    tryparse_tz_format(str) -> Union{TimeZone, Nothing}

Like `parse_tz_format`, but returns either a value of the `TimeZone`, or `nothing` if
the string does not contain a valid format.
"""
function tryparse_tz_format(str::AbstractString)
    x = _parsesub_tz(str)
    if x isa Tuple
        tz, i = x
        return tz
    else
        return nothing
    end
end


text:
using Extents
using Test

@testset "Extents.jl" begin
    ex1 = Extent(X=(1, 2), Y=(3, 4))
    ex2 = Extent(Y=(3, 4), X=(1, 2))
    ex3 = Extent(X=(1, 2), Y=(3, 4), Z=(5.0, 6.0))

    @testset "bounds" begin
        @test bounds(ex1) === (X=(1, 2), Y=(3, 4))
        @test bounds(ex2) === (Y=(3, 4), X=(1, 2))
        @test bounds(ex3) === (X=(1, 2), Y=(3, 4), Z=(5.0, 6.0))
    end

    @testset "extent" begin
        @test extent(ex1) === ex1
    end

    @testset "equality" begin
        @test ex1 == ex2
        @test ex1 != ex3
    end

    @testset "properties" begin
        @test keys(ex1) == (:X, :Y)
        @test values(ex1) == ((1, 2), (3, 4))
    end
end


text:
# TODO bounds checking for all indices

abstract type AbstractList{J} <: ArrowVector{J} end

# default offset type to use
const DefaultOffset = Int32

"""
    List{P<:AbstractPrimitive,J} <: AbstractList{J}

An Arrow formatted array of variable length objects such as strings. The `List` contains "value" data
as well as "offsets" which describe from which elements of the values data an element of `List`
should be constructed.  The offsets are necessarily a `Primitive{Int32}` while the values can
be any `ArrowVector` type (but in most circumstances should be `Primitive`).

## Constructors
    List{J}(len::Integer, offs::Primitive{Int32}, vals::AbstractPrimitive)
    List{J}(offs::Primitive{Int32}, vals::AbstractPrimitive)
    List{J}(data::Vector{UInt8}, offset_idx::Integer, len::Integer, vals::AbstractPrimitive)
    List{J}(data::Vector{UInt8}, offset_idx::Integer, values_idx::Integer, ::Type{C}, x::AbstractVector)
    List{J}(data::Vector{UInt8}, i::Integer, ::Type{C}, x::AbstractVector)
    List(data::Vector{UInt8}, i::Integer, ::Type{C}, x::AbstractVector)
    List{J}(Array, ::Type{C}, x::AbstractVector)
    List(Array, ::Type{C}, x::AbstractVector)
    List(::Type{C}, v::AbstractVector)
    List(v::AbstractVector{<:AbstractString})

Note that by default, `List`s of strings will be encoded in UTF-8.

### Arguments
- `len`: the length of the `List`
- `offs`: a `Primitive{Int32}` containing the offsets data
- `vals`: a `Primitive` containing the values data
- `data`: the data buffer in which the underlying data is stored
- `offset_idx`: location within `data` where the offset data is stored
- `values_idx`: location within `data` where the values data is stored
- `C`: the encoding type (type of values), e.g. for UTF-8 strings this is `UInt8`. defaults to `UInt8`
    if not given explicitly
- `i`: the location in `data` where all data should be stored (offsets, then values)
- `x`, `v`: array to be stored or converted
"""
struct List{J,K<:Integer,P<:AbstractPrimitive} <: AbstractList{J}
    length::Int64
    offsets::Primitive{K}
    values::P
end

# Primitive constructors
function List{J}(len::Integer, offs::Primitive{K}, vals::P) where {J,K<:Integer,P<:AbstractPrimitive}
    List{J,K,P}(len, offs, vals)
end
function List{J}(offs::Primitive{K}, vals::P) where {J,K<:Integer,P<:AbstractPrimitive}
    List{J,K,P}(length(offs)-1, offs, vals)
end

function List{J,K}(data::Vector{UInt8}, offset_idx::Integer, len::Integer,
                   vals::P) where {J,K<:Integer,P}
    List{J,K,P}(offs, vals)
end

# all index constructor
function List{J,K}(data::Vector{UInt8}, offset_idx::Integer, values_idx::Integer, len::Integer,
                   ::Type{C}, values_len::Integer) where {J,K<:Integer,C}
    vals = Primitive{C}(data, values_idx, values_len)
    List{J,K}(data, offset_idx, len, vals)
end

# buffer with location constructors, with values arg
function List{J,K}(data::Vector{UInt8}, offset_idx::Integer, len::Integer, vals::P
                  ) where {K<:Integer,P<:AbstractPrimitive,J}
    offs = Primitive{K}(data, offset_idx, len+1)
    List{J}(len, offs, vals)
end

# buffer with location constructors
function List{J,K}(data::Vector{UInt8}, offset_idx::Integer, values_idx::Integer, ::Type{C},
                   x::AbstractVector) where {K<:Integer,C,J}
    offs = Primitive{K}(data, offset_idx, offsets(K, C, x))
    p = Primitive(data, values_idx, encode(C, x))
    List{J}(offs, p)
end
# this puts offsets first
function List{J,K}(data::Vector{UInt8}, i::Integer, ::Type{C}, x::AbstractVector
                  ) where {K<:Integer,C,J}
    offs = Primitive{K}(data, i, offsets(K, C, x))
    p = Primitive(data, i+totalbytes(offs), encode(C, x))
    List{J}(offs, p)
end

function List{J,K}(data::Vector{UInt8}, i::Integer, x::AbstractVector{<:AbstractString}
                  ) where {K<:Integer,J}
    List{J,K}(data, i, UInt8, x)
end

function List{J,K}(::Type{<:Array}, ::Type{C}, x::AbstractVector) where {K<:Integer,C,J}
    b = Vector{UInt8}(undef, totalbytes(K, C, x))
    List{J,K}(b, 1, C, x)
end

function List{J,K}(::Type{<:Array}, x::AbstractVector{<:AbstractString}) where {J,K<:Integer}
    List{J,K}(Array, UInt8, x)
end
function List(::Type{<:Array}, x::AbstractVector{S}) where S<:AbstractString
    List{S,DefaultOffset}(Array, UInt8, x)
end

function List{J,K}(::Type{C}, v::AbstractVector) where {K<:Integer,J,C}
    offs = Primitive{K}(offsets(K, C, v))
    p = Primitive(encode(C, v))
    List{J}(offs, p)
end
List(v::AbstractVector{<:AbstractString}) = List{String,DefaultOffset}(UInt8, v)

List{J,K,P}(l::List{J,K,P}) where {J,K,P} = List{J,K,P}(l.length, l.offsets, l.values)
List{J}(l::List{J}) where J = List{J}(l.length, l.offsets, l.values)
List{J}(l::List{T}) where {J,T} = List{J}(convert(AbstractVector{J}, l[:]))
List(l::List{J}) where J = List{J}(l)

"""
    NullableList{P<:AbstractPrimitive,J} <: AbstractList{Union{Missing,J}}

An Arrow formatted array of variable length objects such as strings which may be null. The `NullableList`
contains a bit mask specifying which values are null and "offsets" which describe from which elements of
the values data an element of the `NullableList` should be constructed.  The bitmask is contained in a
`Primitive{UInt8}` while the offsets data in a `Primitive{Int32}`. The values can be contained in any
`ArrowVector` type, but in most cases should be `Primitive`.

## Constructors
    NullableList{J}(len::Integer, bmask::Primitive, offs::Primitive, vals::AbstractPrimitive)
    NullableList{J}(bmask::Primitive, offs::Primitive, vals::AbstractPrimitive)
    NullableList{J}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer, len::Integer,
                    vals::AbstractPrimitive)
    NullableList{J}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer, values_idx::Integer,
                    len::Integer, ::Type{C}, values_len::Integer)
    NullableList{J}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer, values_idx::Integer,
                    ::Type{C}, x::AbstractVector)
    NullableList(data::Vector{UInt8}, i::Integer, ::Type{C}, x::AbstractVector)
    NullableList(Array, ::Type{C}, x::AbstractVector)
    NullableList(Array, x::AbstracVector)
    NullableList(::Type{C}, v::AbstractVector)
    NullableList(v::AbstractVector)

If `Array` is given as an argument, a contiguous array will be allocated to store the data.

### Arguments
- `len`: the length of the `NullableList`
- `bmask`: the `Primitive` providing the bit mask
- `offs`: the `Primitive` providing the offsets
- `vals`: the `AbstractPrimitive` providing the values
- `data`: a buffer for storing the data
- `bitmask_idx`: the location in `data` of the bit mask
- `offsets_idx`: the location in `data` of the offsets
- `values_idx`: the location in `data` of the values
- `values_len`: the total length of the values data (i.e. number of elements in the values array)
- `C`: the data type of the values data. defaults to `UInt8` when not provided
- `x`, `v`: array to be stored by the `NullableList`
"""
struct NullableList{J,K<:Integer,P<:AbstractPrimitive} <: AbstractList{Union{Missing,J}}
    length::Int
    bitmask::Primitive{UInt8}
    offsets::Primitive{K}
    values::P
end

# Primitive constructors
function NullableList{J}(len::Integer, bmask::Primitive{UInt8}, offs::Primitive{K},
                         vals::P) where {J,K<:Integer,P}
    NullableList{J,K,P}(len, bmask, offs, vals)
end
function NullableList{J}(bmask::Primitive{UInt8}, offs::Primitive{K}, vals::P) where {J,K<:Integer,P}
    NullableList{J,K,P}(length(offs)-1, bmask, offs, vals)
end

function NullableList{J,K}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer,
                           len::Integer, vals::P) where {J,K<:Integer,P}
    bmask = Primitive{UInt8}(data, bitmask_idx, bytesforbits(len))
    offs = Primitive{K}(data, offset_idx, len+1)
    NullableList{J}(bmask, offs, vals)
end

# all index constructor
function NullableList{J,K}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer,
                           values_idx::Integer, len::Integer, ::Type{C}, values_len::Integer
                          ) where {J,K<:Integer,C}
    vals = Primitive{C}(data, values_idx, values_len)
    NullableList{J,K}(data, bitmask_idx, offset_idx, len, vals)
end

# buffer with location constructors, with values arg
function NullableList{J,K}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer,
                           len::Integer, vals::P) where {J,K<:Integer,P<:AbstractPrimitive}
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskbytes(len))
    offs = Primitive{K}(data, offset_idx, len+1)
    NullableList{J}(bmask, offs, vals)
end

# buffer with location constructors
function NullableList{J,K}(data::Vector{UInt8}, bitmask_idx::Integer, offset_idx::Integer,
                           values_idx::Integer, ::Type{C}, x::AbstractVector
                          ) where {C,K<:Integer,J}
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskpadded(x))
    offs = Primitive{K}(data, offset_idx, offsets(K, x))
    vals = Primitive(data, values_idx, encode(C, x))
    NullableList{J}(bmask, offs, vals)
end
# bitmask, offsets, values
function NullableList{J,K}(data::Vector{UInt8}, i::Integer, ::Type{C}, x::AbstractVector
                          ) where {C,K<:Integer,J}
    bmask = Primitive{UInt8}(data, i, bitmaskpadded(x))
    offs = Primitive{K}(data, i+bitmaskbytes(x), offsets(K, C, x))
    vals = Primitive(data, i+bitmaskbytes(x)+totalbytes(offs), encode(C, x))
    NullableList{J}(bmask, offs, vals)
end

function NullableList{J,K}(data::Vector{UInt8}, i::Integer, x::AbstractVector{J}) where {J,K<:Integer}
    NullableList{J,K}(data, i, UInt8, x)
end
function NullableList{J,K}(data::Vector{UInt8}, i::Integer, x::AbstractVector{Union{J,Missing}}
                          ) where {J,K<:Integer}
    NullableList{J,K}(data, i, UInt8, x)
end

function NullableList{J,K}(::Type{<:Array}, ::Type{C}, x::AbstractVector) where {C,K<:Integer,J}
    b = Vector{UInt8}(undef, minbytes(C, x))
    NullableList{J,K}(b, 1, C, x)
end

function NullableList{J,K}(::Type{<:Array}, x::AbstractVector{J}) where {J,K<:Integer}
    NullableList{J,K}(Array, UInt8, x)
end
function NullableList{J,K}(::Type{<:Array}, x::AbstractVector{Union{J,Missing}}) where {J,K<:Integer}
    NullableList{J,K}(Array, UInt8, x)
end
function NullableList(::Type{<:Array}, x::AbstractVector{J}) where {J<:AbstractString}
    NullableList{J,DefaultOffset}(Array, UInt8, x)
end

function NullableList{J,K}(::Type{C}, v::AbstractVector) where {J,K<:Integer,C}
    bmask = Primitive{UInt8}(bitmaskpadded(v))
    offs = Primitive{K}(offsets(K, C, v))
    vals = Primitive(encode(C, v))
    NullableList{J}(bmask, offs, vals)
end
NullableList(::Type{C}, v::AbstractVector{J}) where {C,J} = NullableList{J,DefaultOffset}(C, v)
function NullableList(::Type{C}, v::AbstractVector{Union{J,Missing}}) where {C,J}
    NullableList{J,DefaultOffset}(C, v)
end
NullableList(v::AbstractVector) = NullableList{String,DefaultOffset}(UInt8, v)

function NullableList{J,K,P}(l::NullableList{J,K,P}) where {J,K,P}
    NullableList{J,K,P}(l.length, l.bitmask, l.offsets, l.values)
end
NullableList{J}(l::NullableList{J}) where J = NullableList{J}(l.length, l.bitmask, l.offsets, l.values)
NullableList{J}(l::NullableList{T}) where {J,T} = NullableList{J}(convert(AbstractVector{J}, p[:]))
NullableList(l::NullableList{J}) where J = NullableList{J}(l)

#====================================================================================================
    common interface
====================================================================================================#
function valuesbytes(::Type{C}, A::AbstractVector) where C
    padding(sum(ismissing(a) ? 0 : length(a)*sizeof(C) for a ∈ A))
end
valuesbytes(A::Union{List{P,J},NullableList{P,J}}) where {P,J} = valuesbytes(A.values)

bitmaskbytes(A::List) = 0
bitmaskbytes(A::NullableList) = bytesforbits(length(A))

offsetsbytes(::Type{K}, len::Integer) where {K<:Integer} = padding((len+1)*sizeof(K))
offsetsbytes(::Type{K}, A::AbstractVector) where {K<:Integer} = offsetsbytes(K, length(A))
offsetsbytes(l::AbstractList) = totalbytes(l.offsets)

function totalbytes(::Type{K}, ::Type{C}, A::AbstractVector) where {K<:Integer,C}
    valuesbytes(C, A) + bitmaskbytes(A) + offsetsbytes(K, A)
end
function totalbytes(::Type{Union{J,Missing}}, ::Type{K}, ::Type{C}, A::AbstractVector
                   ) where {J,K<:Integer,C}
    valuesbytes(C, A) + bitmaskbytes(Union{J,Missing}, A) + offsetsbytes(A)
end
totalbytes(A::AbstractList) = valuesbytes(A) + minbitmaskbytes(A) + offsetsbytes(A)

# helper function for offsets
_offsize(::Type{C}, x) where C = sizeof(x)
_offsize(::Type{C}, x::AbstractString) where {C} = sizeof(C)*ncodeunits(x)

# TODO how to deal with sizeof of Arrow objects such as lists?
# note that this works fine with missings because sizeof(missing) == 0
"""
    offsets(v::AbstractVector)

Construct a `Vector{Int32}` of offsets appropriate for data appearing in `v`.
"""
function offsets(::Type{K}, ::Type{C}, v::AbstractVector) where {K<:Integer,C}
    off = Vector{K}(undef, length(v)+1)
    off[1] = 0
    for i ∈ 2:length(off)
        off[i] = _offsize(C, v[i-1]) + off[i-1]
    end
    off
end
offsets(::Type{K}, v::AbstractVector{C}) where {K<:Integer,C} = offsets(K, C, v)
function offsets(::Type{K}, v::AbstractVector{String}) where {K}
    offsets(K, UInt8, v)
end
function offsets(v::AbstractVector{<:AbstractString})
    throw(ArgumentError("must specify encoding type for computing string offsets"))
end

function check_offset_bounds(l::AbstractList, i::Integer)
    if !(1 ≤ i ≤ length(l)+1)
        throw(ArgumentError("tried to access offset $i from list of length $(length(l))"))
    end
end

rawvalues(p::AbstractList, i) = rawvalues(p.values, i)
rawvalues(p::AbstractList) = rawvalues(p.values)

# note that there are always n+1 offsets
"""
    unsafe_getoffset(l::AbstractList, i::Integer)

Get the offset for element `i`.  Contains a call to `unsafe_load`.
"""
unsafe_getoffset(l::AbstractList, i) = unsafe_getvalue(offsets(l), i)

"""
    unsafe_rawoffsets(p::AbstractList)

Retreive the raw offstets for `p` as a `Vector{UInt8}`.
"""
unsafe_rawoffsets(p::AbstractList) = unsafe_rawpadded(offsetspointer(p), offsetsbytes(p))

"""
    getoffset(l::AbstractList, i::Integer)

Retrieve offset `i` for list `l`.  Note that this retrieves the Arrow formated 0-based indexed raw
numbers!
"""
getoffset(l::AbstractList, i) = l.offsets[i]

"""
    unsafe_setoffset!(l::AbstractList, off::Int32, i::Integer)

Set offset `i` to `off`.  Contains a call to `unsafe_store!`.
"""
unsafe_setoffset!(l::AbstractList, off::Int32, i::Integer) = unsafe_setvalue!(offsets(l), off, i)

setoffset!(l::AbstractList, off::Int32, i::Integer) = setindex!(l.offsets, off, i)

"""
    unsafe_setoffsets!(l::AbstractList, off::Vector{Int32})

Set all offsets to the `Vector{Int32}` `off`.  Contains a call to `unsafe_copy!` which copies the
entirety of `off`.
"""
function unsafe_setoffsets!(l::AbstractList, off::Vector{Int32})
    unsafe_copy!(convert(Ptr{Int32}, l.offsets), pointer(off), length(off))
end

setoffsets!(l::AbstractList, off::AbstractVector{Int32}) = (l.offsets[:] = off)

"""
    unsafe_ellength(l::AbstractList, i::Integer)

Get the length of element `i`. Involves calls to `unsafe_load`.
"""
unsafe_ellength(l::AbstractList, i::Integer) = unsafe_getoffset(l, i+1) - unsafe_getoffset(l, i)

"""
    ellength(l::AbstractList, i::Integer)

Get the length of element `i`.
"""
ellength(l::AbstractList, i::Integer) = getoffset(l, i+1) - getoffset(l, i)

# returns offset, length
@inline function unsafe_elparams(l::AbstractList, i::Integer)
    off = unsafe_getoffset(l, i)
    off, unsafe_getoffset(l, i+1) - off
end

@inline function elparams(l::AbstractList, i::Integer)
    off = getoffset(l, i)
    off, getoffset(l, i+1) - off
end

function unsafe_getvalue(l::Union{List{J},NullableList{J}}, i::Integer) where J
    off, len = unsafe_elparams(l, i)
    unsafe_construct(J, values(l), off+1, len)
end
function unsafe_getvalue(l::AbstractList{T}, idx::AbstractVector{<:Integer}) where T
    T[unsafe_getvalue(l, i) for i ∈ idx]
end
function unsafe_getvalue(l::AbstractList{T}, idx::AbstractVector{Bool}) where T
    T[unsafe_getvalue(l, i) for i ∈ 1:length(l) if idx[i]]
end

function getvalue(l::Union{List{J},NullableList{J}}, i::Integer) where J
    off, len = elparams(l, i)
    construct(J, l.values, off+1, len)
end

function getvalue(l::AbstractList{T}, idx::AbstractVector{<:Integer}) where T
    T[getvalue(l, i) for i ∈ idx]
end
function getvalue(l::AbstractList{T}, idx::AbstractVector{Bool}) where T
    T[getvalue(l, i) for i ∈ 1:length(l) if idx[i]]
end


text:
# Sober Physicists Don't Find Giraffes Hiding In Kitchens
const shells = "SPDFGHIK"

# Electronic configurations of the elements, taken from
# https://math.nist.gov/DFTdata/atomdata/configuration.html
const atom = (
H  = c"1s1",
He = c"[He]",
Li = c"[He] 2s1",
Be = c"[He] 2s2",
B  = c"[He] 2s2 2p1",
C  = c"[He] 2s2 2p2",
N  = c"[He] 2s2 2p3",
O  = c"[He] 2s2 2p4",
F  = c"[He] 2s2 2p5",
Ne = c"[Ne]",
Na = c"[Ne] 3s1",
Mg = c"[Ne] 3s2",
Al = c"[Ne] 3s2 3p1",
Si = c"[Ne] 3s2 3p2",
P  = c"[Ne] 3s2 3p3",
S  = c"[Ne] 3s2 3p4",
Cl = c"[Ne] 3s2 3p5",
Ar = c"[Ar]",
K  = c"[Ar] 4s1",
Ca = c"[Ar] 4s2",
Sc = c"[Ar] 3d1 4s2",
Ti = c"[Ar] 3d2 4s2",
V  = c"[Ar] 3d3 4s2",
Cr = c"[Ar] 3d5 4s1",
Mn = c"[Ar] 3d5 4s2",
Fe = c"[Ar] 3d6 4s2",
Co = c"[Ar] 3d7 4s2",
Ni = c"[Ar] 3d8 4s2",
Cu = c"[Ar] 3d10 4s1",
Zn = c"[Ar] 3d10 4s2",
Ga = c"[Ar] 3d10 4s2 4p1",
Ge = c"[Ar] 3d10 4s2 4p2",
As = c"[Ar] 3d10 4s2 4p3",
Se = c"[Ar] 3d10 4s2 4p4",
Br = c"[Ar] 3d10 4s2 4p5",
Kr = c"[Kr]",
Rb = c"[Kr] 5s1",
Sr = c"[Kr] 5s2",
Y  = c"[Kr] 4d1 5s2",
Zr = c"[Kr] 4d2 5s2",
Nb = c"[Kr] 4d4 5s1",
Mo = c"[Kr] 4d5 5s1",
Tc = c"[Kr] 4d5 5s2",
Ru = c"[Kr] 4d7 5s1",
Rh = c"[Kr] 4d8 5s1",
Pd = c"[Kr] 4d10",
Ag = c"[Kr] 4d10 5s1",
Cd = c"[Kr] 4d10 5s2",
In = c"[Kr] 4d10 5s2 5p1",
Sn = c"[Kr] 4d10 5s2 5p2",
Sb = c"[Kr] 4d10 5s2 5p3",
Te = c"[Kr] 4d10 5s2 5p4",
I  = c"[Kr] 4d10 5s2 5p5",
Xe = c"[Xe]",
Cs = c"[Xe] 6s1",
Ba = c"[Xe] 6s2",
La = c"[Xe] 5d1 6s2",
Ce = c"[Xe] 4f1 5d1 6s2",
Pr = c"[Xe] 4f3 6s2",
Nd = c"[Xe] 4f4 6s2",
Pm = c"[Xe] 4f5 6s2",
Sm = c"[Xe] 4f6 6s2",
Eu = c"[Xe] 4f7 6s2",
Gd = c"[Xe] 4f7 5d1 6s2",
Tb = c"[Xe] 4f9 6s2",
Dy = c"[Xe] 4f10 6s2",
Ho = c"[Xe] 4f11 6s2",
Er = c"[Xe] 4f12 6s2",
Tm = c"[Xe] 4f13 6s2",
Yb = c"[Xe] 4f14 6s2",
Lu = c"[Xe] 4f14 5d1 6s2",
Hf = c"[Xe] 4f14 5d2 6s2",
Ta = c"[Xe] 4f14 5d3 6s2",
W  = c"[Xe] 4f14 5d4 6s2",
Re = c"[Xe] 4f14 5d5 6s2",
Os = c"[Xe] 4f14 5d6 6s2",
Ir = c"[Xe] 4f14 5d7 6s2",
Pt = c"[Xe] 4f14 5d9 6s1",
Au = c"[Xe] 4f14 5d10 6s1",
Hg = c"[Xe] 4f14 5d10 6s2",
Tl = c"[Xe] 4f14 5d10 6s2 6p1",
Pb = c"[Xe] 4f14 5d10 6s2 6p2",
Bi = c"[Xe] 4f14 5d10 6s2 6p3",
Po = c"[Xe] 4f14 5d10 6s2 6p4",
At = c"[Xe] 4f14 5d10 6s2 6p5",
Rn = c"[Rn]",
Fr = c"[Rn] 7s1",
Ra = c"[Rn] 7s2",
Ac = c"[Rn] 6d1 7s2",
Th = c"[Rn] 6d2 7s2",
Pa = c"[Rn] 5f2 6d1 7s2",
 U = c"[Rn] 5f3 6d1 7s2");


text:
abstract type Format end

# input interface

"""
    parse(io, format)
    parse(path, format)
    parse(str, format)
"""
parse(path::AbstractPath, format) = parse(open(path), format)
parse(s::String, format) = parse(IOBuffer(s), format)

function parse(path::AbstractPath)
    format = extensionformat(Val(Symbol(extension(path))))
    return parse(path, format)
end

"""
    extensionformat(Val(Symbol(ext)))

Define a default `Format` for parsing files with extension `ext`.

For example `extensionformat(Val(:html)) == HTML()`, so `parse(p"index.html")`
relays to `parse(p"index.html", HTML())`
"""
function extensionformat(::Val{:ext}) end
function formatextension end

# output interface

"""
    render!(io, doc::XExpr, format, tag = Val(doc.tag))
"""
render!(io, doc::XTree, format) = render!(io, doc, format, Val(doc.tag))

function render!(path::AbstractPath, doc::XTree, format)
    open(path, "w") do f
        render!(f, doc, format)
    end
end

function render(doc, format)
    io = IOBuffer()
    render!(io, doc, format)
    return String(take!(io))
end