text:
import Base.length

abstract PDistribution

check_pd_sum!{N <: Real}(_ps::Array{N,1}) = !isapprox(sum(_ps), 1) ? throw("Probability Distributions must sum to 1 (got $(sum(_ps)))") : true;
check_pd_sum!{N <: Real}(_ps::Array{N,2}) = !isapprox(sum(_ps), size(_ps,1)) ? throw("Conditional Probability Distributions must consist of probability distributions which sum to 1") : true;
check_unique!(s1) = !(length(s1) == length(unique(s1))) ? throw("States must be unique") : true;

type ProbabilityDistributionVector{T <: Real, K} <: PDistribution
    states::Array{K,1}
    ps::Array{T,1}

    function ProbabilityDistributionVector(_states::Array{K,1}, _ps::Array{T,1})
        check_unique!(_states)
        length(_ps) > 0 ? check_pd_sum!(_ps) : nothing;
        new(_states, _ps)
    end
end

type ProbabilityDistributionMatrix{T <: Real, K, Q} <: PDistribution
    states::Array{K,1}
    ps::Array{T,2}
    conditionals::Array{Q,1}

    function ProbabilityDistributionMatrix(_states::Array{K,1}, _ps::Array{T,2}, _conditionals::Array{Q,1})
        check_unique!(_states)
        check_pd_sum!(_ps)
        new(_states, _ps, _conditionals)
    end
end

states(pd::PDistribution) = pd.states
conditionals(pd::ProbabilityDistributionMatrix) = pd.conditionals

unsorted_equality(arr1::Array, arr2::Array) = length(setdiff(Set(arr1),Set(arr2))) == 0

==(pd1::PDistribution, pd2::PDistribution) = unsorted_equality(states(pd1),states(pd2))

ProbabilityDistribution{T,K}(states::Array{K,1}, ps::Array{T,1}) = ProbabilityDistributionVector{T,K}(states, ps)
ProbabilityDistribution{T,K,Q}(states::Array{K,1}, ps::Array{T,2}, conditionals::Array{Q,1}) = ProbabilityDistributionMatrix{T,K,Q}(states, ps, conditionals)
ProbabilityDistribution{K}(states::Array{K,1}) = ProbabilityDistribution(states, Real[])

function probabilities(pd::PDistribution)
    pd.ps
end

function length(pd::PDistribution)
    length(pd.ps)
end

function getindex{R,K}(pd::ProbabilityDistributionVector{R,K}, key::K)
    !(key in states(pd)) ? throw("Invalid key") : nothing

    pd.ps[key .== states(pd)]
end

function getindex{R,K,Q}(pd::ProbabilityDistributionMatrix{R,K,Q}, key::K)
    s = states(pd)
    !(key in s) ? throw("Invalid key") : nothing

    pd.ps[:, key .== states(pd)]
end

function getindex{R,K,Q}(pd::ProbabilityDistributionMatrix{R,K,Q}, keys::Array{K,1})
    s = states(pd)
    for key in keys
        !(key in s) ? throw("Invalid key") : nothing
    end
    indencies = map(x -> findfirst(x,s), keys)

    pd.ps[:, indencies]
end

function getindex{R,K}(pd::ProbabilityDistributionVector{R,K}, keys::Array{K,1})
    s = states(pd)
    for key in keys
        !(key in s) ? throw("Invalid key") : nothing
    end
    indencies = map(x -> findfirst(s,x), keys)

    pd.ps[indencies]
end

## Later...

## function getindex{T,K}(pd::ProbabilityDistribution{T,K}, key1::K, key2::K)
##     if !(key1 in states(pd)[1])
##         throw("Invalid key")
##     end

##     if !(key2 in states(pd)[2])
##         throw("Invalid key")
##     end

##     pd.ps[key1 .== states(pd)[1], key2 .== states(pd)[2]][1]
## end

function lpd{K}(states::K)
    m = length(states)
    pb = Array(Float64, m)
    fill!(pb, 1/m)
    ProbabilityDistribution(pb, states)
end

function unique_states(p1::PDistribution, p2::PDistribution)
    length(setdiff(Set(states(p1)), Set(states(p2)))) == length(states(p1))
end

function *(p1::ProbabilityDistributionVector, p2::ProbabilityDistributionVector)
    if !unique_states(p1,p2)
        throw("State conflict")
    end

    ps_m=probabilities(p1)'.*probabilities(p2)
    s_p1=states(p1)
    s_p2=states(p2)

    combinations=length(s_p1)*length(s_p2)
    # indexes of values to in the ps_m
    ## indencies=find(triu(ones(Bool, size(ps_m,1), size(ps_m,2))))

    common_type=promote_type(eltype(s_p1), eltype(s_p2))
    j_states=Array(Set{common_type}, combinations)

    size_ps_m=size(ps_m)

    for i=1:combinations
        p1_s_index, p2_s_index = single_to_multi_index(i, size_ps_m)
        j_states[i]=Set{common_type}({s_p1[p1_s_index], s_p2[p2_s_index]})
    end

    ProbabilityDistribution(j_states, reshape(ps_m, combinations))
end

function join_sets{K}(x::Set{K}, y)
    c_type=promote_type(K,typeof(y))
    Set{c_type}(c_type[collect(x),y])
end

function join_sets(x, y)
    c_type=promote_type(typeof(x),typeof(y))
    Set{c_type}(c_type[x,y])
end

function zip_sets(arr1::Array, arr2::Array)
    if length(arr1) != length(arr2)
        throw("Length of given arrays must be the same.")
    end
    [ join_sets(arr1[i], arr2[i]) for i=1:length(arr1) ]
end

function *(p1::ProbabilityDistributionMatrix, p2::ProbabilityDistributionVector)
    c_p1 = conditionals(p1)
    s_p2 = states(p2)

    !unsorted_equality(c_p1, s_p2) ? throw("State mismatch") : nothing;

    p1_order=sortperm(c_p1)
    p2_order=sortperm(s_p2)

    sorted_ps_m=probabilities(p1)[p1_order,:] .* probabilities(p2)[p2_order]
    sorted_ps=reshape(sorted_ps_m[p1_order,:], length(sorted_ps_m))

    l_p2_elems = reduce(vcat, map(x -> fill(x, length(c_p1)), states(p1)))
    l_p1_sets = reduce(vcat, map(x -> c_p1, [1:length(states(p1))]))

    ProbabilityDistribution(zip_sets(l_p1_sets, l_p2_elems), sorted_ps)
end

setify(x) = typeof(x) <: Set ? x : Set{typeof(x)}({x});
desetify(arr::Array) = map(desetify, arr);
desetify(x) = typeof(x) <: Set && length(x) == 1 ? collect(x)[1] : x;

function /(p1::ProbabilityDistributionVector, p2::ProbabilityDistributionVector)
    s_p1 = map(setify, states(p1))
    s_p2 = map(setify, states(p2))
    p_p1 = probabilities(p1)
    p_p2 = probabilities(p2)

    unique_s_p2 = reduce(union, s_p2)
    raw_states = map(x -> copy(setdiff(x, unique_s_p2)), s_p1)
    _states = unique(raw_states)
    _u_s = reduce(union, _states)

    c_p1 = map(x -> copy(setdiff(x, _u_s)), s_p1)

    _conditionals = unique(c_p1)

    indencies=Dict(_states, fill(Integer[],length(_states)))
    map(x -> indencies[x[2]] = [indencies[x[2]], x[1]], enumerate(raw_states))

    _ps = zeros(FloatingPoint, length(_conditionals), length(_states))

    for (i,s) in enumerate(_states)
        idxs = indencies[s]
        p_s_p1 = p_p1[idxs][order_perm(c_p1[idxs], _conditionals)]
        p_s_p2 = p2[desetify(_conditionals)]

        _ps[:,i] = p_s_p1 ./ p_s_p2
    end

    ProbabilityDistribution(_states, _ps, _conditionals)
end


text:
module NewtonsMethod

greet() = print("Hello World!")

using LinearAlgebra, ForwardDiff

function fixedpointmap(f, f_prime, x_0; tolerance=1E-12, maxiter=1000)
    # setup the algorithm
#     x_old = iv
    normdiff = Inf
    iter = 1
    if f(x_0) == 0
        return (value = x_0, normdiff=0, iter=0)
    end
    while normdiff > tolerance && iter <= maxiter
        x_new = x_0 - f(x_0)/f_prime(x_0) # use the passed in map
        if norm(x_new - x_0) > tolerance && iter > 100
            println("non-convergence ")
            return #(value = nothing, normdiff=nothing, iter=iter)
        end
        normdiff = norm(x_new - x_0)
        x_0 = x_new
        iter = iter + 1
    end
    return (value = x_0, normdiff=normdiff, iter=iter) # A named tuple
end

D(f) = x -> ForwardDiff.derivative(f, x)
function fixedpointmap(f, x_0; tolerance=1E-12, maxiter=1000)
    f_prime = D(f)
    if nothing == fixedpointmap(f, f_prime, x_0)
        return
    end
    x_0, normdiff, iter = fixedpointmap(f, f_prime, x_0, tolerance=tolerance, maxiter=maxiter)
    return (value = x_0, normdiff=normdiff, iter=iter)
end

export fixedpointmap

end # module

ext: jl


text:
using Documenter, BitsFields

makedocs(
    modules = [BitsFields],
    sitename = "BitsFields",
    pages  = Any[
        "Overview"                 => "index.md",
        "The BitField"             => "thebitfield.md",
        "The BitField, Named"      => "thebitfieldnamed.md",
        "MultiFields"              => "multifields.md",
        "MultiFields, Named"       => "multifieldsnamed.md",
        "IEEE Standard 754-2008"   => "ieeestandard754-2008.md",
        "Worked Use"               => "workeduse.md",
        "References"               => "references.md",
        "Index"                    => "documentindex.md"
        ]
    )

deploydocs(
    repo = "github.com/JeffreySarnoff/BitsFields.jl.git",
    target = "build"
)


text:
const CViewer = Ptr{Cvoid}
struct PetscViewer
    ptr::Ref{CViewer}

    PetscViewer() = new(Ref{CViewer}())
end

# allows us to pass PetscViewer objects directly into CViewer ccall signatures
Base.cconvert(::Type{CViewer}, viewer::PetscViewer) = viewer.ptr[]

"""
    PetscViewerASCIIGetStdout(comm::MPI.Comm, viewer::PetscViewer)

Wrapper for `PetscViewerASCIIGetStdout`
"""
function PetscViewerASCIIGetStdout(comm::MPI.Comm, viewer::PetscViewer)
    error = ccall((:PetscViewerASCIIGetStdout, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr)
    @assert iszero(error)
end

function PetscViewerASCIIGetStdout(comm::MPI.Comm)
    viewer = PetscViewer()
    error = ccall((:PetscViewerASCIIGetStdout, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr)
    @assert iszero(error)
    return viewer
end

PetscViewerStdWorld() = PetscViewerASCIIGetStdout(MPI.COMM_WORLD)

"""
    PetscViewerCreate(comm::MPI.Comm, viewer::PetscViewer)

Wrapper for `PetscViewerCreate`
"""
function PetscViewerCreate(comm::MPI.Comm, viewer::PetscViewer)
    error = ccall((:PetscViewerCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr)
    @assert iszero(error)
end

"""
    PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat)

Wrapper for `PetscViewerPushFormat`
"""
function PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat)
    error = ccall((:PetscViewerPushFormat, libpetsc), PetscErrorCode, (CViewer, PetscViewerFormat), viewer, format)
    @assert iszero(error)
end

"""
    PetscViewerPopFormat(viewer::PetscViewer)

Wrapper for `PetscViewerPopFormat`
"""
function PetscViewerPopFormat(viewer::PetscViewer)
    error = ccall((:PetscViewerPopFormat, libpetsc), PetscErrorCode, (CViewer,), viewer)
    @assert iszero(error)
end

"""
    PetscViewerASCIIOpen(comm::MPI.Comm, filanem, viewer::PetscViewer)

Wrapper for `PetscViewerASCIIOpen`
"""
function PetscViewerASCIIOpen(comm::MPI.Comm, filename, viewer::PetscViewer)
    error = ccall((:PetscViewerASCIIOpen, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Cstring, Ptr{CViewer}), comm, filename, viewer.ptr)
    @assert iszero(error)
end

"""
    PetscViewerDestroy(viewer::PetscViewer)

Wrapper for `PetscViewerDestroy`

Warning : from what I understand, all viewers must not be destroyed explicitely using `PetscViewerDestroy`.

"""
function PetscViewerDestroy(viewer::PetscViewer)
    error = ccall((:PetscViewerDestroy, libpetsc), PetscErrorCode, (Ptr{CViewer}, ), viewer.ptr)
    @assert iszero(error)
end


text:
using WrapperSimplification
using Test

@testset "WrapperSimplification.jl" begin
    # Write your own tests here.
end


text:
function (f::Func)(vars::Variable...)
    f.inputs = [vars...]
    xs = Array{Real}(undef, length(vars))
    xs[1] = vars[1].data
    min_generation = vars[1].generation
    for i in 2:length(vars)
        xs[i] = vars[i].data
        if vars[i].generation < min_generation
            min_generation = vars[i].generation
        end
    end
    ys = forward(f, xs...)
    f.generation = min_generation
    if length(ys) == 1
        outputs = [Variable(ys[1], f, nothing, f.generation - 1, ""),]
        f.outputs = outputs
        return outputs[1]
    else
        outputs = [Variable(y, f, nothing, f.generation - 1, "") for y in ys]
        f.outputs = outputs
        return  outputs
    end
end

function backward!(var::Variable)
    (isnothing(var.grad)) && (var.grad = variable(ones_like(var.data)))
    funcs = PriorityQueue{Func,Int}()
    seen_set = Set{Func}()
    enqueue!(funcs, var.creator, 1)
    push!(seen_set, var.creator)
    while !(isempty(funcs))
        f = dequeue!(funcs)
        gxs = backward(f)
        for (x, gx) in zip(f.inputs, gxs)
            if x.grad === nothing
                x.grad = gx
            else
                x.grad = x.grad + gx
            end
            if (x.creator !== nothing) && (!(x.creator in seen_set))
                push!(seen_set, x.creator)
                enqueue!(funcs, x.creator, x.creator.generation)
            end
        end
    end
end

function variable(data, creator=nothing, grad=nothing, generation=0; name="")
    return Variable(data, creator, grad, generation, name)
end

function get_gys(outputs)
    return [output.grad for output in outputs]
end

function get_data(var::Variable)
    return var.data
end

function get_creator(var::Variable)
    return var.creator
end

function get_grad(var::Variable)
    return var.grad
end

function get_generation(var::Variable)
    return var.generation
end

function get_outputs(f::Func)
    return f.outputs
end

function ones_like(x)
    shape = size(x)
    return isempty(shape) ? 1 : ones(shape)
end

function as_tuple(x)
    typeof(x) <: Tuple ? x : tuple(x)
end

function as_variable(x)
    typeof(x) == Variable ? x : variable(x)
end

function cleargrad!(var::Variable)
    var.grad = nothing
end

function set_name(var::Variable, name)
    var.name = name
end


text:
using Documenter, sem

makedocs(;
    modules=[sem],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/Maximilian-Stefan-Ernst/sem.jl/blob/{commit}{path}#L{line}",
    sitename="sem.jl",
    authors="Maximilian Ernst, Aaron Peikert",
    assets=String[],
)

deploydocs(;
    repo="github.com/Maximilian-Stefan-Ernst/sem.jl",
)


text:
mutable struct FTestResult{N}
    nobs::Int
    ssr::NTuple{N, Float64}
    dof::NTuple{N, Int}
    r2::NTuple{N, Float64}
    fstat::NTuple{N, Float64}
    pval::NTuple{N, Float64}
end

"""A helper function to determine if mod1 is nested in mod2"""
function issubmodel(mod1::LinPredModel, mod2::LinPredModel; atol::Real=0.0)
    mod1.rr.y != mod2.rr.y && return false # Response variables must be equal

    # Test that models are nested
    pred1 = mod1.pp.X
    npreds1 = size(pred1, 2)
    pred2 = mod2.pp.X
    npreds2 = size(pred2, 2)
    # If model 1 has more predictors, it can't possibly be a submodel
    npreds1 > npreds2 && return false
    # Test min norm pred2*B - pred1 ≈ 0
    rtol = Base.rtoldefault(typeof(pred1[1,1]))
    nresp = size(pred2, 1)
    return norm(view(qr(pred2).Q'pred1, npreds2 + 1:nresp, :)) <= max(atol, rtol*norm(pred1))
end

_diffn(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i]-t[i+1], N-1)

_diff(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i+1]-t[i], N-1)

"""
    ftest(mod::LinearModel...; atol::Real=0.0)

For each sequential pair of linear models in `mod...`, perform an F-test to determine if
the one model fits significantly better than the other. Models must have been fitted
on the same data, and be nested either in forward or backward direction.

A table is returned containing consumed degrees of freedom (DOF),
difference in DOF from the preceding model, sum of squared residuals (SSR), difference in
SSR from the preceding model, R², difference in R² from the preceding model, and F-statistic
and p-value for the comparison between the two models.

!!! note
    This function can be used to perform an ANOVA by testing the relative fit of two models
    to the data

Optional keyword argument `atol` controls the numerical tolerance when testing whether
the models are nested.

# Examples

Suppose we want to compare the effects of two or more treatments on some result. Because
this is an ANOVA, our null hypothesis is that `Result ~ 1` fits the data as well as
`Result ~ 1 + Treatment`.

```jldoctest
julia> dat = DataFrame(Result=[1.1, 1.2, 1, 2.2, 1.9, 2, 0.9, 1, 1, 2.2, 2, 2],
                       Treatment=[1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2],
                       Other=categorical([1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 2, 1]));

julia> nullmodel = lm(@formula(Result ~ 1), dat);

julia> model = lm(@formula(Result ~ 1 + Treatment), dat);

julia> bigmodel = lm(@formula(Result ~ 1 + Treatment + Other), dat);

julia> ftest(nullmodel.model, model.model)
F-test: 2 models fitted on 12 observations
─────────────────────────────────────────────────────────────────
     DOF  ΔDOF     SSR     ΔSSR       R²     ΔR²        F*  p(>F)
─────────────────────────────────────────────────────────────────
[1]    2        3.2292           -0.0000
[2]    3     1  0.1283  -3.1008   0.9603  0.9603  241.6234  <1e-7
─────────────────────────────────────────────────────────────────

julia> ftest(nullmodel.model, model.model, bigmodel.model)
F-test: 3 models fitted on 12 observations
──────────────────────────────────────────────────────────────────
     DOF  ΔDOF     SSR     ΔSSR       R²     ΔR²        F*   p(>F)
──────────────────────────────────────────────────────────────────
[1]    2        3.2292           -0.0000
[2]    3     1  0.1283  -3.1008   0.9603  0.9603  241.6234   <1e-7
[3]    5     2  0.1017  -0.0266   0.9685  0.0082    1.0456  0.3950
──────────────────────────────────────────────────────────────────
```
"""
function ftest(mods::LinearModel...; atol::Real=0.0)
    if length(mods) < 2
        Base.depwarn("Passing less than two models to ftest is deprecated", :ftest)
    end
    if !all(==(nobs(mods[1])), nobs.(mods))
        throw(ArgumentError("F test is only valid for models fitted on the same data, " *
                            "but number of observations differ"))
    end
    forward = length(mods) == 1 || dof(mods[1]) <= dof(mods[2])
    if forward
        for i in 2:length(mods)
            if dof(mods[i-1]) >= dof(mods[i]) || !issubmodel(mods[i-1], mods[i], atol=atol)
                throw(ArgumentError("F test is only valid for nested models"))
            end
        end
    else
        for i in 2:length(mods)
            if dof(mods[i]) >= dof(mods[i-1]) || !issubmodel(mods[i], mods[i-1], atol=atol)
                throw(ArgumentError("F test is only valid for nested models"))
            end
        end
    end

    SSR = deviance.(mods)

    df = dof.(mods)
    Δdf = _diff(df)
    dfr = Int.(dof_residual.(mods))
    MSR1 = _diffn(SSR) ./ Δdf
    MSR2 = (SSR ./ dfr)
    if forward
        MSR2 = MSR2[2:end]
        dfr_big = dfr[2:end]
    else
        MSR2 = MSR2[1:end-1]
        dfr_big = dfr[1:end-1]
    end

    fstat = (NaN, (MSR1 ./ MSR2)...)
    pval = (NaN, ccdf.(FDist.(abs.(Δdf), dfr_big), abs.(fstat[2:end]))...)
    return FTestResult(Int(nobs(mods[1])), SSR, df, r2.(mods), fstat, pval)
end

function show(io::IO, ftr::FTestResult{N}) where N
    Δdof = _diff(ftr.dof)
    Δssr = _diff(ftr.ssr)
    ΔR² = _diff(ftr.r2)

    nc = 9
    nr = N
    outrows = Matrix{String}(undef, nr+1, nc)

    outrows[1, :] = ["", "DOF", "ΔDOF", "SSR", "ΔSSR",
                     "R²", "ΔR²", "F*", "p(>F)"]

    outrows[2, :] = ["[1]", @sprintf("%.0d", ftr.dof[1]), " ",
                     @sprintf("%.4f", ftr.ssr[1]), " ",
                     @sprintf("%.4f", ftr.r2[1]), " ", " ", " "]

    for i in 2:nr
        outrows[i+1, :] = ["[$i]",
                           @sprintf("%.0d", ftr.dof[i]), @sprintf("%.0d", Δdof[i-1]),
                           @sprintf("%.4f", ftr.ssr[i]), @sprintf("%.4f", Δssr[i-1]),
                           @sprintf("%.4f", ftr.r2[i]), @sprintf("%.4f", ΔR²[i-1]),
                           @sprintf("%.4f", ftr.fstat[i]), string(PValue(ftr.pval[i])) ]
    end
    colwidths = length.(outrows)
    max_colwidths = [maximum(view(colwidths, :, i)) for i in 1:nc]
    totwidth = sum(max_colwidths) + 2*8

    println(io, "F-test: $N models fitted on $(ftr.nobs) observations")
    println(io, '─'^totwidth)

    for r in 1:nr+1
        for c in 1:nc
            cur_cell = outrows[r, c]
            cur_cell_len = length(cur_cell)

            padding = " "^(max_colwidths[c]-cur_cell_len)
            if c > 1
                padding = "  "*padding
            end

            print(io, padding)
            print(io, cur_cell)
        end
        print(io, "\n")
        r == 1 && println(io, '─'^totwidth)
    end
    print(io, '─'^totwidth)
end


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

@doc raw"""EndpointPort is a tuple that describes a single port.

    IoK8sApiCoreV1EndpointPort(;
        name=nothing,
        port=nothing,
        protocol=nothing,
    )

    - name::String : The name of this port.  This must match the 'name' field in the corresponding ServicePort. Must be a DNS_LABEL. Optional only if one port is defined.
    - port::Int32 : The port number of the endpoint.
    - protocol::String : The IP protocol for this port. Must be UDP, TCP, or SCTP. Default is TCP.
"""
mutable struct IoK8sApiCoreV1EndpointPort <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiCoreV1EndpointPort(;name=nothing, port=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EndpointPort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1EndpointPort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiCoreV1EndpointPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiCoreV1EndpointPort

const _property_map_IoK8sApiCoreV1EndpointPort = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiCoreV1EndpointPort = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("port")=>"Int32", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1EndpointPort }) = collect(keys(_property_map_IoK8sApiCoreV1EndpointPort))
Swagger.property_type(::Type{ IoK8sApiCoreV1EndpointPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EndpointPort[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EndpointPort }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EndpointPort[property_name]

function check_required(o::IoK8sApiCoreV1EndpointPort)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

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


text:
module ConcurrentCollectionsBenchmarks

using BenchmarkTools: Benchmark, BenchmarkGroup

include("utils.jl")
include("bench_dict_histogram.jl")
include("bench_dict_haskey.jl")
include("bench_dict_get_existing.jl")
include("bench_dict_migration.jl")
include("bench_queue_pushpop.jl")
include("bench_queue_hot_potato.jl")

function setup()
    suite = BenchmarkGroup()
    suite["DictHistogram"] = BenchDictHistogram.setup()
    suite["DictGetExisting"] = BenchDictGetExisting.setup()
    suite["DictMigration"] = BenchDictMigration.setup()
    suite["QueuePushPop"] = BenchQueuePushPop.setup()
    suite["HotPotato"] = BenchQueueHotPotato.setup()
    return suite
end

function set_smoke_params!(bench)
    bench.params.seconds = 0.001
    bench.params.evals = 1
    bench.params.samples = 1
    bench.params.gctrial = false
    bench.params.gcsample = false
    return bench
end

foreach_benchmark(f!, bench::Benchmark) = f!(bench)
function foreach_benchmark(f!, group::BenchmarkGroup)
    for x in values(group)
        foreach_benchmark(f!, x)
    end
end

function setup_smoke()
    suite = setup()
    foreach_benchmark(set_smoke_params!, suite)
    return suite
end

function clear()
    BenchDictHistogram.clear()
    BenchDictGetExisting.clear()
    BenchDictMigration.clear()
    BenchQueuePushPop.clear()
    BenchQueueHotPotato.clear()
end

end # module


text:
"""
"""
abstract type SingleFieldFESpace <: FESpace end

"""
"""
function get_dirichlet_dof_ids(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
num_dirichlet_dofs(f::SingleFieldFESpace) = length(get_dirichlet_dof_ids(f))

"""
"""
function zero_dirichlet_values(f::SingleFieldFESpace)
  V = get_vector_type(f)
  allocate_vector(V,num_dirichlet_dofs(f))
end

"""
"""
function num_dirichlet_tags(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function get_dirichlet_dof_tag(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function scatter_free_and_dirichlet_values(f::SingleFieldFESpace,free_values,dirichlet_values)
  @abstractmethod
end

"""
"""
function gather_free_and_dirichlet_values!(free_values, dirichlet_values,fs::SingleFieldFESpace,cell_vals)
  @abstractmethod
end

"""
"""
function get_cell_dof_ids(V::SingleFieldFESpace,trian::Triangulation)
  trian_V = get_triangulation(V)
  if have_compatible_domains(trian_V,trian)
    get_cell_dof_ids(V)
  elseif have_compatible_domains(trian_V,get_background_triangulation(trian))
    get_cell_dof_ids(V,get_cell_to_bgcell(trian))
  elseif have_compatible_domains(
    get_background_triangulation(trian_V),get_background_triangulation(trian))
    cell_to_Vcell = get_cell_to_bgcell(trian,trian_V)
    get_cell_dof_ids(V,cell_to_Vcell)
  elseif have_compatible_domains(
    trian_V,get_background_triangulation(get_background_triangulation(trian)))
    bg_trian = get_background_triangulation(trian)
    bg_cell_dof_ids = get_cell_dof_ids(V,bg_trian)
    cell_to_bgcell = get_cell_to_bgcell(trian)
    get_cell_dof_ids(bg_cell_dof_ids,cell_to_bgcell)
  else
    @unreachable
  end
end

function get_cell_isconstrained(V::SingleFieldFESpace,trian::Triangulation)
  trian_V = get_triangulation(V)
  if have_compatible_domains(trian_V,trian)
    get_cell_isconstrained(V)
  elseif have_compatible_domains(trian_V,get_background_triangulation(trian))
    get_cell_isconstrained(V,get_cell_to_bgcell(trian))
  elseif have_compatible_domains(
    get_background_triangulation(trian_V),get_background_triangulation(trian))
    cell_to_Vcell = get_cell_to_bgcell(trian,trian_V)
    get_cell_isconstrained(V,cell_to_Vcell)
  elseif have_compatible_domains(
    trian_V,get_background_triangulation(get_background_triangulation(trian)))
    bg_trian = get_background_triangulation(trian)
    bg_cell_isconstrained = get_cell_isconstrained(V,bg_trian)
    cell_to_bgcell = get_cell_to_bgcell(trian)
    get_cell_isconstrained(bg_cell_isconstrained,cell_to_bgcell)
  else
    @unreachable
  end
end

function get_cell_constraints(V::SingleFieldFESpace,trian::Triangulation)
  trian_V = get_triangulation(V)
  if have_compatible_domains(trian_V,trian)
    get_cell_constraints(V)
  elseif have_compatible_domains(trian_V,get_background_triangulation(trian))
    get_cell_constraints(V,get_cell_to_bgcell(trian))
  elseif have_compatible_domains(
    get_background_triangulation(trian_V),get_background_triangulation(trian))
    cell_to_Vcell = get_cell_to_bgcell(trian,trian_V)
    get_cell_constraints(V,cell_to_Vcell)
  elseif have_compatible_domains(
    trian_V,get_background_triangulation(get_background_triangulation(trian)))
    bg_trian = get_background_triangulation(trian)
    bg_cell_constraints = get_cell_constraints(V,bg_trian)
    cell_to_bgcell = get_cell_to_bgcell(trian)
    get_cell_constraints(bg_cell_constraints,cell_to_bgcell)
  else
    @unreachable
  end
end

function get_cell_is_dirichlet(V::SingleFieldFESpace,trian::Triangulation)
  trian_V = get_triangulation(V)
  if have_compatible_domains(trian_V,trian)
    get_cell_is_dirichlet(V)
  elseif have_compatible_domains(trian_V,get_background_triangulation(trian))
    get_cell_is_dirichlet(V,get_cell_to_bgcell(trian))
  elseif have_compatible_domains(
    get_background_triangulation(trian_V),get_background_triangulation(trian))
    cell_to_Vcell = get_cell_to_bgcell(trian,trian_V)
    get_cell_is_dirichlet(V,cell_to_Vcell)
  elseif have_compatible_domains(
    trian_V,get_background_triangulation(get_background_triangulation(trian)))
    bg_trian = get_background_triangulation(trian)
    bg_cell_is_dirichlet = get_cell_is_dirichlet(V,bg_trian)
    cell_to_bgcell = get_cell_to_bgcell(trian)
    get_cell_is_dirichlet(bg_cell_is_dirichlet,cell_to_bgcell)
  else
    @unreachable
  end
end

function get_cell_is_dirichlet(V::SingleFieldFESpace)
  trian = get_triangulation(V)
  Fill(true,num_cells(trian))
end

"""
"""
function test_single_field_fe_space(f::SingleFieldFESpace,pred=(==))
  fe_basis = get_cell_shapefuns(f)
  @test isa(fe_basis,CellField)
  test_fe_space(f)
  cell_dofs = get_cell_dof_ids(f)
  dirichlet_values = zero_dirichlet_values(f)
  @test length(dirichlet_values) == num_dirichlet_dofs(f)
  free_values = zero_free_values(f)
  cell_vals = scatter_free_and_dirichlet_values(f,free_values,dirichlet_values)
  fv, dv = gather_free_and_dirichlet_values(f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  fv, dv = gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  fe_function = FEFunction(f,free_values,dirichlet_values)
  @test isa(fe_function, SingleFieldFEFunction)
  test_fe_function(fe_function)
  ddof_to_tag = get_dirichlet_dof_tag(f)
  @test length(ddof_to_tag) == num_dirichlet_dofs(f)
  if length(get_dirichlet_dof_tag(f)) != 0
    @test maximum(get_dirichlet_dof_tag(f)) <= num_dirichlet_tags(f)
  end
  cell_dof_basis = get_cell_dof_basis(f)
  @test isa(cell_dof_basis,CellDof)
end

function test_single_field_fe_space(f,matvecdata,matdata,vecdata,pred=(==))
  test_single_field_fe_space(f,pred)
  test_fe_space(f,matvecdata,matdata,vecdata)
end

# Some default API

"""
"""
function get_dirichlet_values(f::SingleFieldFESpace)
  zero_dirichlet_values(f)
end

"""
"""
function gather_free_and_dirichlet_values(f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  dirichlet_values = zero_dirichlet_values(f)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
end

"""
"""
function gather_dirichlet_values(f::SingleFieldFESpace,cell_vals)
  dirichlet_values = zero_dirichlet_values(f)
  gather_dirichlet_values!(dirichlet_values,f,cell_vals)
  dirichlet_values
end

"""
"""
function gather_dirichlet_values!(dirichlet_values,f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
  dirichlet_values
end

"""
"""
function gather_free_values(f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  gather_free_values!(free_values,f,cell_vals)
  free_values
end

"""
"""
function gather_free_values!(free_values,f::SingleFieldFESpace,cell_vals)
    dirichlet_values = zero_dirichlet_values(f)
    gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
    free_values
end

function CellField(fs::SingleFieldFESpace,cell_vals)
  v = get_cell_shapefuns(fs)
  cell_basis = get_data(v)
  cell_field = lazy_map(linear_combination,cell_vals,cell_basis)
  GenericCellField(cell_field,get_triangulation(v),DomainStyle(v))
end

struct SingleFieldFEFunction{T<:CellField} <: FEFunction
  cell_field::T
  cell_dof_values::AbstractArray{<:AbstractVector{<:Number}}
  free_values::AbstractVector{<:Number}
  dirichlet_values::AbstractVector{<:Number}
  fe_space::SingleFieldFESpace
end

get_data(f::SingleFieldFEFunction) = get_data(f.cell_field)
get_triangulation(f::SingleFieldFEFunction) = get_triangulation(f.cell_field)
DomainStyle(::Type{SingleFieldFEFunction{T}}) where T = DomainStyle(T)

get_free_dof_values(f::SingleFieldFEFunction) = f.free_values
get_cell_dof_values(f::SingleFieldFEFunction) = f.cell_dof_values
get_fe_space(f::SingleFieldFEFunction) = f.fe_space

function get_cell_dof_values(f::SingleFieldFEFunction,trian::Triangulation)
  trian_f = get_triangulation(f)
  if have_compatible_domains(trian_f,trian)
    get_cell_dof_values(f)
  elseif have_compatible_domains(trian_f,get_background_triangulation(trian))
    get_cell_dof_values(f,get_cell_to_bgcell(trian))
  elseif have_compatible_domains(
    get_background_triangulation(trian_f),get_background_triangulation(trian))
    cell_to_Vcell = get_cell_to_bgcell(trian,trian_f)
    get_cell_dof_values(f,cell_to_Vcell)
  elseif have_compatible_domains(
    trian_f,get_background_triangulation(get_background_triangulation(trian)))
    bg_trian = get_background_triangulation(trian)
    bg_cell_dof_values = get_cell_dof_values(f,bg_trian)
    cell_to_bgcell = get_cell_to_bgcell(trian)
    get_cell_dof_values(bg_cell_dof_values,cell_to_bgcell)
  else
    @unreachable
  end
end

"""
    FEFunction(
      fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector)

The resulting FEFunction will be in the space if and only if `dirichlet_values`
are the ones provided by `get_dirichlet_values(fs)`
"""
function FEFunction(
  fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector)
  cell_vals = scatter_free_and_dirichlet_values(fs,free_values,dirichlet_values)
  cell_field = CellField(fs,cell_vals)
  SingleFieldFEFunction(cell_field,cell_vals,free_values,dirichlet_values,fs)
end

function FEFunction(fe::SingleFieldFESpace, free_values)
  diri_values = get_dirichlet_values(fe)
  FEFunction(fe,free_values,diri_values)
end

"""
The resulting FE function is in the space (in particular it fulfills Dirichlet BCs
even in the case that the given cell field does not fulfill them)
"""
function interpolate(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  interpolate!(object, free_values,fs)
end

"""
"""
function interpolate!(object, free_values,fs::SingleFieldFESpace)
    cell_vals = _cell_vals(fs,object)
    gather_free_values!(free_values,fs,cell_vals)
    FEFunction(fs,free_values)
end

function _cell_vals(fs::SingleFieldFESpace,object)
  s = get_cell_dof_basis(fs)
  trian = get_triangulation(s)
  f = CellField(object,trian,DomainStyle(s))
  cell_vals = s(f)
end

"""
like interpolate, but also compute new degrees of freedom for the dirichlet component.
The resulting FEFunction does not necessary belongs to the underlying space
"""
function interpolate_everywhere(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  dirichlet_values = zero_dirichlet_values(fs)
  interpolate_everywhere!(object, free_values,dirichlet_values,fs)
end

"""
"""
function interpolate_everywhere!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
  cell_vals = _cell_vals(fs,object)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,fs,cell_vals)
  FEFunction(fs,free_values,dirichlet_values)
end

"""
"""
function interpolate_dirichlet(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  dirichlet_values = zero_dirichlet_values(fs)
  interpolate_dirichlet!(object, free_values,dirichlet_values,fs)
end

"""
"""
function interpolate_dirichlet!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
  cell_vals = _cell_vals(fs,object)
  gather_dirichlet_values!(dirichlet_values,fs,cell_vals)
  fill!(free_values,zero(eltype(free_values)))
  FEFunction(fs,free_values, dirichlet_values)
end

"""
"""
function compute_dirichlet_values_for_tags(f::SingleFieldFESpace,tag_to_object)
  dirichlet_values = zero_dirichlet_values(f)
  dirichlet_values_scratch = zero_dirichlet_values(f)
  compute_dirichlet_values_for_tags!(dirichlet_values,dirichlet_values_scratch,f,tag_to_object)
end

function compute_dirichlet_values_for_tags!(
  dirichlet_values,
  dirichlet_values_scratch,
  f::SingleFieldFESpace,tag_to_object)

  dirichlet_dof_to_tag = get_dirichlet_dof_tag(f)
  _tag_to_object = _convert_to_collectable(tag_to_object,num_dirichlet_tags(f))
  for (tag, object) in enumerate(_tag_to_object)
    cell_vals = _cell_vals(f,object)
    fill!(dirichlet_values_scratch,zero(eltype(dirichlet_values_scratch)))
    gather_dirichlet_values!(dirichlet_values_scratch,f,cell_vals)
    _fill_dirichlet_values_for_tag!(dirichlet_values,dirichlet_values_scratch,tag,dirichlet_dof_to_tag)
  end
  dirichlet_values
end

function _fill_dirichlet_values_for_tag!(dirichlet_values,dv,tag,dirichlet_dof_to_tag)
  for (dof, v) in enumerate(dv)
    if dirichlet_dof_to_tag[dof] == tag
      dirichlet_values[dof] = v
    end
  end
end

function _convert_to_collectable(object,ntags)
  @assert ntags == length(object) "Incorrect number of dirichlet tags provided"
  object
end

function _convert_to_collectable(object::Function,ntags)
  _convert_to_collectable(Fill(object,ntags),ntags)
end

function _convert_to_collectable(object::Number,ntags)
  _convert_to_collectable(Fill(object,ntags),ntags)
end


text:
import Cartesian
using Base.Test

# @nloops with range expression
for N = 1:4
    @eval begin
        function loopsum{T}(A::StridedArray{T,$N})
            s = zero(eltype(A))
            @inbounds Cartesian.@nloops $N i d->1:size(A,d) begin
                s += Cartesian.@nref $N A i
            end
            s
        end
    end
end

# @nloops with size determined by array
for N = 1:4
    @eval begin
        function loopsum2{T}(A::StridedArray{T,$N})
            s = zero(eltype(A))
            @inbounds Cartesian.@nloops $N i A begin
                s += Cartesian.@nref $N A i
            end
            s
        end
    end
end

A = reshape(1:1000*1001, 1000, 1001)
S = sub(A, 2:999, 2:1000)

As = sum(A)
@assert loopsum(A) == As
@assert loopsum2(A) == As
Ss = sum(S)
@assert loopsum(S) == Ss
@assert loopsum2(S) == Ss

tbase = @elapsed (for k = 1:100; sum(A); end)
t = @elapsed (for k = 1:100; loopsum(A); end)
# @assert t < 2tbase
t = @elapsed (for k = 1:100; loopsum2(A); end)
# @assert t < 2tbase

t = @elapsed (for k = 1:100; loopsum(S); end)
# @assert t < 3tbase
t = @elapsed (for k = 1:100; loopsum2(S); end)
# @assert t < 3tbase

A = reshape(1:1000*1001*20, 1000, 1001, 20)
S = sub(A, 2:999, 2:1000, 2:19)

As = sum(A)
@assert loopsum(A) == As
@assert loopsum2(A) == As
Ss = sum(S)
@assert loopsum(S) == Ss
@assert loopsum2(S) == Ss

tbase = @elapsed sum(A)
t = @elapsed loopsum(A)
# @assert t < 2tbase
t = @elapsed loopsum2(A)
# @assert t < 2tbase
t = @elapsed loopsum(S)
# @assert t < 3tbase
t = @elapsed loopsum2(S)
# @assert t < 3tbase

# @nloops with pre-expression
for N = 1:4
    @eval begin
        function maxoverdims{T}(A::AbstractArray{T,$N}, region)
            szout = [size(A,d) for d = 1:$N]
            szout[[region...]] = 1
            B = fill(typemin(T), szout...)::Array{T,$N}
            Cartesian.@nextract $N szout szout
            Cartesian.@nloops $N i A d->(j_d = szout_d==1 ? 1 : i_d) begin
                (Cartesian.@nref $N B j) = max((Cartesian.@nref $N B j), (Cartesian.@nref $N A i))
            end
            B
        end
    end
end

A = reshape(1:10,5,2)
A1 = maxoverdims(A, 1)
@assert A1 == [5,10]'
A2 = maxoverdims(A, 2)
@assert A2 == reshape(6:10,5,1)
@assert maxoverdims(A, (1,2)) == reshape([10], 1, 1)

# Curly-brace syntax: sum over the upper-triangle
A = reshape(1:16, 4, 4)
s = 0
Cartesian.@nloops 2 i d->d==2?(1:size(A,d)):(1:i_{d+1}) begin
    s += Cartesian.@nref 2 A i
end
@assert s == sum(triu(A))

# @nref, @nrefshift, @nextract, and @nlookup
A = reshape(1:15, 3, 5)
i_1 = 2
i_2 = 3
@assert (Cartesian.@nref 2 A i) == A[i_1, i_2]
j_1 = -1
j_2 = 2
@assert (Cartesian.@nrefshift 2 A i j) == A[i_1+j_1, i_2+j_2]
@assert (Cartesian.@nrefshift 2 A i d->(d==2)?1:0) == A[i_1, i_2+1]
I = ([i_1], [i_2])
Cartesian.@nextract 2 k I
@assert k_1 == [i_1]
@assert k_2 == [i_2]
j_1 = 1
j_2 = 1
@assert (Cartesian.@nlookup 2 A k j) == A[i_1, i_2]

# @nlinear
A = reshape(1:120, 3, 4, 10)
i_1 = 2
i_2 = 2
i_3 = 7
p, index = Cartesian.@nlinear 3 A i
@assert index == A[i_1, i_2, i_3]

# The i_d notation
i_1 = 2
i_2 = -1
pairs = {}
Cartesian.@nloops 2 j d->(1-min(0,i_d):4-max(0,i_d)) begin
    push!(pairs, (j_1,j_2))
end
@assert pairs == {(1,2),(2,2),(1,3),(2,3),(1,4),(2,4)}

# @nall
pairs = {}
Cartesian.@nloops 2 j d->1:4 begin
    if Cartesian.@nall 2 d->(1 <= j_d+i_d <= 4)
        push!(pairs, (j_1,j_2))
    end
end
@assert pairs == {(1,2),(2,2),(1,3),(2,3),(1,4),(2,4)}

# @nexprs
A = reshape(1:20*7, 20, 7)
indexes = (2:5:20,3:7)
strds = strides(A)
i_1 = 2
i_2 = 3
ind = 1
@assert (Cartesian.@nexprs 2 d->(ind += (indexes[d][i_d]-1)*strds[d])) == A[indexes[1][i_1],indexes[2][i_2]]

# @ngenerate
Cartesian.@ngenerate N T function loopsum3{T,N}(A::StridedArray{T,N})
    s = zero(eltype(A))
    @inbounds Cartesian.@nloops(N, i, A, begin
        s += Cartesian.@nref(N, A, i)
    end)
    s
end
A = reshape(1:8, 2, 2, 2)
@test loopsum3(A) == 36
A = reshape(1:32, 2, 2, 2, 2, 2)
@test loopsum3(A) == 528


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 04be2568-92b7-11eb-1d59-97cd595d320d
#Set up the working environment
begin
	using Pkg
	Pkg.activate("..")
	using Plots, PlutoUI, LaTeXStrings, DifferentialEquations, DataFrames, DataFramesMeta, CSV, StatsPlots, Dates, Flux, Optim, DiffEqFlux, DiffEqSensitivity
end

# ╔═╡ 39a8b108-92d0-11eb-2442-1f47a7e9b445
md"""
# Interactive exploration of the "Our World in Data" Covid-19 Dataset

More information: [https://github.com/owid/covid-19-data/tree/master/public/data](https://github.com/owid/covid-19-data/tree/master/public/data)
"""

# ╔═╡ 53479e94-92b7-11eb-2731-c7e03afeb549
# Download dataset
covid_data = DataFrame(CSV.File(download("https://covid.ourworldindata.org/data/owid-covid-data.csv")));

# ╔═╡ cd146c8a-92bb-11eb-2811-c541ca9cef7b
md"""
| Countries     | y-axis      |
|---------------|-------------|
| $(@bind locations MultiSelect(unique(covid_data.location), default=["Austria"])) | $(@bind ydata MultiSelect(names(covid_data), default=["total_cases"])) |
"""

# ╔═╡ 059abcfc-92b8-11eb-0bcc-6f688494a3cb
# Filter data for Austria
D = covid_data[in.(covid_data.location, Ref(locations)), :];

# ╔═╡ bef90d4e-92c8-11eb-2539-19bf3fc70791
md"""
| Plot range   | value                       |
|--------------|-----------------------------|
| xmin         | $(@bind xmn Slider(1:1:length(unique(D[!,"date"]))-1, default=0)) |
| xmax         | $(@bind xmx Slider(1:1:length(unique(D[!,"date"]))-1, default=length(unique(D[!,"date"]))-1)) |
"""

# ╔═╡ 68b6d4b2-92cd-11eb-1d2e-5de1d1be1d95
begin
	xmin = minimum(D[!,"date"]) + Day(xmn);
	xmax = minimum(D[!,"date"]) + Day(xmx);
end;

# ╔═╡ b5dc1850-92c1-11eb-2934-6788cb4a2d6a
begin
	dataplot = @df D plot(cols(Symbol("date")), cols([Symbol(y) for y in ydata]),
			group = (:location),
			xlims = (xmin, xmax),
			labels = reshape(["$l, $y" for y in ydata for l in locations], 1, :),
			fontfamily = "ComputerModern",
			legend = :topleft)

	plot(dataplot)
end

# ╔═╡ Cell order:
# ╟─39a8b108-92d0-11eb-2442-1f47a7e9b445
# ╠═04be2568-92b7-11eb-1d59-97cd595d320d
# ╠═53479e94-92b7-11eb-2731-c7e03afeb549
# ╠═059abcfc-92b8-11eb-0bcc-6f688494a3cb
# ╟─cd146c8a-92bb-11eb-2811-c541ca9cef7b
# ╟─bef90d4e-92c8-11eb-2539-19bf3fc70791
# ╟─68b6d4b2-92cd-11eb-1d2e-5de1d1be1d95
# ╠═b5dc1850-92c1-11eb-2934-6788cb4a2d6a

max_stars_repo_name: el-uhu/intro-to-modelling

max_issues_repo_name: el-uhu/intro-to-modelling

max_forks_repo_name: el-uhu/intro-to-modelling


text:
using Distributions, Plots; pyplot()

L, K, n  = 500, [450, 400, 250, 100, 50], 30
hyperDists = [Hypergeometric(k,L-k,n) for k in K]
xGrid = 0:1:n
pmfs = [ pdf.(dist, xGrid) for dist in hyperDists ]
labels = "Successes = " .* string.(K)

bar( xGrid, pmfs,
	alpha=0.8, c=[:orange :purple :green :red :blue ],
	label=hcat(labels...), ylims=(0,0.25),
	xlabel="x", ylabel="Probability", legend=:top)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Bukdu

# https://github.com/zalandoresearch/fashion-mnist

struct FashionController <: ApplicationController
    conn::Conn
end

function index(c::FashionController)
    FashionClothes  = first.(split("""
👕
👖
🏂
👗
🧥
👡
👔
👟
👜
👢  """, '\n'))
    x = rand(1:10)
    render(HTML, string(x, ' ', FashionClothes[x]))
end

if PROGRAM_FILE == basename(@__FILE__)

routes() do
    get("/", FashionController, index)
end

Bukdu.start(8080)

# Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module FiniteDifference

export upwind,downwind,central

"""
    this function's inputs are w and Δz

        w has 1 row
        the 1st row is h (enthalpy)

        Δz is the space interval

    this function's output is out

        out has the same structure as w

    this function uses upwind finite difference scheme to get the 1st order spacial derivative

    The current function is only good for information tavelling to the right
"""

function upwind(w::Array,Δz::Float64)
    out=zero(deepcopy(w))
    NX=length(w)
    @inbounds for x in 2:NX-1
        out[x] =  (w[x] - w[x-1])/Δz
    end

    out[1]=out[2]
    out[NX]=out[NX-1]

    return out
end

"""
    this function's inputs are w and Δz

        w has 1 row
        the 1st row is h (enthalpy)

        Δz is the space interval

    this function's output is out

        out has the same structure as w

    this function uses downwind finite difference scheme to get the 1st order spacial derivative

    The current function is only good for information tavelling to the left
"""
function downwind(w::Array,Δz::Float64)
    out=zero(deepcopy(w))
    NX=length(w)
    @inbounds for x in 2:NX-1
        out[x] =  (w[x+1] - w[x])/Δz
    end

    out[1]=out[2]
    out[NX]=out[NX-1]

    return out
end

"""
this function's inputs are w and Δz

    w has 1 row
    the 1st row is h (enthalpy)

    Δz is the space interval

this function's output is out

    out has the same structure as w

this function uses central finite difference scheme to get the 2st order spacial derivative
"""

function central(w::Array,Δz::Float64)
    out=zero(deepcopy(w))
    NX=length(w)
    @inbounds for x in 2:NX-1
        out[x] =  (w[x+1] - w[x-1])/(2Δz)
    end

    out[1]=out[2]
    out[NX]=out[NX-1]

    return out
end

end


text:
using SeaPearl

struct OutputDataQueens
    nb_sols ::Int
    indices ::Matrix{Int}
end

"""
    struct MostCenteredVariableSelection{TakeObjective}

VariableSelection heuristic that selects the legal (ie. among the not bounded ones) most centered Queen.
"""
struct MostCenteredVariableSelection{TakeObjective} <: SeaPearl.AbstractVariableSelection{TakeObjective} end
###constructor###
MostCenteredVariableSelection(;take_objective = true) = MostCenteredVariableSelection{take_objective}()

function (::MostCenteredVariableSelection{false})(cpmodel::SeaPearl.CPModel)::SeaPearl.AbstractIntVar
    selectedVar = nothing
    n = length(cpmodel.variables)
    sorted_dict = sort(collect(SeaPearl.branchable_variables(cpmodel)),by = x -> abs(n/2-parse(Int, match(r"[0-9]*$", x[1]).match)),rev=true)
    while !isempty(sorted_dict)
        selectedVar = pop!(sorted_dict)[2]
        if !(selectedVar == cpmodel.objective) && !SeaPearl.isbound(selectedVar)
            break
        end
    end

    if SeaPearl.isnothing(selectedVar) && !SeaPearl.isbound(cpmodel.objective)
        return cpmodel.objective
    end
    return selectedVar
end

function (::MostCenteredVariableSelection{true})(cpmodel::SeaPearl.CPModel)::SeaPearl.AbstractIntVar
    selectedVar = nothing
    n = length(cpmodel.variables)
    sorted_dict = sort(collect(SeaPearl.branchable_variables(cpmodel)),by = x -> abs(n/2-parse(Int, match(r"[0-9]*$", x[1]).match)),rev=true)
    while true
        selectedVar= pop!(sorted_dict)[2]
        if !SeaPearl.isbound(selectedVar)
            break
        end
    end

    return selectedVar
end

"""
    model_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

return the SeaPearl model for to the N-Queens problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent (without solving it)

# Arguments
- `board_size::Int`: dimension of the board
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}() can be set to MostCenteredVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function model_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)

    rows = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows[i] = SeaPearl.IntVar(1, board_size, "row_"*string(i), trailer)
        SeaPearl.addVariable!(model, rows[i]; branchable=true)
    end

    rows_plus = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows_plus[i] = SeaPearl.IntVarViewOffset(rows[i], i, rows[i].id*"+"*string(i))
    end

    rows_minus = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows_minus[i] = SeaPearl.IntVarViewOffset(rows[i], -i, rows[i].id*"-"*string(i))
    end

    push!(model.constraints, SeaPearl.AllDifferent(rows, trailer))
    push!(model.constraints, SeaPearl.AllDifferent(rows_plus, trailer))
    push!(model.constraints, SeaPearl.AllDifferent(rows_minus, trailer))

    return model
end

"""
    solve_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

Solve the SeaPearl model for to the N-Queens problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent, and the function model_queens.

# Arguments
- `board_size::Int`: dimension of the board
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function solve_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    model = model_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    solve_queens(model::SeaPearl.CPModel)

Solve the SeaPearl model for to the N-Queens problem, using an existing model

# Arguments
- `model::SeaPearl.CPModel`: model (from model_queens)
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}().
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic().
"""
function solve_queens(model::SeaPearl.CPModel; variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    outputFromSeaPearl(model::SeaPearl.CPModel; optimality=false)

Shows the results of the N-Queens problem as type OutputDataQueens.

# Arguments
- `model::SeaPearl.CPModel`: needs the model to be already solved (by solve_queens)
"""
function outputFromSeaPearl(model::SeaPearl.CPModel; optimality=false)
    solutions = model.statistics.solutions
    nb_sols = length(solutions)
    n = length(model.variables)
    indices = Matrix{Int}(undef, n, nb_sols)
    for (key,sol) in solutions
        for i in 1:n
            indices[i, key]= sol["row_"*string(i)]
        end
    end
    return OutputDataQueens(nb_sols, indices)
end

"""
    print_queens(model::SeaPearl.CPModel; nb_sols=typemax(Int))

Print at max nb_sols solutions to the N-queens problems.

# Arguments
- `model::SeaPearl.CPModel`: needs the model to be already solved (by solve_queens)
- 'nb_sols::Int' : maximum number of solutions to print
"""
function print_queens(model::SeaPearl.CPModel; nb_sols=typemax(Int))
    variables = model.variables
    solutions = model.statistics.solutions
    board_size = length(model.variables)
    count = 0
    println("The solver found "*string(length(solutions))*" solutions to the "*string(board_size)*"-queens problem. Let's show them.")
    println()
    for key in keys(solutions)
        if(count >= nb_sols)
            break
        end
        println("Solution "*string(count+1))
        count +=1
        sol = solutions[key]
        for i in 1:board_size
            ind_queen = sol["row_"*string(i)]
            for j in 1:board_size
                if (j==ind_queen) print("Q ") else print("_ ") end
            end
            println()
        end
        println()
    end
end

"""
    print_queens(board_size::Int; nb_sols=typemax(Int), benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

Print at max nb_sols solutions to the N-queens problems, N givern as the board_size entry.

# Arguments
- `board_size::Int`: dimension of the board
- 'nb_sols::Int' : maximum number of solutions to print
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function print_queens(board_size::Int; nb_sols=typemax(Int), benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    model = solve_queens(board_size; variableSelection=variableSelection, valueSelection=valueSelection)
    variables = model.variables
    solutions = model.statistics.solutions
    count = 0
    println("The solver found "*string(length(solutions))*" solutions to the "*string(board_size)*"-queens problem. Let's show them.")
    println()
    for key in keys(solutions)
        if(count >= nb_sols)
            break
        end
        println("Solution "*string(count+1))
        count +=1
        sol = solutions[key]
        for i in 1:board_size
            ind_queen = sol["row_"*string(i)]
            for j in 1:board_size
                if (j==ind_queen) print("Q ") else print("_ ") end
            end
            println()
        end
        println()
    end
end

"""
    nb_solutions_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())::Int

return the number of solutions for to the N-Queens problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent.

# Arguments
- `board_size::Int`: dimension of the board
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function nb_solutions_queens(board_size::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())::Int
    return(length(solve_queens(board_size; variableSelection=variableSelection, valueSelection=valueSelection).statistics.solutions))
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
# Benchmark the counting of alanine residues in a PDB file

using MIToS.PDB

pdb_filepath = "pdbs/1AKE.pdb"
struc = read(pdb_filepath, PDBFile)

count() = length(findobjects(struc, Is(:name, "ALA")))

# Run to JIT compile
count()

elapsed = @elapsed count()

println(elapsed)

max_stars_repo_name: chemlove/pdb-benchmarks

max_issues_repo_name: chemlove/pdb-benchmarks

max_forks_repo_name: chemlove/pdb-benchmarks


text:
module KentFinancialEconomics

using HTTP
using CSV
using JSON
using DataFrames
using TimeSeries
using ZipFile
using Dates

export fundamentals, events, metadata, stockprices, fundprices

include("data_sources/quandl.jl")

end

max_stars_repo_name: alfredjmduncan/KentFinancialEconomics

max_issues_repo_name: alfredjmduncan/KentFinancialEconomics

max_forks_repo_name: alfredjmduncan/KentFinancialEconomics


text:
using DiffEqSensitivity, OrdinaryDiffEq, Zygote, Test

function get_param(breakpoints, values, t)
    for (i, tᵢ) in enumerate(breakpoints)
        if t <= tᵢ
            return values[i]
        end
    end
    return values[end]
end

function fiip(du, u, p, t)
    a = get_param([1., 2., 3.], p[1:4],  t)
    du[1] = dx =  a * u[1] - u[1] * u[2]
    du[2] = dy = -a * u[2] + u[1] * u[2]
end

p = [1., 1., 1., 1.]; u0 = [1.0;1.0]
prob = ODEProblem(fiip, u0, (0.0, 4.0), p);

dp1 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, sensealg = ForwardDiffSensitivity(), saveat = 0.1, abstol=1e-12, reltol=1e-12)), p)
dp2 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, sensealg = ForwardDiffSensitivity(convert_tspan=true), saveat = 0.1, abstol=1e-12, reltol=1e-12)), p)
dp3 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, sensealg = ForwardSensitivity(), saveat = 0.1, abstol=1e-12, reltol=1e-12)), p)
dp4 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, saveat = 0.1, abstol=1e-12, reltol=1e-12)), p)
dp5 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, saveat = 0.1, abstol=1e-12, reltol=1e-12, sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))), p)
dp6 = Zygote.gradient(p->sum(solve(prob, Tsit5(), u0=u0, p=p, saveat = 0.1, abstol=1e-12, reltol=1e-12, sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))), p)

@test dp1[1] ≈ dp2[1]
@test dp1[1] ≈ dp3[1]
@test dp1[1] ≈ dp4[1]
@test dp1[1] ≈ dp5[1]
@test sum(dp5[1]) ≈ sum(dp6[1])
@test all(dp6[1][1:3] .== 0)


text:
struct AsymptoticBlockSKCETest{K<:Kernel,E,S,Z} <: HypothesisTests.ZTest
    """Kernel."""
    kernel::K
    """Number of observations per block."""
    blocksize::Int
    """Number of blocks of observations."""
    nblocks::Int
    """Calibration error estimate (average of evaluations of blocks of observations)."""
    estimate::E
    """Standard error of evaluations of blocks of observations."""
    stderr::S
    """z-statistic."""
    z::Z
end

function AsymptoticBlockSKCETest(
    kernel::Kernel, blocksize::Int, predictions::AbstractVector, targets::AbstractVector
)
    # obtain number of samples
    nsamples = length(predictions)
    nsamples ≥ blocksize || error("there must be at least ", blocksize, " samples")

    # compute number of blocks
    nblocks = nsamples ÷ blocksize
    nblocks ≥ 2 || error("there must be at least 2 blocks")

    # evaluate U-statistic of the first block
    istart = 1
    iend = blocksize
    x = UnbiasedSKCE(kernel)(view(predictions, istart:iend), view(targets, istart:iend))

    # initialize the estimate and the sum of squares
    estimate = x / 1
    S = zero(estimate)^2

    # for all other blocks
    for b in 2:nblocks
        # evaluate U-statistic
        istart += blocksize
        iend += blocksize
        x = UnbiasedSKCE(kernel)(view(predictions, istart:iend), view(targets, istart:iend))

        # update the estimate
        Δestimate = x - estimate
        estimate += Δestimate / b
        S += Δestimate * (x - estimate)
    end

    # compute standard error and z-statistic
    stderr = sqrt(S / (nblocks * (nblocks - 1)))
    z = estimate / stderr

    return AsymptoticBlockSKCETest(kernel, blocksize, nblocks, estimate, stderr, z)
end

# HypothesisTests interface

HypothesisTests.default_tail(::AsymptoticBlockSKCETest) = :right

## have to specify and check keyword arguments in `pvalue` and `confint` to
## force `tail = :right` due to the default implementation in HypothesisTests

function HypothesisTests.pvalue(test::AsymptoticBlockSKCETest; tail=:right)
    if tail === :right
        StatsFuns.normccdf(test.z)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

# confidence interval by inversion
function HypothesisTests.confint(test::AsymptoticBlockSKCETest; level=0.95, tail=:right)
    HypothesisTests.check_level(level)

    if tail === :right
        q = StatsFuns.norminvcdf(level)
        lowerbound = test.estimate - q * test.stderr
        (max(zero(lowerbound), lowerbound), oftype(lowerbound, Inf))
    else
        throw(ArgumentError("tail = $(tail) is invalid"))
    end
end

HypothesisTests.testname(test::AsymptoticBlockSKCETest) = "Asymptotic block SKCE test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::AsymptoticBlockSKCETest)
    return "SKCE", zero(test.estimate), test.estimate
end

function HypothesisTests.show_params(io::IO, test::AsymptoticBlockSKCETest, ident="")
    println(io, ident, "number of observations per block: ", test.blocksize)
    println(io, ident, "number of blocks of observations: ", test.nblocks)
    println(io, ident, "z-statistic: ", test.z)
    return println(
        io, ident, "standard error of evaluations of pairs of observations: ", test.stderr
    )
end


text:
using DemoPackageMCF
using Test

@testset "DemoPackageMCF.jl" begin
    # Write your tests here.
end


text:
using Random, Test, HDF5

import HDF5.datatype

struct Simple
    a::Float64
    b::Int
end

function datatype(::Type{Simple})
    dtype = HDF5.API.h5t_create(HDF5.API.H5T_COMPOUND, sizeof(Simple))
    HDF5.API.h5t_insert(dtype, "a", fieldoffset(Simple, 1), datatype(Float64))
    HDF5.API.h5t_insert(dtype, "b", fieldoffset(Simple, 2), datatype(Int))
    HDF5.Datatype(dtype)
end

@testset "custom" begin
    N = 5
    v = [Simple(rand(Float64), rand(Int)) for i in 1:N, j in 1:N]

    fn = tempname()
    h5open(fn, "w") do h5f
        dtype = datatype(Simple)
        dspace = dataspace(v)
        dset = create_dataset(h5f, "data", dtype, dspace)
        write_dataset(dset, dtype, v)
    end

    h5open(fn, "r") do h5f
        dset = h5f["data"]
        @test_throws ErrorException read(dset, Float64)
        @test_throws ErrorException read(dset, Union{Float64,Int})

        v_read = read(dset, Simple)
        @test v_read == v

        indices = (1, 1:2:4)
        v_read = read(dset, Simple, indices...)
        @test v_read == v[indices...]

        v_read = read(h5f, "data"=>Simple)
        @test v_read == v
    end

    v_read = h5read(fn, "data"=>Simple)
    @test v_read == v
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

"""
# Vectors
"""
has_subvector(vec::Vector, sub::Vector) =
  any([sub == vec[i:(i+length(sub)-1)] for i=1:(length(vec)-length(sub)+1)])

"""
# Arrays
"""
dim(value::Any) = length(size(value))
function reduce(Val::Vector,Prb::Vector)
    idx = sortperm(Val)
    sPrb, sVal = Prb[idx], Val[idx]
    val = unique(sVal)
    prb = [sum(sPrb[sorted_range(sVal,nv)]) for nv in val]

    return val, prb
end
function reduce(Val::Array, Prb::Vector)
    idx = sortperm(collect(eachrow(Val)))
    sVal, sPrb = Val[idx,:], Prb[idx]
    val = unique(sVal, dims=1)
    cVal = collect(eachrow(sVal))
    prb = [sum(sPrb[sorted_range(cVal,nv)]) for nv in eachrow(val)]

    return val, prb
end
sorted_range(sVal,nv) = searchsortedfirst(sVal,nv):searchsortedlast(sVal,nv)

"""
# KWARGS
"""
# functions
isa_svr(kwargs::Iterators.Pairs) =
    prod([isa(nv,Single) || isa(nv,Vector) || isa(nv,UnitRange)
        for nv in values(kwargs)])
isa_sm(kwargs::Iterators.Pairs) =
    prod([isa(nv,Single) || isa(nv,Matrix) for nv in values(kwargs)])
function equal_size(kwargs::Iterators.Pairs)
    values_length = [size(value) for value in values(kwargs)
                                 if !isa(value,Single)]
    empty = isempty(values_length)
    equal = all(x->x==first(values_length),values_length)
    # test
    empty || equal || return false
end
function test(kwargs::Iterators.Pairs)
    (isa_svr(kwargs) || isa_sm(kwargs)) || return false                         # TODO Memento logger
    equal_size(kwargs) || return false                                          # TODO Memento logger
end

indices_of(kwargs::Iterators.Pairs) =
    CartesianIndices(kwargs[findfirst(x->!isa(x,Single),values(kwargs))])
reduce(kwargs::Iterators.Pairs, nc; excl=[]) =
    PropDict(key => isa(value,Single) ? value : value[nc] for (key,value) in kwargs
                                                          if !in(key,excl))

"""
# UNITS
"""
Base.log(a::Quantity{T,D,U}) where T <: Real where D where U = log(a.val)unit(a)
Base.exp(a::Quantity{T,D,U}) where T <: Real where D where U = exp(a.val)unit(a)


text:
"""
Provides a (possibly incomplete) description of a permutation representation.
"""
type PermutationRepresentationDefinition{T} <: RepresentationDefinition
  alphabet::Vector{T}
end

immutable PermutationRepresentation{T} <: Representation
  """
  The alphabet of values (which may take any form, as specified by the type parameter, T)
  whose permutations are to be searched.
  """
  alphabet::Vector{T}

  """
  The number of values in the alphabet for this representation.
  """
  size::Int
end

"""
Composes a permutation representation from a provided definition.
"""
compose!{T}(def::PermutationRepresentationDefinition{T}) =
  PermutationRepresentation{T}(def.alphabet, length(def.alphabet))

"""
Builds a permutation representation using an alphabet defined by a given
unit range of values.
"""
permutation(range::UnitRange) =
  permutation(collect(range))

"""
Builds a permutation representation using an alphabet supplied as a vector
(i.e. a one-dimensional array) of values.
"""
permutation{T}(alphabet::Vector{T}) =
  PermutationRepresentationDefinition(alphabet)

"""
Returns the chromosome used by this permutation representation (which is always
implemented as a vector of some type, T, given by the permutation template
parameters).
"""
chromosome{T}(rep::PermutationRepresentation{T}) = Vector{T}

"""
Returns the codon type used by this permutation, given by its template
parameter, T.
"""
codon_type{T}(rep::PermutationRepresentation{T}) = T

#"""
#Produces a human-readable string representation of an instance of this
#permutation representation.
#"""
#describe{T}(::PermutationRepresentationDefinition{T}, ins::Vector{T}) =

"""
Produces a random permutation using a provided permutation representation.
"""
generate{T}(r::PermutationRepresentation{T}) =
  shuffle(r.alphabet)


text:
module MinimaxCalibratedEBayes

using Reexport

@reexport using StatsBase,
                Distributions,
                EBayes

using ApproxFun:Interval, Fun
using DiffResults
using Expectations
using ExponentialFamilies
using ForwardDiff
using JuMP
using KernelDensity
using LaTeXStrings
using LinearAlgebra
using MathOptInterface
using Optim
using OrthogonalPolynomialsQuasi
using OscillatoryIntegrals
using Parameters
using Plots
using RecipesBase
using Roots
using Setfield
using SpecialFunctions:erfi
using QuadGK

import StatsBase:confint,
				 fit,
                 Histogram,
                 binindex,
                 midpoints,
				 nobs,
                 response

import Statistics:var

import Distributions:cf,
                     estimate,
                     pdf,
                     location,
					 loglikelihood,
					 support,
					 DiscreteNonParametric

import Expectations:expectation

import Base:eltype,
            extrema,
			step,
            first,
            last,
            length,
            zero,
			getindex,
			lastindex

import Base.Broadcast: broadcastable

include("bias_adjusted_ci.jl")
include("marginal_binning.jl")
include("inference_targets.jl")
include("normal_rules.jl")
include("marginal_kde.jl")
include("logspline_g_new.jl")
include("butucea_comte.jl")
include("prior_convex_class.jl")
include("sinkhorn.jl")
include("npmle.jl")
include("hermite.jl")
include("main_mceb.jl")
include("load_datasets.jl")
include("helper_plots.jl")

export MCEBHistogram,
	   DiscretizedStandardNormalSamples,
       DiscretizedAffineEstimator,
       marginalize,
       SincKernel,
       DeLaValleePoussinKernel,
       certainty_banded_KDE,
	   KDEInfinityBandOptions,
	   KDEInfinityBand,
       EBayesTarget,
       MarginalDensityTarget,
       PriorDensityTarget,
	   PriorTailProbability,
       riesz_representer,
	   PosteriorTarget,
       PosteriorMeanNumerator,
       LFSRNumerator,
       PosteriorMean,
       LFSR,
       ButuceaComte,
       GaussianMixturePriorClass,
	   HermitePriorClass,
       worst_case_bias,
       SteinMinimaxEstimator,
       steinminimaxplot,
       bias_adjusted_gaussian_ci,
	   confint,
	   confintplot,
	   confintplot!,
	   set_neighborhood,
	   MinimaxCalibratorSetup,
	   ExponentialFamilyDeconvolutionMLE,
	   MinimaxCalibratorOptions

end # module


text:
"""
    stack(df::AbstractDataFrame, [measure_vars], [id_vars];
          variable_name=:variable, value_name=:value,
          view::Bool=false, variable_eltype::Type=String)

Stack a data frame `df`, i.e. convert it from wide to long format.

Return the long-format `DataFrame` with: columns for each of the `id_vars`,
column `value_name` (`:value` by default)
holding the values of the stacked columns (`measure_vars`), and
column `variable_name` (`:variable` by default) a vector holding
the name of the corresponding `measure_vars` variable.

If `view=true` then return a stacked view of a data frame (long format).
The result is a view because the columns are special `AbstractVectors`
that return views into the original data frame.

# Arguments
- `df` : the AbstractDataFrame to be stacked
- `measure_vars` : the columns to be stacked (the measurement variables),
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  If neither `measure_vars` or `id_vars` are given, `measure_vars`
  defaults to all floating point columns.
- `id_vars` : the identifier columns that are repeated during stacking,
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  Defaults to all variables that are not `measure_vars`
- `variable_name` : the name (`Symbol` or string) of the new stacked column that
  shall hold the names of each of `measure_vars`
- `value_name` : the name (`Symbol` or string) of the new stacked column containing
  the values from each of `measure_vars`
- `view` : whether the stacked data frame should be a view rather than contain
  freshly allocated vectors.
- `variable_eltype` : determines the element type of column `variable_name`.
  By default a `PooledArray{String}` is created.
  If `variable_eltype=Symbol` a `PooledVector{Symbol}` is created,
  and if `variable_eltype=CategoricalValue{String}`
  a `CategoricalArray{String}` is produced (call `using CategoricalArrays` first if needed)
  Passing any other type `T` will produce a `PooledVector{T}` column
  as long as it supports conversion from `String`.
  When `view=true`, a `RepeatedVector{T}` is produced.

# Examples
```julia
d1 = DataFrame(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stack(d1, [:c, :d])
d1s2 = stack(d1, [:c, :d], [:a])
d1m = stack(d1, Not([:a, :b, :e]))
d1s_name = stack(d1, Not([:a, :b, :e]), variable_name=:somemeasure)
```
"""
function stack(df::AbstractDataFrame,
               measure_vars = findall(col -> eltype(col) <: Union{AbstractFloat, Missing},
                                      eachcol(df)),
               id_vars = Not(measure_vars);
               variable_name::SymbolOrString=:variable,
               value_name::SymbolOrString=:value, view::Bool=false,
               variable_eltype::Type=String)
    variable_name_s = Symbol(variable_name)
    value_name_s = Symbol(value_name)
    # getindex from index returns either Int or AbstractVector{Int}
    mv_tmp = index(df)[measure_vars]
    ints_measure_vars = mv_tmp isa Int ? [mv_tmp] : mv_tmp
    idv_tmp = index(df)[id_vars]
    ints_id_vars = idv_tmp isa Int ? [idv_tmp] : idv_tmp
    if view
        return _stackview(df, ints_measure_vars, ints_id_vars,
                          variable_name=variable_name_s,
                          value_name=value_name_s,
                          variable_eltype=variable_eltype)
    end
    N = length(ints_measure_vars)
    cnames = _names(df)[ints_id_vars]
    push!(cnames, variable_name_s)
    push!(cnames, value_name_s)
    if variable_eltype === Symbol
        catnms = PooledArray(_names(df)[ints_measure_vars])
    elseif variable_eltype === String
        catnms = PooledArray(names(df, ints_measure_vars))
    else
        # this covers CategoricalArray{String} in particular
        # (note that copyto! inserts levels in their order of appearance)
        nms = names(df, ints_measure_vars)
        simnms = similar(nms, variable_eltype)
        catnms = simnms isa Vector ? PooledArray(catnms) : simnms
        copyto!(catnms, nms)
    end
    return DataFrame(AbstractVector[[repeat(df[!, c], outer=N) for c in ints_id_vars]..., # id_var columns
                                    repeat(catnms, inner=nrow(df)),                       # variable
                                    vcat([df[!, c] for c in ints_measure_vars]...)],      # value
                     cnames, copycols=false)
end

function _stackview(df::AbstractDataFrame, measure_vars::AbstractVector{Int},
                    id_vars::AbstractVector{Int}; variable_name::Symbol,
                    value_name::Symbol, variable_eltype::Type)
    N = length(measure_vars)
    cnames = _names(df)[id_vars]
    push!(cnames, variable_name)
    push!(cnames, value_name)
    if variable_eltype === Symbol
        catnms = _names(df)[measure_vars]
    elseif variable_eltype === String
        catnms = names(df, measure_vars)
    else
        # this covers CategoricalArray{String} in particular,
        # as copyto! inserts levels in their order of appearance
        nms = names(df, measure_vars)
        catnms = copyto!(similar(nms, variable_eltype), nms)
    end
    return DataFrame(AbstractVector[[RepeatedVector(df[!, c], 1, N) for c in id_vars]..., # id_var columns
                                    RepeatedVector(catnms, nrow(df), 1),                  # variable
                                    StackedVector(Any[df[!, c] for c in measure_vars])],  # value
                     cnames, copycols=false)
end

"""
    unstack(df::AbstractDataFrame, rowkeys, colkey, value; renamecols::Function=identity)
    unstack(df::AbstractDataFrame, colkey, value; renamecols::Function=identity)
    unstack(df::AbstractDataFrame; renamecols::Function=identity)

Unstack data frame `df`, i.e. convert it from long to wide format.

If `colkey` contains `missing` values then they will be skipped and a warning
will be printed.

If combination of `rowkeys` and `colkey` contains duplicate entries then last
`value` will be retained and a warning will be printed.

# Arguments
- `df` : the AbstractDataFrame to be unstacked
- `rowkeys` : the columns with a unique key for each row, if not given,
  find a key by grouping on anything not a `colkey` or `value`.
  Can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
- `colkey` : the column ($COLUMNINDEX_STR) holding the column names in wide format,
  defaults to `:variable`
- `value` : the value column ($COLUMNINDEX_STR), defaults to `:value`
- `renamecols` : a function called on each unique value in `colkey` which must
                 return the name of the column to be created (typically as a string
                 or a `Symbol`). Duplicate names are not allowed.

# Examples
```julia
wide = DataFrame(id = 1:12,
                 a  = repeat([1:3;], inner = [4]),
                 b  = repeat([1:4;], inner = [3]),
                 c  = randn(12),
                 d  = randn(12))

long = stack(wide)
wide0 = unstack(long)
wide1 = unstack(long, :variable, :value)
wide2 = unstack(long, :id, :variable, :value)
wide3 = unstack(long, [:id, :a], :variable, :value)
wide4 = unstack(long, :id, :variable, :value, renamecols=x->Symbol(:_, x))
```
Note that there are some differences between the widened results above.
"""
function unstack(df::AbstractDataFrame, rowkey::ColumnIndex, colkey::ColumnIndex,
                 value::ColumnIndex; renamecols::Function=identity)
    refkeycol = categorical(df[!, rowkey])
    droplevels!(refkeycol)
    keycol = categorical(df[!, colkey])
    droplevels!(keycol)
    valuecol = df[!, value]
    return _unstack(df, index(df)[rowkey], index(df)[colkey],
                    keycol, valuecol, refkeycol, renamecols)
end

function _unstack(df::AbstractDataFrame, rowkey::Int, colkey::Int,
                  keycol::CategoricalVector, valuecol::AbstractVector,
                  refkeycol::CategoricalVector, renamecols::Function)
    Nrow = length(refkeycol.pool)
    Ncol = length(keycol.pool)
    unstacked_val = [similar_missing(valuecol, Nrow) for i in 1:Ncol]
    hadmissing = false # have we encountered missing in refkeycol
    mask_filled = falses(Nrow+1, Ncol) # has a given [row,col] entry been filled?
    warned_dup = false # have we already printed duplicate entries warning?
    warned_missing = false # have we already printed missing in keycol warning?
    for k in 1:nrow(df)
        kref = keycol.refs[k]
        if kref <= 0 # we have found missing in colkey
            if !warned_missing
                @warn("Missing value in variable :$(_names(df)[colkey]) at row $k. Skipping.")
                warned_missing = true
            end
            continue # skip processing it
        end
        refkref = refkeycol.refs[k]
        if refkref <= 0 # we have found missing in rowkey
            if !hadmissing # if it is the first time we have to add a new row
                hadmissing = true
                # we use the fact that missing is greater than anything
                for i in eachindex(unstacked_val)
                    push!(unstacked_val[i], missing)
                end
            end
            i = length(unstacked_val[1])
        else
            i = refkref
        end
        if !warned_dup && mask_filled[i, kref]
            @warn("Duplicate entries in unstack at row $k for key "*
                  "$(refkeycol[k]) and variable $(keycol[k]).")
            warned_dup = true
        end
        unstacked_val[kref][i] = valuecol[k]
        mask_filled[i, kref] = true
    end
    levs = levels(refkeycol)
    # we have to handle a case with missings in refkeycol as levs will skip missing
    col = similar(df[!, rowkey], length(levs) + hadmissing)
    copyto!(col, levs)
    hadmissing && (col[end] = missing)
    df2 = DataFrame(unstacked_val, Symbol.(renamecols.(levels(keycol))), copycols=false)
    return insertcols!(df2, 1, _names(df)[rowkey] => col)
end

function unstack(df::AbstractDataFrame, rowkeys, colkey::ColumnIndex,
                 value::ColumnIndex; renamecols::Function=identity)
    rowkey_ints = index(df)[rowkeys]
    @assert rowkey_ints isa AbstractVector{Int}
    length(rowkey_ints) == 0 && throw(ArgumentError("No key column found"))
    length(rowkey_ints) == 1 && return unstack(df, rowkey_ints[1], colkey, value,
                                               renamecols=renamecols)
    g = groupby(df, rowkey_ints, sort=true)
    keycol = categorical(df[!, colkey])
    droplevels!(keycol)
    valuecol = df[!, value]
    return _unstack(df, rowkey_ints, index(df)[colkey], keycol, valuecol, g, renamecols)
end

function unstack(df::AbstractDataFrame, colkey::ColumnIndex, value::ColumnIndex;
                 renamecols::Function=identity)
    colkey_int = index(df)[colkey]
    value_int = index(df)[value]
    return unstack(df, Not(colkey_int, value_int), colkey_int, value_int,
            renamecols=renamecols)
end

unstack(df::AbstractDataFrame; renamecols::Function=identity) =
    unstack(df, :variable, :value, renamecols=renamecols)

function _unstack(df::AbstractDataFrame, rowkeys::AbstractVector{Int},
                  colkey::Int, keycol::CategoricalVector,
                  valuecol::AbstractVector, g::GroupedDataFrame,
                  renamecols::Function)
    idx, starts, ends = g.idx, g.starts, g.ends
    groupidxs = [idx[starts[i]:ends[i]] for i in 1:length(starts)]
    rowkey = zeros(Int, size(df, 1))
    for i in 1:length(groupidxs)
        rowkey[groupidxs[i]] .= i
    end
    df1 = df[idx[starts], g.cols]
    Nrow = length(g)
    Ncol = length(levels(keycol))
    unstacked_val = [similar_missing(valuecol, Nrow) for i in 1:Ncol]
    mask_filled = falses(Nrow, Ncol)
    warned_dup = false
    warned_missing = false
    for k in 1:nrow(df)
        kref = keycol.refs[k]
        if kref <= 0
            if !warned_missing
                @warn("Missing value in variable :$(_names(df)[colkey]) at row $k. Skipping.")
                warned_missing = true
            end
            continue
        end
        i = rowkey[k]
        if !warned_dup && mask_filled[i, kref]
            @warn("Duplicate entries in unstack at row $k for key "*
                 "$(tuple((df[k,s] for s in rowkeys)...)) and variable $(keycol[k]).")
            warned_dup = true
        end
        unstacked_val[kref][i] = valuecol[k]
        mask_filled[i, kref] = true
    end
    df2 = DataFrame(unstacked_val, Symbol.(renamecols.(levels(keycol))), copycols=false)
    hcat(df1, df2, copycols=false)
end

"""
    StackedVector <: AbstractVector

An `AbstractVector` that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

# Constructor
```julia
StackedVector(d::AbstractVector)
```

# Arguments
- `d...` : one or more AbstractVectors

# Examples
```julia
StackedVector(Any[[1,2], [9,10], [11,12]])  # [1,2,9,10,11,12]
```
"""
struct StackedVector{T} <: AbstractVector{T}
    components::Vector{Any}
end

StackedVector(d::AbstractVector) =
    StackedVector{promote_type(map(eltype, d)...)}(d)

function Base.getindex(v::StackedVector{T}, i::Int)::T where T
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    return v.components[j][k]
end

Base.IndexStyle(::Type{StackedVector}) = Base.IndexLinear()
Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.eltype(v::Type{StackedVector{T}}) where {T} = T
Base.similar(v::StackedVector, T::Type, dims::Union{Integer, AbstractUnitRange}...) =
    similar(v.components[1], T, dims...)

CategoricalArrays.CategoricalArray(v::StackedVector) =
    CategoricalArray(v[:]) # could be more efficient

"""
    RepeatedVector{T} <: AbstractVector{T}

An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

# Constructor
```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

# Arguments
- `parent` : the AbstractVector that's repeated
- `inner` : the numer of times each element is repeated
- `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

# Examples
```julia
RepeatedVector([1,2], 3, 1)   # [1,1,1,2,2,2]
RepeatedVector([1,2], 1, 3)   # [1,2,1,2,1,2]
RepeatedVector([1,2], 2, 2)   # [1,2,1,2,1,2,1,2]
```
"""
struct RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

Base.parent(v::RepeatedVector) = v.parent
DataAPI.levels(v::RepeatedVector) = levels(parent(v))
CategoricalArrays.isordered(v::RepeatedVector{<:Union{CategoricalValue, Missing}}) =
    isordered(parent(v))

function Base.getindex(v::RepeatedVector, i::Int)
    N = length(parent(v))
    idx = Base.fld1(mod1(i,v.inner*N),v.inner)
    parent(v)[idx]
end

Base.IndexStyle(::Type{<:RepeatedVector}) = Base.IndexLinear()
Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(parent(v))
Base.eltype(v::Type{RepeatedVector{T}}) where {T} = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(parent(v)), v.inner, v.outer)
Base.similar(v::RepeatedVector, T::Type, dims::Dims) = similar(parent(v), T, dims)
Base.unique(v::RepeatedVector) = unique(parent(v))

function CategoricalArrays.CategoricalArray(v::RepeatedVector)
    res = CategoricalArray(parent(v), levels=levels(parent(v)))
    res.refs = repeat(res.refs, inner = [v.inner], outer = [v.outer])
    res
end


text:
@testset "rounding mode" begin
  @test ArbNumerics.match_rounding_mode(RoundNearest) == 4
end


text:
@testset "basic MortalityTable" begin
    @testset "basic structure" begin
        q1 = UltimateMortality([i for i = 0:19], start_age = 0)
        @test q1[0] == 0
        @test q1[1] == 1
        @test q1[0:1] == [0, 1]
        @test omega(q1) == 19
        @test_throws BoundsError q1[omega(q1)+1]

        # non-zero start age
        q2 = UltimateMortality([i for i = 0:9], start_age = 5)
        @test_throws BoundsError q2[4]
        @test q2[5] == 0

        # select strucutre
        s = [ia + d for ia = 0:5, d = 0:4]

        q3 = SelectMortality(s, q1, start_age = 0)
        @test q3[0][0] == 0
        @test q3[0][0:1] == [0, 1]
        @test omega(q3[0]) == 19
        @test_throws BoundsError q3[omega(q3[0])+1]
        @test q3[0][19] == 19

        @test q3[5][5] == 5

        mt1 = MortalityTable(q3, q1)

        @test mt1.select[0][1] == 1
        @test mt1.ultimate[1] == 1

        mt2 = MortalityTable(q1)

        @test mt2.ultimate[0] == 0
        @test mt2[0] == 0

    end

    @testset "off-aligned select and ult" begin
        select_matrix = [(i + j - 1) / 100 for i = 0:10, j = 1:20]
        ult = UltimateMortality([i / 100 for i = 18:100], start_age = 18)

        select = SelectMortality(select_matrix, ult)

        for issue_age = 0:10
            for dur = 1:30
                @test select[issue_age][issue_age+dur-1] == (issue_age + dur - 1) / 100
            end
        end
    end

    # test time zero accumlated force
    @testset "accumulated force" begin
        q4 = UltimateMortality([0.1, 0.3, 0.6, 1])

        @test survival(q4, 0) ≈ 1
        @test decrement(q4, 0) ≈ 0

        @test survival(q4, 1) ≈ 0.9
        @test decrement(q4, 1) ≈ 0.1

        @test survival(q4, 1, 1) ≈ 1.0
        @test survival(q4, 1, 2) ≈ 0.7
        @test decrement(q4, 1, 1) ≈ 0.0
        @test decrement(q4, 1, 2) ≈ 0.3

        @test survival(q4, 1, 4) ≈ 0.0
        @test decrement(q4, 1, 4) ≈ 1.0

        @test survival(q4, -1) ≈ 1.0
        @test survival(q4, 4, -1) ≈ 1.0
        @test decrement(q4, -1) ≈ 0.0
        @test decrement(q4, 4, -1) ≈ 0.0
    end

    @testset "Metadata" begin
        d = TableMetaData()

        @test isnothing(d.name)

        d = TableMetaData(name = "test")
        @test d.name == "test"
    end

    @testset "mortality_vector" begin
        v = [i for i = 3:10]
        q = mortality_vector(v, start_age = 3)
        @test q[3] == 3
        @test q[10] == 10

        q = mortality_vector(collect(0:5))
        @test q[0] == 0
        @test q[5] == 5
    end

    @testset "utility functions" begin
        @test MortalityTables._decrement(1.0, 0.05) ≈ 1.0 * (1 - 0.05)
    end

end


text:
using WeavePynb
using Mustache

mmd(fname) = mmd_to_html(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")
## uncomment to generate just .md files
mmd(fname) = mmd_to_md(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")

fnames = [
          "odes",
          "euler"
          ]

function process_file(nm, twice=false)
    include("$nm.jl")
    mmd_to_md("$nm.mmd")
    markdownToHTML("$nm.md")
    twice && markdownToHTML("$nm.md")
end

process_files(twice=false) = [process_file(nm, twice) for nm in fnames]

"""
## TODO ODEs

"""


text:
#######################################################
"""
Copyright 2021 Alliance for Sustainable Energy and other
NAERM Project Developers. See the top-level COPYRIGHT file for details.

Author: Surya Chandan Dhulipala
Email: suryachandan.dhulipala@nrel.gov
"""
# November 2021
# EGRET --> SIIP Linkage Module
# EGRET JSON --> CSV files formatted according to SIIP Tabular Data Parser Requirements --> SIIP PSY System
#######################################################
module EGRET2SIIP
#################################################################################
# Exports
#################################################################################
export parse_EGRET_JSON
export parse_tabular_data
export EGRET_TO_PSY
#################################################################################
# Imports
#################################################################################
import JSON
import DataFrames
import CSV
import Dates
import PowerSystems
#################################################################################
# Includes
#################################################################################
include("parsers/egret_json_to_csv.jl")
include("main/parse_tabular_data.jl")
include("main/egret_to_psy.jl")
end


text:
#---------------------------------------------------------------------
#*	only liquid-phase
#*--------------------------------------------------------------------
type stoic_liq
	stoic_liq()=begin
		new(
			tank_liq(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of reactions",
				:Default=>1
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Stoichiometric matrix",
				:Symbol=>"\\nu"
			)),(NComp,NReac)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Overall component rate of reaction"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fractional conversion of component",
				:Symbol=>"X",
				:Default=>0
			)),(NComp)),
			[
				:(Outlet.F*Outlet.z = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z + rate*_base_1._base_1.Tank.V),
				:(Outlet.P = _base_1._base_1.Outletm.P),
				:(Outlet.F*Outlet.h = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.h),
				:(Outlet.F = _base_1._base_1.Outletm.F),
				:(Outlet.F*Outlet.z([1:_base_1._base_1.NComp]) = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z([1:_base_1._base_1.NComp])*(1 - conv([1:_base_1._base_1.NComp]))),
				:(conv([1:_base_1._base_1.NComp]) = 1),
				:(conv([1:_base_1._base_1.NComp]) = 0),
			],
			[
				"Outlet stream","Mechanical equilibrium","Energy balance","Steady-state","Molar conversion","Molar conversion","Molar conversion",
			],
			[:NReac,:stoic,],
			[:Outlet,:rate,:conv,]
		)
	end
	_base_1::tank_liq
	NReac::DanaInteger
	stoic::Array{DanaReal}
	Outlet::liquid_stream
	rate::Array{reaction_mol}
	conv::Array{DanaReal}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export stoic_liq
function setEquationFlow(in::stoic_liq)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	if (Outletm.z(i) > 1e-16)
		addEquation(5)
	else
		if (Outlet.z(i) > 0)
			addEquation(6)
			# ?

		else
			addEquation(7)
			# ?

		end
	end
end
function atributes(in::stoic_liq,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Basic model for a liquid-phase stoichiometric CSTR"
	fields[:Info]="
== Assumptions ==
* only liquid-phase
* steady-state
"
	drive!(fields,_)
	return fields
end
stoic_liq(_::Dict{Symbol,Any})=begin
	newModel=stoic_liq()
	newModel.attributes=atributes(newModel,_)
	newModel
end