text:
##########
# Traces #
##########

"""
    Trace{T}

Abstract type for a trace of a generative function.
"""
abstract type Trace end

"""
    get_args(trace)

Return the argument tuple for a given execution.

Example:
```julia
args::Tuple = get_args(trace)
```
"""
function get_args end

"""
    get_retval(trace)

Return the return value of the given execution.

Example for generative function with return type `T`:
```julia
retval::T = get_retval(trace)
```
"""
function get_retval end

"""
    get_choices(trace)

Return a value implementing the assignment interface

Note that the value of any non-addressed randomness is not externally accessible.

Example:

```julia
choices::ChoiceMap = get_choices(trace)
z_val = choices[:z]
```
"""
function get_choices end

"""
    get_score(trace)

Return:
```math
\\log \\frac{p(t, r; x)}{q(r; x, t)}
```

When there is no non-addressed randomness, this simplifies to the log probability `\$\\log p(t; x)\$.
"""
function get_score end

"""
    gen_fn::GenerativeFunction = get_gen_fn(trace)

Return the generative function that produced the given trace.
"""
function get_gen_fn end

"""
    value = trace[addr]

Get the value of the random choice at address `addr`.
"""
function Base.getindex(trace::Trace, addr)
    get_choices(trace)[addr]
end

export get_args
export get_retval
export get_choices
export get_score
export get_gen_fn

######################
# GenerativeFunction #
######################

"""
    GenerativeFunction{T,U <: Trace}

Abstract type for a generative function with return value type T and trace type U.
"""
abstract type GenerativeFunction{T,U <: Trace} end

get_return_type(::GenerativeFunction{T,U}) where {T,U} = T
get_trace_type(::GenerativeFunction{T,U}) where {T,U} = U

"""
    bools::Tuple = has_argument_grads(gen_fn::Union{GenerativeFunction,Distribution})

Return a tuple of booleans indicating whether a gradient is available for each of its arguments.
"""
function has_argument_grads end

"""
    req::Bool = accepts_output_grad(gen_fn::GenerativeFunction)

Return a boolean indicating whether the return value is dependent on any of the *gradient source elements* for any trace.

The gradient source elements are:

- Any argument whose position is true in `has_argument_grads`

- Any trainable parameter

- Random choices made at a set of addresses that are selectable by `choice_gradients`.
"""
function accepts_output_grad end

"""
    get_params(gen_fn::GenerativeFunction)

Return an iterable over the trainable parameters of the generative function.
"""
get_params(::GenerativeFunction) = ()

"""
    trace = simulate(gen_fn, args)

Execute the generative function and return the trace.

Given arguments (`args`), sample \$t \\sim p(\\cdot; x)\$ and \$r \\sim p(\\cdot; x,
t)\$, and return a trace with choice map \$t\$.
"""
function simulate(::GenerativeFunction, ::Tuple)
    error("Not implemented")
end

"""
    (trace::U, weight) = generate(gen_fn::GenerativeFunction{T,U}, args::Tuple)

Return a trace of a generative function.

    (trace::U, weight) = generate(gen_fn::GenerativeFunction{T,U}, args::Tuple,
                                    constraints::ChoiceMap)

Return a trace of a generative function that is consistent with the given
constraints on the random choices.

Given arguments \$x\$ (`args`) and assignment \$u\$ (`constraints`) (which is empty for the first form), sample \$t \\sim
q(\\cdot; u, x)\$ and \$r \\sim q(\\cdot; x, t)\$, and return the trace \$(x, t, r)\$ (`trace`).
Also return the weight (`weight`):
```math
\\log \\frac{p(t, r; x)}{q(t; u, x) q(r; x, t)}
```

Example without constraints:
```julia
(trace, weight) = generate(foo, (2, 4))
```

Example with constraint that address `:z` takes value `true`.
```julia
(trace, weight) = generate(foo, (2, 4), choicemap((:z, true))
```
"""
function generate(::GenerativeFunction, ::Tuple, ::ChoiceMap)
    error("Not implemented")
end

function generate(gen_fn::GenerativeFunction, args::Tuple)
    generate(gen_fn, args, EmptyChoiceMap())
end

"""
    weight = project(trace::U, selection::Selection)

Estimate the probability that the selected choices take the values they do in a
trace.

Given a trace \$(x, t, r)\$ (`trace`) and a set of addresses \$A\$ (`selection`),
let \$u\$ denote the restriction of \$t\$ to \$A\$. Return the weight
(`weight`):
```math
\\log \\frac{p(r, t; x)}{q(t; u, x) q(r; x, t)}
```
"""
function project(trace, selection::Selection)
    error("Not implemented")
end

"""
    (choices, weight, retval) = propose(gen_fn::GenerativeFunction, args::Tuple)

Sample an assignment and compute the probability of proposing that assignment.

Given arguments (`args`), sample \$t \\sim p(\\cdot; x)\$ and \$r \\sim p(\\cdot; x,
t)\$, and return \$t\$
(`choices`) and the weight (`weight`):
```math
\\log \\frac{p(r, t; x)}{q(r; x, t)}
```
"""
function propose(gen_fn::GenerativeFunction, args::Tuple)
    trace = simulate(gen_fn, args)
    weight = get_score(trace)
    (get_choices(trace), weight, get_retval(trace))
end

"""
    (weight, retval) = assess(gen_fn::GenerativeFunction, args::Tuple, choices::ChoiceMap)

Return the probability of proposing an assignment

Given arguments \$x\$ (`args`) and an assignment \$t\$ (`choices`) such that
\$p(t; x) > 0\$, sample \$r \\sim q(\\cdot; x, t)\$ and
return the weight (`weight`):
```math
\\log \\frac{p(r, t; x)}{q(r; x, t)}
```
It is an error if \$p(t; x) = 0\$.
"""
function assess(gen_fn::GenerativeFunction, args::Tuple, choices::ChoiceMap)
    (trace, weight) = generate(gen_fn, args, choices)
    (weight, get_retval(trace))
end

"""
    (new_trace, weight, retdiff, discard) = update(trace, args::Tuple, argdiffs::Tuple,
                                                   constraints::ChoiceMap)

Update a trace by changing the arguments and/or providing new values for some
existing random choice(s) and values for some newly introduced random choice(s).

Given a previous trace \$(x, t, r)\$ (`trace`), new arguments \$x'\$ (`args`), and
a map \$u\$ (`constraints`), return a new trace \$(x', t', r')\$ (`new_trace`)
that is consistent with \$u\$.  The values of choices in \$t'\$ are
either copied from \$t\$ or from \$u\$ (with \$u\$ taking
precedence) or are sampled from the internal proposal distribution.  All choices in \$u\$ must appear in \$t'\$.  Also return an
assignment \$v\$ (`discard`) containing the choices in \$t\$ that were
overwritten by values from \$u\$, and any choices in \$t\$ whose address does
not appear in \$t'\$. Sample \$t' \\sim q(\\cdot; x', t + u)\$, and \$r' \\sim
q(\\cdot; x', t')\$, where \$t + u\$ is the choice map obtained by merging
\$t\$ and \$u\$ with \$u\$ taking precedence for overlapping addresses.  Also
return a weight (`weight`):
```math
\\log \\frac{p(r', t'; x') q(r; x, t)}{p(r, t; x) q(r'; x', t') q(t'; x', t + u)}
```
"""
function update(trace, ::Tuple, argdiffs::Tuple, ::ChoiceMap)
    error("Not implemented")
end

"""
    (new_trace, weight, retdiff) = regenerate(trace, args::Tuple, argdiffs::Tuple,
                                              selection::Selection)

Update a trace by changing the arguments and/or randomly sampling new values
for selected random choices using the internal proposal distribution family.

Given a previous trace \$(x, t, r)\$ (`trace`), new arguments \$x'\$ (`args`), and
a set of addresses \$A\$ (`selection`), return a new trace \$(x', t')\$
(`new_trace`) such that \$t'\$ agrees with \$t\$ on all addresses not in \$A\$
(\$t\$ and \$t'\$ may have different sets of addresses).  Let \$u\$ denote the
restriction of \$t\$ to the complement of \$A\$.  Sample \$t' \\sim Q(\\cdot;
u, x')\$ and sample \$r' \\sim Q(\\cdot; x', t')\$.
Return the new trace \$(x', t', r')\$ (`new_trace`) and the weight
(`weight`):
```math
\\log \\frac{p(r', t'; x') q(t; u', x) q(r; x, t)}{p(r, t; x) q(t'; u, x') q(r'; x', t')}
```
where \$u'\$ is the restriction of \$t'\$ to the complement of \$A\$.
"""
function regenerate(trace, args::Tuple, argdiffs::Tuple, selection::Selection)
    error("Not implemented")
end

"""
    arg_grads = accumulate_param_gradients!(trace, retgrad, scale_factor=1.)

Increment gradient accumulators for parameters by the gradient of the
log-probability of the trace, optionally scaled, and return the gradient with
respect to the arguments (not scaled).

Given a previous trace \$(x, t)\$ (`trace`) and a gradient with respect to the
return value \$∇_y J\$ (`retgrad`), return the following gradient (`arg_grads`)
with respect to the arguments \$x\$:
```math
∇_x \\left( \\log P(t; x) + J \\right)
```
Also increment the gradient accumulators for the trainable parameters \$Θ\$ of
the function by:
```math
∇_Θ \\left( \\log P(t; x) + J \\right)
```
"""
function accumulate_param_gradients!(trace, retgrad, scale_factor)
    error("Not implemented")
end

accumulate_param_gradients!(trace, retgrad) = accumulate_param_gradients!(trace, retgrad, 1.)

"""
    (arg_grads, choice_values, choice_grads) = choice_gradients(trace, selection::Selection,
                                                                retgrad)

Given a previous trace \$(x, t)\$ (`trace`) and a gradient with respect to the
return value \$∇_y J\$ (`retgrad`), return the following gradient (`arg_grads`)
with respect to the arguments \$x\$:
```math
∇_x \\left( \\log P(t; x) + J \\right)
```
Also given a set of addresses \$A\$ (`selection`) that are continuous-valued
random choices, return the folowing gradient (`choice_grads`) with respect to
the values of these choices:
```math
∇_A \\left( \\log P(t; x) + J \\right)
```
Also return the assignment (`choice_values`) that is the restriction of \$t\$ to \$A\$.
"""
function choice_gradients(trace, selection::Selection, retgrad)
    error("Not implemented")
end

export GenerativeFunction
export has_argument_grads
export accepts_output_grad
export get_params
export simulate
export generate
export project
export propose
export assess
export update
export regenerate
export accumulate_param_gradients!
export choice_gradients

max_stars_repo_name: cameronraysmith/Gen

max_issues_repo_name: cameronraysmith/Gen

max_forks_repo_name: cameronraysmith/Gen


text:
module ExamplePkg

function square(x::Number)
    return x*x
end

end # module

ext: jl


text:
using UnivariateFunctions: evaluate, years_between, years_from_global_base, PE_Function, Sum_Of_Functions, change_base_of_PE_Function, derivative, indefinite_integral, right_integral, left_integral, evaluate_integral
using Dates
tol = 10*eps()

today = Date(2000,1,1)
pe_func = PE_Function(1.0,2.0,today, 3)
(pe_func.base_ - years_from_global_base(today))   < tol
date_in_2020 = Date(2020,1,1)
pe_func2 = PE_Function(1.0,2.0,date_in_2020, 3)
(pe_func2.base_ - years_from_global_base(date_in_2020))   < tol
abs(evaluate(pe_func, date_in_2020) - evaluate(pe_func, years_from_global_base(date_in_2020)) ) < tol

#Sum of functions
sum_func = Sum_Of_Functions([pe_func, PE_Function(2.0,2.5,today, 3) ])
abs(evaluate(sum_func, date_in_2020) - evaluate(sum_func, years_from_global_base(date_in_2020)) ) < tol

# left and right integrals
l_int = left_integral(pe_func, today)
(evaluate(l_int, date_in_2020) - evaluate_integral(pe_func,today,date_in_2020)) < tol
r_int = right_integral(pe_func, date_in_2020)
(evaluate(r_int, today) - evaluate_integral(pe_func,today,date_in_2020)) < tol

# With DateTimes
today_time = DateTime(2000,1,1, 10, 10, 1)
pe_func = PE_Function(1.0,2.0,today_time, 3)

# left and right integrals
later_today = DateTime(2000,1,1, 10, 14, 1)
l_int = left_integral(pe_func, today)
(evaluate(l_int, later_today) - evaluate_integral(pe_func,today,later_today)) < tol
r_int = right_integral(pe_func, later_today)
(evaluate(r_int, today) - evaluate_integral(pe_func,today,later_today)) < tol


text:
using NonlinearBandits
using Documenter

DocMeta.setdocmeta!(
    NonlinearBandits, :DocTestSetup, :(using NonlinearBandits); recursive=true
)

makedocs(;
    modules=[NonlinearBandits],
    authors="Doug Corbin <dfcorbin98@gmail.com> and contributors",
    repo="https://github.com/dfcorbin/NonlinearBandits.jl/blob/{commit}{path}#{line}",
    sitename="NonlinearBandits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dfcorbin.github.io/NonlinearBandits.jl",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Tutorials" => ["Bandits Tutorial" => "bandit_tutorial.md"],
        "API" => ["Bandits" => "bandits_api.md", "Models" => "model_api.md"],
    ],
)

deploydocs(; repo="github.com/dfcorbin/NonlinearBandits.jl", devbranch="main")


text:
using DataDepsGenerators
using Test

@test remove_cite_version(citation_text("10.3732/ajb.1000481")) ==
    "Smith, S. A., Beaulieu, J. M., Stamatakis, A., & Donoghue, M. J. (2011). Understanding angiosperm diversification using small and large phylogenetic trees. American Journal of Botany, 98(3), 404–414. doi:10.3732/ajb.1000481"

@test remove_cite_version(citation_text("https://doi.org/10.5061/dryad.8790")) ==
    "Smith, S. A., Beaulieu, J. M., Stamatakis, A., & Donoghue, M. J. (2011). Data from: Understanding angiosperm diversification using small and large phylogenetic trees [Data set]. Dryad Digital Repository. https://doi.org/10.5061/dryad.8790"


text:
for op = [:add, # :sub, :max, :min,
          ]
    @eval function $(Symbol("scatter_", op, "!"))(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{Int}) where T
        function kernel!(ys, us, xs)
            li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if li <= length(xs) && i <= size(ys, 1)
                ind = Tuple(CartesianIndices(xs)[li])
                CUDAnative.$(Symbol("atomic_", op, "!"))(
                    pointer(ys,
                            Base._to_linear_index(ys,
                                                  i, xs[li]
                                                  )
                            ),
                    us[i, ind...]
                )
            end

            return
        end

        max_threads = 256
        thread_x = min(max_threads, size(ys, 1))
        thread_y = min(max_threads ÷ thread_x, length(xs))
        threads = (thread_x, thread_y)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end

    @eval function $(Symbol("scatter_", op, "!"))(ys::CuArray{T}, us::CuArray{T}, xs::CuArray{<:Tuple}) where T
        function kernel!(ys, us, xs)
            li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if li <= length(xs) && i <= size(ys, 1)
                ind = Tuple(CartesianIndices(xs)[li])
                CUDAnative.$(Symbol("atomic_", op, "!"))(
                    pointer(ys,
                            Base._to_linear_index(ys,
                                                  i, xs[li]...
                                                  )
                            ),
                    us[i, ind...]
                )
            end

            return
        end

        max_threads = 256
        thread_x = min(max_threads, size(ys, 1))
        thread_y = min(max_threads ÷ thread_x, length(xs))
        threads = (thread_x, thread_y)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end
end


text:
"""
    Nlplex{T <: Real} = Nilplexification{T}

A nilplex number is a nilplexification of the real numbers.
"""
const Nlplex{T <: Real} = Nilplexification{T}

function Nlplex(a::T, b::T) where T <: Real
    Nlplex{T}(a, b)
end

function Nlplex(a::T) where T <: Real
    Nlplex{T}(a, zero(T))
end

function Nlplex(a::Real, b::Real)
    Nlplex(promote(a, b)...)
end

"""
    unreal(z::Nlplex)

The unreal part of a nilplex number is a 1-vector.
"""
function unreal(z::Nlplex{T}) where T <: Real
    T[z.r]
end

function show(io::IO, z::Nlplex)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", A: ")
    print(io, z.r)
    print(io, "]")
end

function random(::Type{Nlplex{T}}) where T <: Real
    Nlplex{T}(rand(T), rand(T))
end

function random(::Type{Nlplex{BigInt}})
    Nlplex{BigInt}(randomBigInt(), randomBigInt())
end

function random(::Type{Nlplex{BigFloat}})
    Nlplex{BigFloat}(randomBigFloat(), randomBigFloat())
end

function random(::Type{Nlplex{Rational{BigInt}}})
    Nlplex{Rational{BigInt}}(randomBigRational(), randomBigRational())
end

function abs2(z::Nlplex)
    abs2(z.l)
end

function barA(z::Nlplex)
    Nlplex(z.l, -(z.r))
end

function iszerodivisor(z::Nlplex)
    iszero(z.l)
end

function inv(z::Nlplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    barA(z) / abs2(z)
end


text:
function u_pack(x, box::Box, cl::CellList; parallel=false)
    cl = UpdateCellList!(x, box, cl; parallel=parallel)
    u = map_pairwise!(
        (x, y, i, j, d2, u) -> begin
            u += (sqrt(d2) - box.cutoff)^2 # objective function
            return u
        end,
        0.0,
        box,
        cl;
        parallel=parallel,
    )
    return u
end

function gradient_descent!(x::Vector{T}, f, g!; tol=1e-5, maxtrial=50_000) where {T}
    gnorm(x) = maximum(norm(v) for v in x)
    itrial = 0
    step = 1.0
    xtrial = similar(x)
    g = fill!(similar(x), zero(T))
    fx = f(x)
    g = g!(g, x)
    while (gnorm(g) > tol) && (itrial < maxtrial)
        @. xtrial = x - step * g
        ftrial = f(xtrial)
        if ftrial >= fx
            step = step / 2
        else
            x .= xtrial
            fx = ftrial
            g = g!(g, x)
            step = step * 2
        end
        itrial += 1
    end
    return x
end

function fpair_cl(x, y, i, j, d2, f, box::Box)
    Δv = y - x
    d = sqrt(d2)
    fₓ = 2.0 * (d - box.cutoff) * (Δv / d)
    f[i] += fₓ
    f[j] -= fₓ
    return f
end

function forces_cl!(f::Vector{T}, x, box::Box, cl::CellList, fpair::F) where {T,F}
    fill!(f, zero(T))
    cl = UpdateCellList!(x, box, cl; parallel=false)
    map_pairwise!(
        (x, y, i, j, d2, f) -> fpair(x, y, i, j, d2, f, box), f, box, cl; parallel=false
    )
    return f
end

function packpositions(positions, box::Box; parallel=false)
    box_pack = Box(box.unit_cell_max, 0.5)
    cl_pack = CellList(positions, box_pack; parallel=parallel)
    xpos = gradient_descent!(
        copy(positions),
        (x) -> u_pack(x, box_pack, cl_pack),
        (g, x) -> -forces_cl!(g, x, box_pack, cl_pack, fpair_cl),
    )

    return xpos, box_pack, cl_pack
end

function packsystem!(system::System, cl::CellList, uij; parallel=false)
    uenergy = map_pairwise!(uij, 0.0, system.box, cl; parallel=parallel) / system.npart
    println("initial energy $(uenergy)")

    (pack_pos, boxpack, clpack) = packpositions(copy(system.xpos), system.box)
    cl = UpdateCellList!(pack_pos, system.box, cl; parallel=parallel)
    upacked = u_pack(pack_pos, boxpack, clpack; parallel=parallel)
    system.xpos = copy(pack_pos)
    @show upacked

    uenergy = map_pairwise!(uij, 0.0, system.box, cl; parallel=parallel) / system.npart
    println("initial energy $(uenergy)")

    return nothing
end


text:
"""
   DiscretizeUniformWidth(alg::Symbol)
If `nbins` is a symbol, automatically determine the number of bins to use
"""
struct DiscretizeUniformWidth <: DiscretizationAlgorithm
    nbins::Union{Int,Symbol}
end

function get_nbins(alg::Symbol, data::AbstractArray{N}) where {N<:Real}

    n = length(data)

    if alg == :sqrt
        # Square root (of data size) estimator, used by Excel and other programs for its speed and simplicity.
        nbins = ceil(Int, sqrt(n))
    elseif alg == :sturges
        # R’s default method, only accounts for data size.
        # Only optimal for gaussian data and underestimates number of bins for large non-gaussian datasets.
        # It implicitly bases the bin sizes on the range of the data and can perform poorly if n < 30,
        # because the number of bins will be small—less than seven—and unlikely to show trends in the data well.
        # It may also perform poorly if the data are not normally distributed.
        nbins = ceil(Int, log(2,n)) + 1
    elseif alg == :rice
        # Estimator does not take variability into account, only data size. Commonly overestimates number of bins required.
        nbins = ceil(Int, 2cbrt(n))
    elseif alg == :doane
        # An improved version of Sturges’ estimator that works better with non-normal datasets.
        g₁ = moment(data, 3)
        σ = sqrt((6*(n-2))/((n+1)*(n+3)))
        nbins = ceil(Int, 1 + log(2,n) + log(2, 1+abs(g₁)/σ))
    elseif alg == :scott
        # Less robust estimator that that takes into account data variability and data size.
        σ = std(data)
        binwidth = 3.5σ/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil(Int, (hi - lo)/binwidth)
    elseif alg == :fd # Freedman Diaconis Estimator
        # Robust (resilient to outliers) estimator that takes into account data variability and data size
        binwidth = 2iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil(Int, (hi - lo)/binwidth)
    else # alg == :auto
        # Maximum of the ‘sturges’ and ‘fd’ estimators. Provides good all round performance

        binwidth = 2iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins_fd = ceil(Int, (hi - lo)/binwidth)
        nbins_sturges = ceil(Int, log(2,n)) + 1
        nbins = max(nbins_fd, nbins_sturges)
    end

    nbins
end

function binedges(alg::DiscretizeUniformWidth, data::AbstractArray{N}) where {N<:AbstractFloat}
    lo, hi = extrema(data)
    @assert(hi > lo)

    nbins = (isa(alg.nbins, Symbol) ? get_nbins(alg.nbins, data) : alg.nbins)::Int

    convert(Vector{N}, collect(range(lo, stop=hi, length=nbins+1)))
end
function binedges(alg::DiscretizeUniformWidth, data::AbstractArray{N}) where {N<:Integer}
    lo, hi = extrema(data)
    @assert(hi > lo)
    collect(range(lo, stop=hi, length=alg.nbins+1))
end


text:
"""
# module ReadInput

# Examples

```jldoctest
julia>
```
"""
module ReadInput

using Parsso.CoreDataStructure

const NAMELIST_IDENTIFIERS = Dict(
    program => map(name -> string("&", name), names) for (program, names) in ALLOWED_NAMELISTS
)

const CARD_IDENTIFIERS = ALLOWED_CARDS

end


text:
include("RectangularRoom.jl")
using Plots

# Let's define the room size along x, y, and z, in meters
Lx = 5.0
Ly = 4.0
Lz = 3.0

# Here we define the step of the grid over which we want to calculate the mode
# shape. It will be one point every rx, ry, rz meters along x, y, and z
# respectively.
rx = 0.1
ry = 0.1
rz = 0.1

# We identify the mode we want to calculate with these three indices, the mode
# numbers.
nx = 1
ny = 0
nz = 0

# Here we assign an amplitude to the mode, as if we were plotting the pressure
# field in the room as driven by a sine wave source of appropriate power to
# sustain the target SPL, at steady state regime. dB re 20uPa
SPL = 80.0

# Finally, we set the height in the room at which we want to plot the result.
lh = 1.80

# Computations below.

x, y, z = roomAxes(rx, ry, rz, Lx, Ly, Lz)
X, Y, Z = meshGrid(x, y, z)
li      = argmin(abs.(z .- lh))

P       = 20e-6 * exp10(SPL / 20.0)

eigenV  = P * mode.(X, Y, Z, nx, ny, nz, Lx, Ly, Lz)
eigenF  = modeFrequency(nx, ny, nz, Lx, Ly, Lz)

clibrary(:colorbrewer)

levs    = range(-P, stop = P, length = 128)

# Let's make a little plot
modePlot = contourf(x, y, eigenV[:, :, li],
        levels = levs,
        xlabel = "x [m]",
        ylabel = "y [m]",
        colorbar_title = "Pressure Disturbance [Pa]",
        title = "Resonance mode at $eigenF Hz, z = $(z[li]) m"
        )

### Can make a little animation too - Comment this section in case of troubles
eigenT  = 1.0 / eigenF

Tsteps  = 64

nT      = 1

aL      = nT * eigenT

Fs      = Tsteps * eigenF

nF      = floor(Integer, aL * Fs)

anim = @animate for i = 1:nF
    contourf(x, y, eigenV[:, :, li] * cospi(2 * eigenF * (i - 1) / Fs),
        levels = levs,
        xlabel = "x [m]",
        ylabel = "y [m]",
        colorbar_title = "Pressure Disturbance [Pa]",
        title = "Resonance mode at $eigenF Hz, z = $(z[li]) m"
        )
end

fps     = 12

gif(anim, "/tmp/rmode.gif", fps = fps)

###

# Let's show the first plot:
display(modePlot)

max_stars_repo_name: CrocoDuckoDucks/Acoustics-Models

max_issues_repo_name: CrocoDuckoDucks/RectangularRoom

max_forks_repo_name: CrocoDuckoDucks/RectangularRoom


text:
module HistoricTests

include("utils.jl")
include("samples.jl")
include("test_samples.jl")
include("test_logging.jl")
include("test_defaultlogger.jl")
include("test_id.jl")
include("test_uniqueid.jl")

end  # module HistoricTests


text:
export ## Types
       Field,
       MicrowaveControl,
       QuadratureControl,
       ## Methods
       amplitude,
       label,
       load_sequence!

#AWG channels are controls

#A double-balanced mixer driven by an RF/microwave souce and a single AWG channel
type MicrowaveControl <: Control
    label::String
    freq::Float64
    phase::Float64
    timeStep::Float64
    sequence::InterpGrid
end
MicrowaveControl(label, freq; phase=0., timeStep=1/1.2, sequence=InterpGrid([0.0], BCnearest, InterpNearest)) = MicrowaveControl(label, freq, phase, timeStep, sequence)

# Methods to create the interpolated sequence object
function load_sequence!{T}(mc::MicrowaveControl, sequence::Vector{T})
    mc.sequence = InterpGrid(sequence, BCnearest, InterpNearest)
end

function load_sequence!(mc::MicrowaveControl, seqDict::Dict, n::Int)
    #For double balanced mixer need one AWG channel
    #Hack around off the expected BBNAPSx-xx style
    chan = label(mc)[end]
    APSLabel = label(mc)[1:7]
    mc.sequence = InterpGrid(seqDict[APSLabel]["chan_"*string(chan)][n], BCnearest, InterpNearest)
    return nothing
end
function amplitude(mc::MicrowaveControl, t::Float64)
    #Scale time by the timestep before interpolating
    return mc.sequence[t/mc.timeStep]*cos(2*pi*mc.freq*t + mc.phase)
end

#A pair of AWG channels driving an IQ mixer with a microwave source at a given frequency
type QuadratureControl <: Control
    label::String
    freq::Float64
    phase::Float64
    timeStep::Float64
    sequence_I::InterpGrid
    sequence_Q::InterpGrid
end
QuadratureControl(label, freq; phase=0., timeStep=1/1.2, sequence_I=InterpGrid([0.0], BCnearest, InterpNearest), sequence_Q=InterpGrid([0.0], BCnearest, InterpNearest)) = QuadratureControl(label, freq, phase, timeStep, sequence_I, sequence_Q)

#Create the interpolated object
function load_sequence!(qc::QuadratureControl, seqDict::Dict, n::Int)
    #For quadrature controls we need two channels
    #Hack around off the expected BBNAPSx-xx style
    chan_I = label(qc)[end-1]
    chan_Q = label(qc)[end]
    APSLabel = label(qc)[1:7]
    qc.sequence_I = InterpGrid(seqDict[APSLabel]["chan_"*string(chan_I)][n], BCnearest, InterpNearest)
    qc.sequence_Q = InterpGrid(seqDict[APSLabel]["chan_"*string(chan_Q)][n], BCnearest, InterpNearest)
    return nothing
end
function amplitude(qc::QuadratureControl, t::Float64)
    #Create the complex I/Q pair
    #Scale time by the timestep
    phasor = qc.sequence_I[t/qc.timeStep] - 1im*qc.sequence_Q[t/qc.timeStep]
    return abs(phasor)*cos(2*pi*qc.freq*t + qc.phase + angle(phasor))
end

#AWG controls get connected, possibly through a transfer function, to fields
type Field
    control::Control
end
#Pull in amplitude from controls
#[todo] - allow for transfer function
amplitude(f::Field, t::Float64) = amplitude(f.control, t)


text:
module TestParameter

using ModiaLang
using DifferentialEquations
@usingModiaPlot
include("$(ModiaLang.path)/models/Blocks.jl")

inputSignal(t) = sin(t)

FirstOrder = Model(
    T = parameter | Var(value=0.2),
    x = Var(init=0.3),
    equations = :[u = inputSignal(time/u"s"),
                  T * der(x) + x = u,
                  y = 2*x]
)

firstOrder = @instantiateModel(FirstOrder)

simulate!(firstOrder, Tsit5(), stopTime = 10, log=false, logParameters=false, logEvaluatedParameters=false, requiredFinalStates = [-0.3617373025974107])

plot(firstOrder, ["u", "x", "der(x)", "y"], figure=1)

# Second order system:
#     w = 1.0
#     D = 0.1
#     k = 2.0
#     der(x1) = x2
#     der(x2) = -w^2*x1 - 2*D*w*x2 + w^2*u
#          y  = k*x1
SecondOrder = Model(
    w = 20.0,
    D =  0.1,
    k =  2.0,
    sys = StateSpace | Map(A = parameter | :([  0        1;
                                              -w^2  -2*D*w]),
                           B = parameter | :([0; w^2]),
                           C = parameter | :([k 0]),
                           D = parameter | :(zeros(1)),
                           x = Var(init = zeros(2))),
    equations = :[sys.u = 1.0]
)

using Test
@test_skip begin
secondOrder = @instantiateModel(SecondOrder, unitless=true)

simulate!(secondOrder, merge=Map(D=0.3), logParameters=true, logEvaluatedParameters=true,
          requiredFinalStates = [0.9974089572681231, 0.011808652820321723])
plot(secondOrder, [("sys.u", "sys.y"); "sys.x"], figure = 2)

simulate!(secondOrder, merge=Map(w=100.0, D=0.1, k=1.1), logParameters=true, logEvaluatedParameters=true,
          requiredFinalStates = [0.9999806306477742, -0.00391695258236727])
plot(secondOrder, [("sys.u", "sys.y"); "sys.x"], figure = 3)
end
end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
@testset "90.subsets-ii.jl" begin
    @test subsets_with_dup([1, 2, 2]) == [Int64[], [1], [2], [1, 2], [2, 2], [1, 2, 2]]
    @test subsets_with_dup([1, 2, 2, 2]) ==
          [Int64[], [1], [2], [1, 2], [2, 2], [1, 2, 2], [2, 2, 2], [1, 2, 2, 2]]
end


text:
using Missings
using CSV
using DataStreams, WeakRefStrings, CategoricalArrays
using DataFrames
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end
if VERSION < v"0.7.0-DEV.2575"
    using Base.Dates
else
    using Dates
end

const dir = joinpath(dirname(@__FILE__),"test_files/")
# dir = joinpath(Pkg.dir("CSV"), "test/test_files")

@testset "CSV" begin

include("parsefields.jl")
include("io.jl")

include("source.jl")
include("sink.jl")
include("multistream.jl")
include("validate.jl")

end


text:
module Prinz

using ModelingToolkit
using ..NeuronBuilder

import ..channel_dynamics
include("channels.jl")

export channel_dynamics

end


text:
#!/usr/bin/env julia

"""
Check integer size of a blas library.
Seems like all BLAS routines return just positive indices.
"""

const blas_prec = (:ilaprec_64_, BLAS.libblas)
ccall(blas_prec, Clong, (Cstring,), "X")

const libblas = "/usr/lib/libblas.so.3"
const openblas = "/usr/lib/libopenblas.so.0"

using Base.LinAlg.BlasInt

A = [1.0 0.0; 0.0 -1.0]
lda = max(1,stride(A,2))
infot = Ref{BlasInt}()
uplo = 'U'
ccall((:dpotrf_64_, Base.LinAlg.BLAS.libblas), Void,
      (Ptr{UInt8}, Ptr{BlasInt}, Ptr{eltype(A)}, Ptr{BlasInt}, Ptr{BlasInt}),
      &uplo, &size(A,1), A, &lda, info)

ccall((:dpotrf_, openblas), Void,
      (Ptr{UInt8}, Ptr{BlasInt}, Ptr{eltype(A)}, Ptr{BlasInt}, Ptr{BlasInt}),
      &uplo, &size(A,1), A, &lda, infot)


text:
module TestInteractiveCodeSearch

include("preamble.jl")
using InteractiveCodeSearch:
    Shallow, Recursive, list_locatables, module_methods, choose_method,
    read_stdout, parse_loc, single_macrocall, isliteral
using Base: find_source_file

function with_config(f; kwargs...)
    config = deepcopy(InteractiveCodeSearch.CONFIG)
    try
        for (name, value) in kwargs
            setfield!(InteractiveCodeSearch.CONFIG, name, value)
        end
        f()
    finally
        for name in fieldnames(typeof(config))
            value = getfield(config, name)
            setfield!(InteractiveCodeSearch.CONFIG, name, value)
        end
    end
end

@testset "read_stdout" begin
    @test strip(String(read_stdout("spam", `cat`))) == "spam"
    @test strip(String(read_stdout(io -> write(io, "egg"), `cat`))) == "egg"
end

@testset "parse_loc" begin
    @test parse_loc("@test(ex::ANY, kws...) in Base.Test at test.jl:249") ==
        ("test.jl", 249)
end

module ModuleA
    module ModuleB
        b() = nothing
    end
    a() = nothing
end

@testset "list_locatables" begin
    locs = Set(list_locatables(Shallow(), InteractiveCodeSearch))
    @test Set([
        InteractiveCodeSearch.list_locatables,
        InteractiveCodeSearch.module_methods,
        InteractiveCodeSearch.read_stdout,
        InteractiveCodeSearch.parse_loc,
        InteractiveCodeSearch.search_methods,
        getfield(InteractiveCodeSearch, Symbol("@search")),
    ]) <= locs

    locs = Set(list_locatables(Shallow(), ModuleA))
    @test locs >= Set(Any[ModuleA.a])
    locs = Set(list_locatables(Shallow(), ModuleA.ModuleB))
    @test locs >= Set(Any[ModuleA.ModuleB.b])
    if VERSION >= v"0.7-"
        locs = Set(list_locatables(Recursive(), ModuleA))
        @test locs >= Set(Any[ModuleA.a, ModuleA.ModuleB.b])
    end
end

@testset "module_methods" begin
    @testset for
            p in [Shallow(), Recursive()],
            m in [InteractiveCodeSearch, Base.Filesystem]
        @test_nothrow module_methods(p, m)
    end
end

@testset "choose_method" begin
    single_method = module_methods  # has only one method

    with_config(
        open = (_...) -> error("must not be called"),
        interactive_matcher = `echo " at test.jl:249"`,
        auto_open = true,
    ) do
        # when function has only one method, `auto_open` has to kick-in:
        path, line = choose_method(methods(single_method))
        @test path isa String
        @test path != "test.jl"
        @test line isa Integer
        @test line != 249

        # `code_search` has several methods, so the matcher has to be called:
        path, line = choose_method(methods(InteractiveCodeSearch.code_search))
        @test (path, line) == ("test.jl", 249)

        # "Not opening uninteresting location"
        chosen = choose_method(methods(InteractiveCodeSearch._Dummy, (Any,)))
        @eval @test ($chosen) === nothing
    end

    with_config(
        open = (_...) -> error("must not be called"),
        interactive_matcher = `echo " at test.jl:249"`,
        auto_open = false,
    ) do
        # When `auto_open = false`, the matcher has to be called
        path, line = choose_method(methods(single_method))
        @test (path, line) == ("test.jl", 249)
    end
end

@testset "find_source_file" begin
    m = first(methods(Pkg.clone))
    file = string(m.file)
    found = InteractiveCodeSearch.find_source_file(file)
    @eval @test isfile($found)
end

@testset "isliteral" begin
    @test isliteral(1)
    @test isliteral(1.0)
    @test isliteral("1")
    @test ! isliteral(:a)
    @test ! isliteral(:(a + b))
end

@testset "single_macrocall" begin
    @test single_macrocall(:(@search)) == Symbol("@search")
    @test single_macrocall(quote @search end) == Symbol("@search")
    @test single_macrocall(:f) == nothing
    @test single_macrocall(quote f end) == nothing
end

@testset "patched" begin
    open_args = []
    dummy_openline(args...) = push!(open_args, args)

    with_config(
        interactive_matcher = `echo " at test.jl:249"`,
        open = dummy_openline,
        auto_open = false,
    ) do
        @test_nothrow @eval @search read_stdout
        @test_nothrow @eval @search read_stdout(`cat`, "hello")
        @test_nothrow @eval @search read_stdout(::Cmd, ::String)
        @test_nothrow @eval @search read_stdout(`cat`::Cmd, "hello"::String)
        @test_nothrow @eval @search @search
        @test_nothrow @eval @search @search(read_stdout)
        @test_nothrow @eval @search ""
        @test_nothrow @eval @search InteractiveCodeSearch
        @test_nothrow @eval @searchmethods im
        @test_nothrow @eval @searchmethods ::Complex
        @test_nothrow @eval @searchmethods c::Complex
        @test open_args == repeat([(find_source_file("test.jl"), 249)],
                                  outer=11)

        # @show open_args
    end

    with_config(
        interactive_matcher = `true`,
        open = (args...) -> error("open must not be called"),
        auto_open = false,
    ) do
        @test_nothrow @eval @search read_stdout
    end
end

if v"0.7-" <= VERSION < v"1.2-"
    include("test_return.jl")
end

end  # module


text:
function reconstruct_U(U, parameters::GeneralParams)
    @unpack N1, N2, M1, M2 = parameters

    reconstr = reshape(U, (N1 + 1) * M2, 3)
    datah = reshape(reconstr[:,1], M2, N1 + 1)'
    dataq = reshape(reconstr[:,2], M2, N1 + 1)'
    datap = reshape(reconstr[:,3], M2, N1 + 1)'

    datah   = ffti(datah * M1 * M2, parameters)
    dataq   = ffti(dataq * M1 * M2, parameters)
    datap   = ffti(datap * M1 * M2, parameters)

    return datah, dataq, datap
end

function reconstruct_velocity_fields(U, kx, kz, ScaleH, ScaleU, parameters::GeneralParams)

    @unpack M1, M2, N1, N2, h0 = parameters

    sizeU = size([1, :]);

    # assumption: simplified model
    reconstru = reshape(U[1, :], (N1 + 1) * M2, 3)
    dataS1 = zeros(size(reconstru, 1), 1)
    dataS2 = dataS1;
    dataR1 = dataS2;
    dataR2 = dataS2;

    datah = reconstru[:,1];
    dataq = reconstru[:,2];
    datap = reconstru[:,3];

    datah = reshape(datah, M2, N1 + 1)'
    dataq = reshape(dataq, M2, N1 + 1)'
    datap = reshape(datap, M2, N1 + 1)'
    dataS1 = reshape(dataS1, M2, N1 + 1)'
    dataS2 = reshape(dataS2, M2, N1 + 1)'
    dataR1 = reshape(dataR1, M2, N1 + 1)'
    dataR2 = reshape(dataR2, M2, N1 + 1)'

    datah[1,1] = h0 + imag(datah[1,1])

    kx2 = kx .* kx;
    kz2 = kz .* kz;
    kxz = kx .* kz;
    k2 = kx2 + kz2;
    M1M2 = M1 * M2;

    # refactor: derivatives are duplicated in derivative.jl
    # create a re-useable derivatives function

    # x-derivatives
    datahx = kx .* imag(datah) - kx .* real(datah) * 1im;
    dataqx = kx .* imag(dataq) - kx .* real(dataq) * 1im;
    datapx = kx .* imag(datap) - kx .* real(datap) * 1im;
    dataS1x = kx .* imag(dataS1) - kx .* real(dataS1) * 1im;
    dataS2x = kx .* imag(dataS2) - kx .* real(dataS2) * 1im;
    dataR1x = kx .* imag(dataR1) - kx .* real(dataR1) * 1im;
    dataR2x = kx .* imag(dataR2) - kx .* real(dataR2) * 1im;

    hx  = real(ffti(datahx * M1M2, parameters))
    px  = real(ffti(datapx * M1M2, parameters));
    qx  = real(ffti(dataqx * M1M2, parameters));
    s1x  = real(ffti(dataS1x * M1M2, parameters));
    s2x  = real(ffti(dataS2x * M1M2, parameters));
    r1x  = real(ffti(dataR1x * M1M2, parameters));
    r2x  = real(ffti(dataR2x * M1M2, parameters));

    # z-derivativess
    datahz = kz .* imag(datah) - kz .* real(datah) * 1im;
    dataqz = kz .* imag(dataq) - kz .* real(dataq) * 1im;
    datapz = kz .* imag(datap) - kz .* real(datap) * 1im;
    dataS1z = kz .* imag(dataS1) - kz .* real(dataS1) * 1im;
    dataS2z = kz .* imag(dataS2) - kz .* real(dataS2) * 1im;
    dataR1z = kz .* imag(dataR1) - kz .* real(dataR1) * 1im;
    dataR2z = kz .* imag(dataR2) - kz .* real(dataR2) * 1im;

    hZ  = real(ffti(datahz * M1M2, parameters));
    pz  = real(ffti(datapz * M1M2, parameters));
    qz  = real(ffti(dataqz * M1M2, parameters));
    s1z  = real(ffti(dataS1z * M1M2, parameters));
    s2z  = real(ffti(dataS2z * M1M2, parameters));
    r1z  = real(ffti(dataR1z * M1M2, parameters));
    r2z  = real(ffti(dataR2z * M1M2, parameters));

    # double x-derivatives
    datahxx = -kx2 .* real(datah) - kx2 .* imag(datah) * 1im;
    dataqxx = -kx2 .* real(dataq) - kx2 .* imag(dataq) * 1im;
    datapxx = -kx2 .* real(datap) - kx2 .* imag(datap) * 1im;

    # TODO: check accuracy of ffti here versus matlab
    # have seen small deviations
    hxx = real(ffti(datahxx * M1M2, parameters));
    pxx = real(ffti(datapxx * M1M2, parameters));
    qxx = real(ffti(dataqxx * M1M2, parameters));

    # double z-derivatives
    datahzz = -kz2 .* real(datah) - kz2 .* imag(datah) * 1im;
    dataqzz = -kz2 .* real(dataq) - kz2 .* imag(dataq) * 1im;
    datapzz = -kz2 .* real(datap) - kz2 .* imag(datap) * 1im;

    hzz = real(ffti(datahzz * M1M2, parameters));
    pzz = real(ffti(datapzz * M1M2, parameters));
    qzz = real(ffti(dataqzz * M1M2, parameters));

    # cross derivatives
    datahxz = -kxz .* real(datah) - kxz .* imag(datah) * 1im;
    dataqxz = -kxz .* real(dataq) - kxz .* imag(dataq) * 1im;
    datapxz = -kxz .* real(datap) - kxz .* imag(datap) * 1im;

    # transform variables back to real space
    hxz = real(ffti(datahxz * M1M2, parameters));
    pxz = real(ffti(datapxz * M1M2, parameters));
    qxz = real(ffti(dataqxz * M1M2, parameters));

    # gradient of the laplacian of h (=curvature K)
    datahLx = -(kx .* k2 .* imag(datah) - kx .* k2 .* real(datah) * 1im);
    datahLz = -(kz .* k2 .* imag(datah) - kz .* k2 .* real(datah) * 1im);

    dhdt = -real(dataqx) - real(datapz) - imag(dataqx) * 1im - imag(datapz) * 1im;

    h   = real(ffti(datah * M1M2, parameters));
    P   = real(ffti(datap * M1M2, parameters));
    q   = real(ffti(dataq * M1M2, parameters));
    s1   = real(ffti(dataS1 * M1M2, parameters));
    s2   = real(ffti(dataS2 * M1M2, parameters));
    r1   = real(ffti(dataR1 * M1M2, parameters));
    r2   = real(ffti(dataR2 * M1M2, parameters));

    hmax = maximum(h)

    # slightly increase upper bound compared to matlab
    # so that upperbound is included

    # cf.
    # julia: collect(0:0 .011178042476774046:1 .3413650972128854) -> 120 elements
    # matlab: 0:0 .011178042476774046:1 .3413650972128854 -> 121 elements
    y = reshape(collect(0:0.01 * hmax:hmax * 1.21), 1, :)
    y = y * ScaleH

    # determine yBar with datah in real space
    yBar = zeros(M1, M2, size(y, 2))
    for j = 1:1:M1
        for k = 1:1:M2
            yBar[j, k, :] = (y ./ h[j, k]);
        end
    end

    # calculate functions F0-F3 and G0-G3

    F0 = yBar - 1 / 2 * yBar.^2;
    F1 = yBar - 17 / 6 * yBar.^2 + 7 / 3 * yBar.^3 - 7 / 12 * yBar.^4;
    F2 = yBar - 13 / 2 * yBar.^2 + 57 / 4 * yBar.^3 - 111 / 8 * yBar.^4 + 99 / 16 * yBar.^5 - 33 / 32 * yBar.^6;
    F3 = yBar - 531 / 62 * yBar.^2 + 2871 / 124 * yBar.^3 - 6369 / 248 * yBar.^4 + 29601 / 2480 * yBar.^5 - 9867 / 4960 * yBar.^6;

    k, l = size(q)

    uField = similar(yBar)
    wField = similar(yBar)
    vField = similar(yBar)

    for k = 1:1:l
        F0Star = F0[:, k, :]
        F1Star = F1[:, k, :]
        F2Star = F2[:, k, :]
        j, m = size(y)
        for j = 1:1:m
            uField[:, k, j] = 3 * (q[:, k] - s1[:, k] - s2[:, k]) ./ h[:, k] .* F0Star[:, j] + 45 * s1[:, k] ./ h[:, k] .* F1Star[:, j] + 210 * s2[:, k] ./ h[:, k] .* F2Star[:, j]
            wField[:, k, j] = 3 * (P[:, k] - r1[:, k] - r2[:, k]) ./ h[:, k] .* F0Star[:, j] + 45 * r1[:, k] ./ h[:, k] .* F1Star[:, j] + 210 * r2[:, k] ./ h[:, k] .* F2Star[:, j]
            vField[:, k, j] = (y[:, j].^2 .* (-24 .* h[:, k].^6 .* (pz[:, k] + qx[:, k] + 14 .* r1z[:, k] + 69 .* r2z[:, k] + 14 .* s1x[:, k] + 69 .* s2x[:, k]) - 3465 .* (hZ[:, k] .* r2[:, k] + hx[:, k] .* s2[:, k]) .* y[:, j].^5 + 495 .* h[:, k] .* y[:, j].^4 .* (42 .* hZ[:, k] .* r2[:, k] + 42 .* hx[:, k] .* s2[:, k] + r2z[:, k] .* y[:, j] + s2x[:, k] .* y[:, j])
        - 105 .* h[:, k].^2 .* y[:, j].^3 .* (4 .* hZ[:, k] .* r1[:, k] + 444 .* hZ[:, k] .* r2[:, k] + 4 .* hx[:, k] .* s1[:, k] + 444 .* hx[:, k] .* s2[:, k] + 33 .* r2z[:, k] .* y[:, j] + 33 .* s2x[:, k] .* y[:, j])
        + 84 .* h[:, k].^3 .* y[:, j].^2 .* (20 .* hZ[:, k] .* r1[:, k] + 570 .* hZ[:, k] .* r2[:, k] + 20 .* hx[:, k] .* s1[:, k] + 570 .* hx[:, k] .* s2[:, k] + r1z[:, k] .* y[:, j] + 111 .* r2z[:, k] .* y[:, j] + s1x[:, k] .* y[:, j] + 111 .* s2x[:, k] .* y[:, j]) + 8 .* h[:, k].^5 .* (6 .* hZ[:, k] .* P[:, k] + 6 .* hx[:, k] .* q[:, k] + 84 .* hZ[:, k] .* r1[:, k] + 414 .* hZ[:, k] .* r2[:, k] + 84 .* hx[:, k] .* s1[:, k] + 414 .* hx[:, k] .* s2[:, k] + pz[:, k] .* y[:, j] + qx[:, k] .* y[:, j] + 84 .* r1z[:, k] .* y[:, j] + 909 .* r2z[:, k] .* y[:, j] + 84 .* s1x[:, k] .* y[:, j] + 909 .* s2x[:, k] .* y[:, j])
        - 6 .* h[:, k].^4 .* y[:, j] .* (4 .* hZ[:, k] .* P[:, k] + 4 .* hx[:, k] .* q[:, k] + 336 .* hZ[:, k] .* r1[:, k] + 3636 .* hZ[:, k] .* r2[:, k] + 336 .* hx[:, k] .* s1[:, k] + 3636 .* hx[:, k] .* s2[:, k] + 70 .* r1z[:, k] .* y[:, j] + 1995 .* r2z[:, k] .* y[:, j] + 70 .* s1x[:, k] .* y[:, j] + 1995 .* s2x[:, k] .* y[:, j]))) ./ (16 .* h[:, k].^8)

        end
    end

    uField = uField * ScaleU;
    vField = vField * ScaleU;
    wField = wField * ScaleU;
    h = h * ScaleH;

    k, l, m = size(uField)

    alpha = similar(yBar)
    for j = 1:1:M1
        for k = 1:1:M2
            for l = 1:1:m
                if y[l] > h[j, k]
                    uField[j, k, l] = NaN;
                    vField[j, k, l] = NaN;
                    wField[j, k, l] = NaN;
                    alpha[j, k, l] = 0;
                else
                    alpha[j, k, l] = 1;
                end
            end
        end
    end

    yBar = 1
    F0 = yBar - 1 / 2 * yBar.^2
    F1 = yBar - 17 / 6 .* yBar.^2 + 7 / 3 .* yBar.^3 - 7 / 12 .* yBar.^4
    F2 = yBar - 13 / 2 .* yBar.^2 + 57 / 4 * yBar.^3 - 111 / 8 .* yBar.^4 + 99 / 16 .* yBar.^5 - 33 / 32 .* yBar.^6
    F3 = yBar - 531 / 62 .* yBar.^2 + 2871 / 124 * yBar.^3 - 6369 / 248 .* yBar.^4 + 29601 / 2480 .* yBar.^5 - 9867 / 4960 .* yBar.^6

    ac, l = size(q)
    uSurf = zeros(ac, l)
    wSurf = similar(uSurf)
    vSurf = similar(uSurf)
    cWave = similar(uSurf)

    for k = 1:1:l
        F0Star = F0;
        F1Star = F1;
        F2Star = F2;

        j, m = size(y)

        j = 1

    # matlab -> julia
    # F0Star[:, j] -> F0Star
    # F1Star[:, j] -> F0Star
    # F2Star[:, j] -> F0Star
        uSurf[:,k] = 3 * (q[:, k] - s1[:, k] - s2[:, k]) ./ h[:, k] .* F0Star + 45 * s1[:, k] ./ h[:, k] .* F1Star + 210 * s2[:, k] ./ h[:, k] .* F2Star
        wSurf[:, k] = 3 * (P[:, k] - r1[:, k] - r2[:, k]) ./ h[:, k] .* F0Star + 45 * r1[:, k] ./ h[:, k] .* F1Star + 210 * r2[:, k] ./ h[:, k] .* F2Star
        vSurf[:,k] = (h[:,k].^2 .* (-24 .* h[:,k].^6 .* (pz[:,k] + qx[:,k] + 14 .* r1z[:,k] + 69 .* r2z[:,k] + 14 .* s1x[:,k] + 69 .* s2x[:,k]) - 3465 .* (hZ[:,k] .* r2[:,k] + hx[:,k] .* s2[:,k]) .* h[:,k].^5 + 495 .* h[:,k] .* h[:,k].^4 .* (42 .* hZ[:,k] .* r2[:,k] + 42 .* hx[:,k] .* s2[:,k] + r2z[:,k] .* h[:,k] + s2x[:,k] .* h[:,k])
    - 105 .* h[:,k].^2 .* h[:,k].^3 .* (4 .* hZ[:,k] .* r1[:,k] + 444 .* hZ[:,k] .* r2[:,k] + 4 .* hx[:,k] .* s1[:,k] + 444 .* hx[:,k] .* s2[:,k] + 33 .* r2z[:,k] .* h[:,k] + 33 .* s2x[:,k] .* h[:,k])
    + 84 .* h[:,k].^3 .* h[:,k].^2 .* (20 .* hZ[:,k] .* r1[:,k] + 570 .* hZ[:,k] .* r2[:,k] + 20 .* hx[:,k] .* s1[:,k] + 570 .* hx[:,k] .* s2[:,k] + r1z[:,k] .* h[:,k] + 111 .* r2z[:,k] .* h[:,k] + s1x[:,k] .* h[:,k] + 111 .* s2x[:,k] .* h[:,k])
    + 8 .* h[:,k].^5 .* (6 .* hZ[:,k] .* P[:,k] + 6 .* hx[:,k] .* q[:,k] + 84 .* hZ[:,k] .* r1[:,k] + 414 .* hZ[:,k] .* r2[:,k] + 84 .* hx[:,k] .* s1[:,k] + 414 .* hx[:,k] .* s2[:,k] + pz[:,k] .* h[:,k] + qx[:,k] .* h[:,k] + 84 .* r1z[:,k] .* h[:,k] + 909 .* r2z[:,k] .* h[:,k] + 84 .* s1x[:,k] .* h[:,k] + 909 .* s2x[:,k] .* h[:,k])
    - 6 .* h[:,k].^4 .* h[:,k] .* (4 .* hZ[:,k] .* P[:,k] + 4 .* hx[:,k] .* q[:,k] + 336 .* hZ[:,k] .* r1[:,k] + 3636 .* hZ[:,k] .* r2[:,k] + 336 .* hx[:,k] .* s1[:,k] + 3636 .* hx[:,k] .* s2[:,k] + 70 .* r1z[:,k] .* h[:,k] + 1995 .* r2z[:,k] .* h[:,k] + 70 .* s1x[:,k] .* h[:,k] + 1995 .* s2x[:,k] .* h[:,k]))) ./ (16 .* h[:,k].^8)

        for ls = 1:ac
            if (hx[ls, k] > 0.005)
                cWave[ls, k] = (hx[ls, k] .* uSurf[ls, k] - vSurf[ls, k]) ./ hx[ls, k];
            else
                cWave[ls, k] = NaN;
            end
        end
    end

    return uField, vField, wField, uSurf, vSurf, hx, y, cWave, h, alpha
end


text:
using Survival
using Test
using CSV
using DataFrames
using CategoricalArrays
using Distributions
using LinearAlgebra
using StatsBase
using StatsModels

@testset "Event times" begin
    @test isevent(EventTime{Int}(44, true))
    @test !isevent(EventTime(3.2, false))
    @test isevent(EventTime(2.5f0))

    @test !iscensored(EventTime(3))
    @test iscensored(EventTime(2.1, false))

    @test eltype(EventTime(1)) == Int

    @test sprint(show, EventTime(1, true)) == "1"
    @test sprint(show, EventTime(1, false)) == "1+"

    @test convert(Int, EventTime(1)) == 1
    @test convert(Int, EventTime(1, false)) == 1
    @test convert(Float64, EventTime(1)) == 1.0
    @test convert(EventTime, 1) == EventTime(1)

    @test isless(EventTime(1), EventTime(1, false))
    @test !isless(EventTime(2), EventTime(1, false))

    let x = [EventTime(2, false), EventTime(1), EventTime(2)]
        @test sort(x) == [EventTime(1), EventTime(2), EventTime(2, false)]
    end
end

@testset "Kaplan-Meier" begin
    t = [
        310, 361, 654, 728,  61,  81, 520, 473, 107, 122, 965, 731, 153, 433, 145,  95, 765,
        735,   5, 687, 345, 444,  60, 208, 821, 305, 226, 426, 705, 363, 167, 641, 740, 245,
        588, 166, 559, 450, 529, 351, 201, 524, 199, 550, 551, 543, 293, 511, 511, 371, 201,
         62, 356, 340, 315, 182, 364, 376, 384, 268, 266, 194, 348, 382, 296, 186, 145, 269,
        350, 272, 292, 332, 285, 243, 276,  79, 240, 202, 235, 224, 239, 173, 252,  92, 192,
        211, 175, 203, 105, 177,
    ]
    s = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
        1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1,
        0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0,
    ]

    # Results generated from R 3.3.2 using survival 2.41
    # survfit(Surv(t, s) ~ 1, conf.type = "log-log")
    # $surv
    r_surv = [
        0.9888889, 0.9777778, 0.9666667, 0.9555556, 0.9444444,  0.9333333,  0.9333333,  0.9220884,
        0.9220884, 0.9107045, 0.8993207, 0.8765531, 0.8651693,  0.8537855,  0.8424017,  0.8424017,
        0.8424017, 0.8424017, 0.8305369, 0.8186721, 0.8186721,  0.8066328,  0.7945935,  0.7705149,
        0.7705149, 0.7705149, 0.7580872, 0.7580872, 0.7580872,  0.7452383,  0.7452383,  0.7321639,
        0.7321639, 0.7321639, 0.7186054, 0.7186054, 0.7186054,  0.7045151,  0.7045151,  0.7045151,
        0.7045151, 0.6895254, 0.6895254, 0.6742026, 0.6742026,  0.6585235,  0.6428443,  0.6428443,
        0.6428443, 0.6263611, 0.609878,  0.5933948, 0.5769116,  0.5604284,  0.5604284,  0.5434457,
        0.526463,  0.526463,  0.5089143, 0.5089143, 0.5089143,  0.5089143,  0.4893406,  0.469767,
        0.4501934, 0.4306198, 0.4110461, 0.4110461, 0.3894121,  0.3677781,  0.3677781,  0.3677781,
        0.3432596, 0.3432596, 0.3432596, 0.3432596, 0.3120542,  0.2808487,  0.2496433,  0.2184379,
        0.1872325, 0.1560271, 0.1248217, 0.1248217, 0.08321444, 0.08321444, 0.08321444,
    ]
    # $n.risk
    r_risk = [
        90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68,
        67, 66, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45,
        44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23,
        22, 21, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10,  9,  8,  7,  6,  5,  4,  3,  2,  1
    ]
    # $std.err
    r_stderr = [
        0.01117336, 0.01589104, 0.01957401, 0.02273314, 0.02556550, 0.02817181, 0.02817181,
        0.03066888, 0.03066888, 0.03308929, 0.03539956, 0.03977328, 0.04186640, 0.04391166,
        0.04591747, 0.04591747, 0.04591747, 0.04591747, 0.04805852, 0.05016633, 0.05016633,
        0.05230824, 0.05442696, 0.05861552, 0.05861552, 0.05861552, 0.06082918, 0.06082918,
        0.06082918, 0.06318557, 0.06318557, 0.06561783, 0.06561783, 0.06561783, 0.06822833,
        0.06822833, 0.06822833, 0.07104408, 0.07104408, 0.07104408, 0.07104408, 0.07422800,
        0.07422800, 0.07755545, 0.07755545, 0.08104663, 0.08455340, 0.08455340, 0.08455340,
        0.08845361, 0.09238657, 0.09636405, 0.10039761, 0.10449888, 0.10449888, 0.10893570,
        0.11346828, 0.11346828, 0.11842498, 0.11842498, 0.11842498, 0.11842498, 0.12475150,
        0.13126159, 0.13798985, 0.14497408, 0.15225631, 0.15225631, 0.16157339, 0.17138826,
        0.17138826, 0.17138826, 0.18475887, 0.18475887, 0.18475887, 0.18475887, 0.20791044,
        0.23310483, 0.26120251, 0.29340057, 0.33150176, 0.37845310, 0.43957565, 0.43957565,
        0.59991117, 0.59991117, 0.59991117,
    ]
    # $lower
    r_lower = [
        0.92374348, 0.91406011, 0.90021666, 0.88590944, 0.87167275, 0.85762000, 0.85762000,
        0.84350521, 0.84350521, 0.82935113, 0.81542216, 0.78813835, 0.77474465, 0.76149619,
        0.74838092, 0.74838092, 0.74838092, 0.74838092, 0.73463920, 0.72104592, 0.72104592,
        0.70732360, 0.69373775, 0.66693587, 0.66693587, 0.66693587, 0.65314492, 0.65314492,
        0.65314492, 0.63887126, 0.63887126, 0.62441211, 0.62441211, 0.62441211, 0.60940174,
        0.60940174, 0.60940174, 0.59378857, 0.59378857, 0.59378857, 0.59378857, 0.57705901,
        0.57705901, 0.56006894, 0.56006894, 0.54279401, 0.52574922, 0.52574922, 0.52574922,
        0.50779307, 0.49008837, 0.47261930, 0.45537265, 0.43833741, 0.43833741, 0.42084761,
        0.40357982, 0.40357982, 0.38579388, 0.38579388, 0.38579388, 0.38579388, 0.36559209,
        0.34575866, 0.32627779, 0.30713724, 0.28832788, 0.28832788, 0.26728290, 0.24672627,
        0.24672627, 0.24672627, 0.22307271, 0.22307271, 0.22307271, 0.22307271, 0.19157693,
        0.16205660, 0.13440887, 0.10859921, 0.08465965, 0.06269905, 0.04293072, 0.04293072,
        0.01850525, 0.01850525, 0.01850525,
    ]
    # $upper
    r_upper = [
        0.9984273, 0.9943955, 0.9891262, 0.9830833, 0.9764926, 0.9694844, 0.9694844,
        0.9620777, 0.9620777, 0.9543175, 0.9463099, 0.9296781, 0.9211003, 0.9123712,
        0.9035043, 0.9035043, 0.9035043, 0.9035043, 0.8942180, 0.8848013, 0.8848013,
        0.8751547, 0.8653903, 0.8455370, 0.8455370, 0.8455370, 0.8352056, 0.8352056,
        0.8352056, 0.8244966, 0.8244966, 0.8135325, 0.8135325, 0.8135325, 0.8021438,
        0.8021438, 0.8021438, 0.7902942, 0.7902942, 0.7902942, 0.7902942, 0.7777438,
        0.7777438, 0.7648356, 0.7648356, 0.7515503, 0.7381102, 0.7381102, 0.7381102,
        0.7240036, 0.7097282, 0.6952904, 0.6806956, 0.6659482, 0.6659482, 0.6507130,
        0.6353162, 0.6353162, 0.6193746, 0.6193746, 0.6193746, 0.6193746, 0.6019235,
        0.5842205, 0.5662707, 0.5480775, 0.5296427, 0.5296427, 0.5096032, 0.4892180,
        0.4892180, 0.4892180, 0.4666889, 0.4666889, 0.4666889, 0.4666889, 0.4401063,
        0.4122125, 0.3830497, 0.3526092, 0.3208327, 0.2876052, 0.2527363, 0.2527363,
        0.2123639, 0.2123639, 0.2123639,
    ]

    km = fit(KaplanMeier, t, s)
    jl_surv = km.survival

    @test length(jl_surv) == length(r_surv)
    @test jl_surv ≈ r_surv atol=1e-6

    @test km.times == sort!(unique(t))
    @test km.natrisk == r_risk
    @test km.nevents == [sum(s[t .== tᵢ]) for tᵢ in sort!(unique(t))]
    @test km.ncensor == [sum(iszero, s[t .== tᵢ]) for tᵢ in sort!(unique(t))]
    @test km.stderr ≈ r_stderr atol=1e-6

    conf = confint(km)
    jl_lower = first.(conf)
    jl_upper = last.(conf)

    @test jl_lower ≈ r_lower atol=1e-6
    @test jl_upper ≈ r_upper atol=1e-6

    @test_throws DimensionMismatch fit(KaplanMeier, [1, 2], [true])
    @test_throws ArgumentError fit(KaplanMeier, Float64[], Bool[])

    km_et = fit(KaplanMeier, EventTime.(t, Bool.(s)))
    @test all(f->getfield(km, f) ≈ getfield(km_et, f), fieldnames(KaplanMeier))

    @test_throws ArgumentError fit(KaplanMeier, EventTime{Int}[])
end

@testset "Nelson-Aalen" begin
    t = [
        310, 361, 654, 728,  61,  81, 520, 473, 107, 122, 965, 731, 153, 433, 145,  95, 765,
        735,   5, 687, 345, 444,  60, 208, 821, 305, 226, 426, 705, 363, 167, 641, 740, 245,
        588, 166, 559, 450, 529, 351, 201, 524, 199, 550, 551, 543, 293, 511, 511, 371, 201,
         62, 356, 340, 315, 182, 364, 376, 384, 268, 266, 194, 348, 382, 296, 186, 145, 269,
        350, 272, 292, 332, 285, 243, 276,  79, 240, 202, 235, 224, 239, 173, 252,  92, 192,
        211, 175, 203, 105, 177,
    ]
    s = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
        1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1,
        0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0,
    ]

    na = fit(NelsonAalen, t, s)
    km = fit(KaplanMeier, t, s)

    @test na.times == km.times
    @test na.nevents == km.nevents
    @test na.ncensor == km.ncensor
    @test na.natrisk == km.natrisk
    @test exp.(-na.chaz[1:50]) ≈ km.survival[1:50] rtol=1e-2
    @test na.stderr[1:50] ≈ km.stderr[1:50] rtol=2e-2
    na_conf = confint(na)
    na_lower, na_upper = getindex.(na_conf, 1), getindex.(na_conf, 2)
    @test cdf.(Normal.(na.chaz, na.stderr), na_lower) ≈ fill(0.025, length(na.chaz)) rtol=1e-8
    @test cdf.(Normal.(na.chaz, na.stderr), na_upper) ≈ fill(0.975, length(na.chaz)) rtol=1e-8
    na_conf = confint(na, 0.01)
    na_lower, na_upper = getindex.(na_conf, 1), getindex.(na_conf, 2)
    @test cdf.(Normal.(na.chaz, na.stderr), na_lower) ≈ fill(0.005, length(na.chaz)) rtol=1e-8
    @test cdf.(Normal.(na.chaz, na.stderr), na_upper) ≈ fill(0.995, length(na.chaz)) rtol=1e-8
end

@testset "Cox" begin
    rossi = CSV.read(joinpath(@__DIR__, "data", "rossi.csv"), DataFrame)
    rossi.event = EventTime.(rossi.week, rossi.arrest .== 1)

    outcome = coxph(@formula(event ~ fin + age + race + wexp + mar + paro + prio), rossi; tol=1e-8)
    outcome_coefmat = coeftable(outcome)

    regressor_matrix = Matrix(rossi[!, [:fin, :age, :race, :wexp, :mar, :paro, :prio]])
    event_vector = rossi.event

    outcome_without_formula = coxph(regressor_matrix, event_vector)

    @test sprint(show, outcome_without_formula) == chomp("""
        CoxModel{Float64}

        Coefficients:
        ──────────────────────────────────────────────
              Estimate  Std.Error    z value  Pr(>|z|)
        ──────────────────────────────────────────────
        x1  -0.379416   0.191379   -1.98253     0.0474
        x2  -0.0574299  0.0219988  -2.61059     0.0090
        x3   0.31392    0.307995    1.01924     0.3081
        x4  -0.14981    0.212226   -0.705898    0.4803
        x5  -0.433724   0.38187    -1.13579     0.2560
        x6  -0.0848615  0.195756   -0.433505    0.6646
        x7   0.091521   0.0286469   3.1948      0.0014
        ──────────────────────────────────────────────
        """)

    coef_matrix = ModelMatrix(ModelFrame(@formula(event ~ 0 + fin + age + race + wexp + mar + paro + prio), rossi)).m
    outcome_from_matrix     = coxph(coef_matrix, rossi.event; tol=1e-8, l2_cost=0)
    outcome_from_matrix32   = coxph(Float32.(coef_matrix), rossi.event; tol=1e-5)
    outcome_from_matrix_int = coxph(Int64.(coef_matrix), rossi.event; tol=1e-6, l2_cost=0.0)

    expected_coefs = [
        -0.379422   0.191379   -1.98256   0.0474;
        -0.0574377  0.0219995  -2.61087   0.0090;
         0.3139     0.307993    1.01918   0.3081;
        -0.149796   0.212224   -0.705837  0.4803;
        -0.433704   0.381868   -1.13574   0.2561;
        -0.0848711  0.195757   -0.433554  0.6646;
         0.0914971  0.0286485   3.19378   0.0014
    ]

    @test coef(outcome_from_matrix) ≈ coef(outcome) atol=1e-5
    @test coef(outcome_from_matrix) ≈ coef(outcome_from_matrix32) atol=1e-4
    @test coef(outcome_from_matrix) ≈ coef(outcome_from_matrix_int) atol=1e-5
    @test nobs(outcome) == size(rossi, 1)
    @test dof(outcome) == 7
    @test loglikelihood(outcome) > nullloglikelihood(outcome)
    @test all(x->x > 0, eigen(outcome.model.fischer_info).values)
    @test outcome.model.fischer_info * vcov(outcome) ≈ I atol=1e-10
    @test norm(outcome.model.score) < 1e-5
    @test hcat(outcome_coefmat.cols[1:3]...) ≈ expected_coefs[:,1:3] atol=1e-5

    outcome_fin = coxph(@formula(event ~ fin), rossi; tol=1e-8)
    @test coeftable(outcome_fin).rownms == ["fin"]
    outcome_finrace = coxph(@formula(event ~ fin * race), rossi; tol=1e-8)
    @test coeftable(outcome_finrace).rownms == ["fin", "race","fin & race"]
    transform!(rossi, :fin => categorical, renamecols = false)
    outcome_fincat = coxph(@formula(event ~ fin), rossi; tol=1e-8)
    @test coeftable(outcome_fincat).rownms == ["fin: 1"]
    @test coef(outcome_fin) ≈ coef(outcome_fincat) atol=1e-8
    outcome_fincatrace = coxph(@formula(event ~ fin * race), rossi; tol=1e-8)
    @test coeftable(outcome_fincatrace).rownms == ["fin: 1", "race","fin: 1 & race"]
    @test coef(outcome_fincatrace) ≈ coef(outcome_finrace) atol=1e-8
    transform!(rossi, :race => categorical, renamecols = false)
    outcome_fincatracecat = coxph(@formula(event ~ fin * race), rossi; tol=1e-8)
    @test coeftable(outcome_fincatracecat).rownms == ["fin: 1", "race: 1","fin: 1 & race: 1"]
    @test coef(outcome_fincatracecat) ≈ coef(outcome_finrace) atol=1e-8
end

@testset "Newton-Raphson" begin
    function fgh!(x, grad, hes, compute_ders)
        if compute_ders
            grad[1] = 2exp(x[1])*(exp(x[1]) - 1)
            hes[1,1] = 2exp(x[1])*(2exp(x[1]) - 1)
        end
        (exp(x[1]) - 1)^2
    end
    x, y, grad, hes = Survival.newton_raphson(fgh!, [2.2], tol=1e-5)
    @test x ≈ [0] atol=1e-5
    @test y ≈ 0 atol=1e-5
    @test grad ≈ [0] atol=1e-5
    @test hes ≈ [2] atol=1e-5
    @test_throws ConvergenceException Survival.newton_raphson(fgh!, [2.2], max_iter=2)

    function wrong_fgh!(x, grad, hes, compute_ders)
        if compute_ders
            grad[1] = -2exp(x[1])*(exp(x[1]) - 1) # wrong sign
            hes[1,1] = 2exp(x[1])*(2exp(x[1]) - 1)
        end
        (exp(x[1]) - 1)^2
    end
    @test_throws ErrorException Survival.newton_raphson(wrong_fgh!, [2.2])
end


text:
# # API

# ## Warnings

"""
    warn_on_recipe_aliases!(plt, plotattributes, recipe_type, args...)

Warn if an alias is dedected in `plotattributes` after a recipe of type `recipe_type` is
applied to 'args'. `recipe_type` is either `:user`, `:type`, `:plot` or `:series`.
"""
function warn_on_recipe_aliases!(plt, plotattributes, recipe_type, args...) end

# ## Grouping

"""
    splittable_attribute(plt, key, val, len)

Returns `true` if the attribute `key` with the value `val` can be split into groups with
group provided as a vector of length `len`, `false` otherwise.
"""
splittable_attribute(plt, key, val, len) = false
splittable_attribute(plt, key, val::AbstractArray, len) =
    !(key in (:group, :color_palette)) && length(axes(val, 1)) == len
splittable_attribute(plt, key, val::Tuple, len) =
    all(v -> splittable_attribute(plt, key, v, len), val)

"""
    split_attribute(plt, key, val, indices)

Select the proper indices from `val` for attribute `key`.
"""
split_attribute(plt, key, val::AbstractArray, indices) =
    val[indices, fill(Colon(), ndims(val) - 1)...]
split_attribute(plt, key, val::Tuple, indices) =
    Tuple(split_attribute(plt, key, v, indices) for v in val)

# ## Preprocessing attributes

"""
    preprocess_attributes!(plt, plotattributes)

Any plotting package specific preprocessing of user or recipe input happens here.
For example, Plots replaces aliases and expands magic arguments.
"""
function preprocess_attributes!(plt, plotattributes) end

# TODO: should the Plots version be defined as fallback in RecipesPipeline?
"""
    is_subplot_attribute(plt, attr)

Returns `true` if `attr` is a subplot attribute, otherwise `false`.
"""
is_subplot_attribute(plt, attr) = false

# TODO: should the Plots version be defined as fallback in RecipesPipeline?
"""
    is_axis_attribute(plt, attr)

Returns `true` if `attr` is an axis attribute, i.e. it applies to `xattr`, `yattr` and
`zattr`, otherwise `false`.
"""
is_axis_attribute(plt, attr) = false

# ## User recipes

"""
    process_userrecipe!(plt, attributes_list, attributes)

Do plotting package specific post-processing and add series attributes to attributes_list.
For example, Plots increases the number of series in `plt`, sets `:series_plotindex` in
attributes and possible adds new series attributes for errorbars or smooth.
"""
function process_userrecipe!(plt, attributes_list, attributes)
    push!(attributes_list, attributes)
end

"""
    get_axis_limits(plt, letter)

Get the limits for the axis specified by `letter` (`:x`, `:y` or `:z`) in `plt`. If it
errors, `tryrange` from PlotUtils is used.
"""
get_axis_limits(plt, letter) = throw(ErrorException("Axis limits not defined."))

# ## Plot recipes

"""
    type_alias(plt, st)

Return the seriestype alias for `st`.
"""
type_alias(plt, st) = st

# ## Plot setup

"""
    plot_setup!(plt, plotattributes, kw_list)

Setup plot, subplots and layouts.
For example, Plots creates the backend figure, initializes subplots, expands extrema and
links subplot axes.
"""
function plot_setup!(plt, plotattributes, kw_list) end

# ## Series recipes

"""
    slice_series_attributes!(plt, kw_list, kw)

For attributes given as vector with one element per series, only select the value for
current series.
"""
function slice_series_attributes!(plt, kw_list, kw) end

"""
    process_sliced_series_attributes!(plt, kw_list)

All series attributes are now properly resolved. Any change of the `kw_list` before the application of recipes must come here.
"""
function process_sliced_series_attributes!(plt, kw_list) end

"""
    series_defaults(plt)

Returns a `Dict` storing the defaults for series attributes.
"""
series_defaults(plt) = Dict{Symbol, Any}()

# TODO: Add a more sensible fallback including e.g. path, scatter, ...

"""
    is_seriestype_supported(plt, st)

Check if the plotting package natively supports the seriestype `st`.
"""
is_seriestype_supported(plt, st) = false

"""
    is_key_supported(key)

Check if the plotting package natively supports the attribute `key`
"""
RecipesBase.is_key_supported(key) = true

# ## Finalizer

"""
    add_series!(plt, kw)

Adds the series defined by `kw` to the plot object.
For example Plots updates the current subplot arguments, expands extrema and pushes the
the series to the series_list of `plt`.
"""
function add_series!(plt, kw) end


text:
#=
$ mkdir ~/tmp
$ cd ~/tmp
$ git clone https://github.com/avik-pal/Flux.jl.git
$ cd Flux.jl
$ git checkout -b depthwiseconv origin depthwiseconv
replace ~/.julia/packages/Flux/UHjNa with all of them under ~/tmp/Flux.jl
=#

using Flux
using Flux:Chain
using Flux:@treelike

struct NaiveBottleNeck
    layers::Chain
    stride::Int
end

function NaiveBottleNeck(in_ch::Int, middle_ch::Int, out_ch::Int; stride=1)
    chain = Chain(Conv((1,1), in_ch=>middle_ch, relu),
                  Conv((3,3), middle_ch=>middle_ch, relu, stride=stride, pad=1),
                  Conv((1,1), middle_ch=>out_ch),
                  )
    NaiveBottleNeck(chain, stride)
end

@treelike NaiveBottleNeck

function (m::NaiveBottleNeck)(x)
    h = m.layers(x)
    if size(h)==size(x)
        relu.(h + x)
    else
        relu.(h)
    end
end

#= test
out=NaiveRes(5=>5)(rand(23,30,5,2))
println(size(out))
=#

struct NaiveResNet
    layers::Chain
end

forward() = Chain(Conv((7,7), 3=>32, relu, stride=2, pad=3),
                    NaiveBottleNeck(32, 16, 64),
                    NaiveBottleNeck(64,32,64,stride=2),
                    NaiveBottleNeck(64,32,128),
                    NaiveBottleNeck(128,64,128,stride=2),
                    NaiveBottleNeck(128,64,256),
                    NaiveBottleNeck(256,128,256,stride=2),
                    NaiveBottleNeck(256,128,256),
                    MaxPool((2,2)),
                    x -> reshape(x, :, size(x, 4)),
                    Dense(12544, 2048),
                    Dense(2048, 1024),
                    Dense(1024, 101),
                    softmax
                    )

NaiveResNet() = NaiveResNet(forward())

@treelike NaiveResNet

(model::NaiveResNet)(x) = model.layers(x)

using Flux
using CuArrays

function test_func()
    model = NaiveResNet()
    @show size(model(rand(224, 224, 3, 2)))
    model = model |> gpu
    @show size(model(rand(224, 224, 3, 32)|> gpu))
    batchsize=64
    loss(x,y)= Flux.crossentropy(model(x), y)
    total_loss=0
    for i in 1:500
        X = rand(224,224,3,batchsize) |> gpu
        Y = Flux.onehotbatch([rand(1:101) for i in 1:batchsize], 1:101) |> gpu
        Flux.testmode!(model)
        total_loss += loss(X, Y).data
    end
end

#test_func()


text:
module MakieLayout

using AbstractPlotting
using AbstractPlotting: Rect2D
import AbstractPlotting: IRect2D
using AbstractPlotting.Keyboard
using AbstractPlotting.Mouse
using AbstractPlotting: ispressed, is_mouseinside
using Observables: onany
import Observables
import Formatting
using Match
import Animations
import PlotUtils
using GridLayoutBase
import Showoff

include("types.jl")
include("helpers.jl")
include("mousestatemachine.jl")
include("ticklocators/linear.jl")
include("ticklocators/wilkinson.jl")
include("defaultattributes.jl")
include("lineaxis.jl")
include("lobjects/laxis.jl")
include("lobjects/lcolorbar.jl")
include("lobjects/ltext.jl")
include("lobjects/lslider.jl")
include("lobjects/lbutton.jl")
include("lobjects/lrect.jl")
include("lobjects/ltoggle.jl")
include("lobjects/llegend.jl")
include("lobjects/lobject.jl")
include("lobjects/lscene.jl")

export LAxis
export LSlider
export LButton
export LColorbar
export LText
export LRect
export LToggle
export LLegend
export LegendEntry, MarkerElement, PolyElement, LineElement, LegendElement
export LScene
export linkxaxes!, linkyaxes!, linkaxes!
export AxisAspect, DataAspect
export autolimits!
export AutoLinearTicks, ManualTicks, CustomTicks, WilkinsonTicks
export layoutscene
export set_close_to!
# TODO: These functions are exported but are not in the docs and have no docstrings
export hidexdecorations!, hideydecorations!, hidedecorations!
export tight_xticklabel_spacing!, tight_yticklabel_spacing!, tight_ticklabel_spacing!, tightlimits!

# from GridLayoutBase
export GridLayout, GridPosition
export GridLayoutSpec
export BBox
export LayoutObservables
export Inside, Outside, Mixed
export Fixed, Auto, Relative, Aspect
export width, height, top, bottom, left, right
export with_updates_suspended
export appendcols!, appendrows!, prependcols!, prependrows!, deletecol!, deleterow!, trim!
export gridnest!
export AxisAspect, DataAspect
export colsize!, rowsize!, colgap!, rowgap!
export Left, Right, Top, Bottom, TopLeft, BottomLeft, TopRight, BottomRight
export grid!, hbox!, vbox!
export swap!
export ncols, nrows

const FPS = Node(30)
const COLOR_ACCENT = Ref(RGBf0(((79, 122, 214) ./ 255)...))
const COLOR_ACCENT_DIMMED = Ref(RGBf0(((174, 192, 230) ./ 255)...))

end # module


text:
struct ConstFactor <: AbstractFactor{UInt8}
    name::String
    level::String
    length::Integer
end

Base.length(factor::ConstFactor) = factor.length

getlevels(factor::ConstFactor) = [factor.level]

function ConstFactor(length::Integer)
    ConstFactor("Intercept", "Intercept", length)
end

function slice(factor::ConstFactor, fromobs::Integer, toobs::Integer, slicelength::Integer)
    if fromobs > toobs
        EmptySeq{AbstractVector{UInt8}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(factor))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        buffer = ones(UInt8, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{UInt8,1,Array{UInt8,1},Tuple{UnitRange{Int64}},true}) do rng
            view(buffer, 1:(rng[2] - rng[1] + 1))
        end
    end
end

function Base.map(factor::ConstFactor, dataframe::AbstractDataFrame)
    ConstFactor(factor.name, factor.level, length(dataframe))
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

const alg = MethodOfSteps(BS3())
const prob = DDEProblemLibrary.prob_dde_constant_1delay_ip

const dde_int = init(prob, alg)
const sol = solve!(dde_int)

@testset "reference" begin
  @test sol.errors[:l∞] < 3.0e-5
  @test sol.errors[:final] < 2.1e-5
  @test sol.errors[:l2] < 1.2e-5
end

@testset "constant lag function" begin
  # constant delay specified as lag function
  prob2 = remake(prob; constant_lags = nothing, dependent_lags = ((u, p, t) -> 1,))
  dde_int2 = init(prob2, alg)
  sol2 = solve!(dde_int2)

  @test getfield.(dde_int.tracked_discontinuities, :t) ≈
    getfield.(dde_int2.tracked_discontinuities, :t)
  @test getfield.(dde_int.tracked_discontinuities, :order) ==
    getfield.(dde_int2.tracked_discontinuities, :order)

  # worse than results above with constant delays specified as scalars
  @test sol2.errors[:l∞] < 3.2e-3
  @test sol2.errors[:final] < 1.2e-4
  @test sol2.errors[:l2] < 1.2e-3

  # simple convergence tests
  @testset "convergence" begin
    sol3 = solve(prob2, alg; abstol=1e-9, reltol=1e-6)

    @test sol3.errors[:l∞] < 3.0e-6
    @test sol3.errors[:final] < 1.4e-7
    @test sol3.errors[:l2] < 8.4e-7

    sol4 = solve(prob2, alg; abstol=1e-13, reltol=1e-13)

    @test sol4.errors[:l∞] < 6.0e-11
    @test sol4.errors[:final] < 4.7e-11
    @test sol4.errors[:l2] < 5.9e-12
  end
end

# without any delays specified is worse
@testset "without delays" begin
  prob2 = remake(prob; constant_lags = nothing)
  sol2 = solve(prob2, alg)

  @test sol2.errors[:l∞] > 1.0e-2
  @test sol2.errors[:final] > 6.3e-6
  @test sol2.errors[:l2] > 4.6e-3
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Sources
export Source, CompositeSource

using ....OpticSim
using ...Emitters
using ...Geometry
using ..Spectrum
using ..Directions
using ..Origins
using ..AngularPower
using LinearAlgebra
using Distributions

abstract type AbstractSource{T} <: OpticSim.OpticalRayGenerator{T} end

Base.getindex(s::AbstractSource, index) = generate(s, index)
Base.firstindex(s::AbstractSource) = 0
Base.lastindex(s::AbstractSource) = length(s) - 1
Base.copy(a::AbstractSource) = a # most don't have any heap allocated stuff so don't really need copying

"""
    Source{T<:Real, Tr<:Transform{T}, S<:Spectrum.AbstractSpectrum{T}, O<:Origins.AbstractOriginDistribution{T}, D<:Directions.AbstractDirectionDistribution{T}, P<:AngularPower.AbstractAngularPowerDistribution{T}} <: AbstractSource{T}

This data-type represents the basic emitter (Source), which is a combination of a Spectrum, Angular Power Distribution, Origins and Directions distibution and a 3D Transform.

```julia
Source(::Type{T} = Float64;
       transform::Tr = Transform(),
       spectrum::S = Spectrum.Uniform(),
       origins::O = Origins.Point(),
       directions::D = Directions.Constant(),
       power::P = AngularPower.Lambertian(),
       sourcenum::Integer = 0) where {
            Tr<:Transform,
            S<:Spectrum.AbstractSpectrum,
            O<:Origins.AbstractOriginDistribution,
            D<:Directions.AbstractDirectionDistribution,
            P<:AngularPower.AbstractAngularPowerDistribution,
            T<:Real}

Source(transform::Tr, spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64; sourcenum::Integer = 0) where {
            Tr<:Transform,
            S<:Spectrum.AbstractSpectrum,
            O<:Origins.AbstractOriginDistribution,
            D<:Directions.AbstractDirectionDistribution,
            P<:AngularPower.AbstractAngularPowerDistribution,
            T<:Real}
```
"""
struct Source{
    T<:Real,
    Tr<:Transform{T},
    S<:Spectrum.AbstractSpectrum{T},
    O<:Origins.AbstractOriginDistribution{T},
    D<:Directions.AbstractDirectionDistribution{T},
    P<:AngularPower.AbstractAngularPowerDistribution{T}
} <: AbstractSource{T}
    transform::Tr
    spectrum::S
    origins::O
    directions::D
    power_distribution::P
    sourcenum::Integer

    function Source(
        ::Type{T} = Float64;
        transform::Tr = Transform(T),
        spectrum::S = Spectrum.Uniform(T),
        origins::O = Origins.Point(T),
        directions::D = Directions.Constant(T),
        power::P = AngularPower.Lambertian(T),
        sourcenum::Integer = 0
    ) where {
        Tr<:Transform,
        S<:Spectrum.AbstractSpectrum,
        O<:Origins.AbstractOriginDistribution,
        D<:Directions.AbstractDirectionDistribution,
        P<:AngularPower.AbstractAngularPowerDistribution,
        T<:Real
    }
        new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum)
    end

    function Source(
        transform::Tr,
        spectrum::S,
        origins::O,
        directions::D,
        power::P,
        ::Type{T} = Float64;
        sourcenum::Integer = 0
    ) where {
        Tr<:Transform,
        S<:Spectrum.AbstractSpectrum,
        O<:Origins.AbstractOriginDistribution,
        D<:Directions.AbstractDirectionDistribution,
        P<:AngularPower.AbstractAngularPowerDistribution,
        T<:Real
    }
         new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum)
    end
end

# used to not generate new origin points if we can use last points - mainly to keep consistency of origin generation when randomness is involved
struct SourceGenerationState{T<:Real}
    n::Integer
    last_index_used::Integer
    last_point_generated::Vec3{T}
end

Base.length(s::Source) = length(s.origins) * length(s.directions)
Base.iterate(s::Source{T}) where {T<:Real} = iterate(s, SourceGenerationState(1, -1, zeros(Vec3{T})))
Base.iterate(s::Source, state::SourceGenerationState) = state.n > length(s) ? nothing : generate(s, state)

function Emitters.generate(s::Source{T}, n::Integer) where {T<:Real}
    return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1]
end

function Emitters.generate(s::Source{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real}
    # @info "New State Generation"
    n::Integer = state.n - 1
    origin_index = floor(Integer, (n / length(s.directions)))
    direction_index = mod(n, length(s.directions))

    if (origin_index == state.last_index_used)
        o = state.last_point_generated
    else
        o = generate(s.origins, origin_index)
        o = s.transform * o
    end

    # o = generate(s.origins, origin_index)
    d = generate(s.directions, direction_index)

    # rotate the direction according to the orientation of the ray-originator
    d = rotation(s.transform) * d

    ray = OpticSim.Ray(o, d)
    power, wavelength = generate(s.spectrum)

    power = apply(s.power_distribution, s.transform, power, ray)

    # return (OpticSim.Ray(o, d), SourceGenerationState(state.n+1, origin_index, o))
    return (OpticSim.OpticalRay(ray, power, wavelength; sourcenum=s.sourcenum), SourceGenerationState(state.n+1, origin_index, o))
end

"""
    CompositeSource{T} <: AbstractSource{T}

This data-type represents the composite emitter (Source) which is constructed with a list of basic or composite emitters and a 3D Transform.

```julia
CompositeSource(transform::Transform{T}, sources::Vector{<:AbstractSource}) where {T<:Real}
```
"""
struct CompositeSource{T} <: AbstractSource{T}
    transform::Transform{T}
    sources::Vector{<:AbstractSource}

    uniform_length::Int
    total_length::Int
    start_indexes::Vector{Int}

    function CompositeSource(transform::Transform{T}, sources::Vector{<:AbstractSource}) where {T<:Real}
        lens = [length(src) for src in sources]
        size1 = findmin(lens)[1]
        size2 = findmax(lens)[1]
        if (size1 == size2)
            uniform_length = size1
            start_indexes = Vector{Int}()
            total_length = length(sources) * uniform_length
            # @info "Uniform Length $size1 total=$total_length"
        else
            uniform_length = -1

            start_indexes = Vector{Int}()
            start = 0
            for s in sources
                push!(start_indexes, start)
                start = start + length(s)
            end
            push!(start_indexes, start)
            # @info start_indexes
            total_length = start
        end

        new{T}(transform, sources, uniform_length, total_length, start_indexes)
    end
end

Base.length(s::CompositeSource) = s.total_length
Base.iterate(s::CompositeSource{T}) where {T<:Real} = iterate(s, SourceGenerationState(1, -1, zeros(Vec3{T})))
Base.iterate(s::CompositeSource, state::SourceGenerationState) = state.n > length(s) ? nothing : generate(s, state)

function Emitters.generate(s::CompositeSource{T}, n::Integer) where {T<:Real}
    return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1]
end

function Emitters.generate(s::CompositeSource{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real}
    # @info "New Composite State Generation"
    n::Integer = state.n - 1

    if (s.uniform_length != -1)
        source_index = floor(Integer, (n / s.uniform_length))
        ray_index = mod(n, s.uniform_length)
        # @info "New Composite State Generation: source=$source_index  ray=$ray_index  ($(length(s.sources)))"
    else
        n = mod(n, s.total_length)

        ## perform a binary search to find out which source can generate the spesific ray
        low = 1
        high = length(s.start_indexes) - 1    # last cell in this vector is a dummy cell containing the index after the last ray
        source_index = -1
        while (low <= high)
            mid = floor(Integer, (low + high) / 2)

            if (n < s.start_indexes[mid])
                high = mid - 1
            elseif (n >= s.start_indexes[mid+1])
                low = mid + 1
            else
                source_index = mid
                break;
            end
        end
        @assert source_index != -1
        # @info "Bin Search n=$n source_index=$source_index"
        source_index = source_index - 1
        ray_index = n - s.start_indexes[source_index+1]
    end

    ray, new_state = generate(s.sources[source_index+1], SourceGenerationState(ray_index+1, state.last_index_used, state.last_point_generated))
    ray = s.transform * ray
    return (ray, SourceGenerationState(state.n+1, new_state.last_index_used, new_state.last_point_generated))
end

end # module Sources


text:
function juliatype(x::AbstractVector)
    Union{map(juliatype, x)...}
end

# iterate a Julia Union{...} type, producing an array of unioned types
Base.@pure function eachunion(T, elems)
    if elems === nothing
        elems = Type[]
    end
    if T isa Union
        if T.a isa Union
            push!(elems, T.b)
            return eachunion(T.a, elems)
        else
            push!(elems, T.a)
            return eachunion(T.b, elems)
        end
    else
        push!(elems, T)
        return Tuple{elems...}
    end
end

function schematype(::StructTypes.Struct, U::Union)
    T = eachunion(U, nothing)
    sch = Vector{Union{String, SchemaType, LogicalType}}(undef, fieldcount(T))
    for i = 1:fieldcount(T)
        sch[i] = schematype(fieldtype(T, i))
    end
    return sch
end

function unionindex(T, UT)
    for i = 1:length(UT)
        T == UT[i] && return (i - 1)
    end
    error("value with avro schema $T not valid for union type $UT")
end

function writevalue(B::Binary, UT::UnionType, x::T, buf, pos, len, opts) where {T}
    sch = schematype(T)
    i = unionindex(sch, UT)
    pos = writevalue(B, long, i, buf, pos, len, opts)
    return writevalue(B, sch, x, buf, pos, len, opts)
end

function nbytes(UT::UnionType, x::T) where {T}
    sch = schematype(T)
    i = unionindex(sch, UT)
    n = nbytes(long, i)
    return n + nbytes(sch, x)
end

function readvalue(B::Binary, UT::UnionType, ::Type{T}, buf, pos, len, opts) where {T}
    i, pos = readvalue(B, long, Int, buf, pos, len, opts)
    sch = UT[i + 1]
    S = fieldtype(eachunion(T, nothing), i + 1)
    x, pos = readvalue(B, sch, S, buf, pos, len, opts)
    return x, pos
end

function skipvalue(B::Binary, UT::UnionType, ::Type{T}, buf, pos, len, opts) where {T}
    i, pos = readvalue(B, long, Int, buf, pos, len, opts)
    sch = UT[i + 1]
    S = fieldtype(eachunion(T, nothing), i + 1)
    pos = skipvalue(B, sch, S, buf, pos, len, opts)
    return pos
end