text:
"""
    LeakyReLU{T, N}(α::T) <: Bijector{N}

Defines the invertible mapping

    x ↦ x if x ≥ 0 else αx

where α > 0.
"""
struct LeakyReLU{T, N} <: Bijector{N}
    α::T
end

LeakyReLU(α::T; dim::Val{N} = Val(0)) where {T<:Real, N} = LeakyReLU{T, N}(α)
LeakyReLU(α::T; dim::Val{N} = Val(D)) where {D, T<:AbstractArray{<:Real, D}, N} = LeakyReLU{T, N}(α)

up1(b::LeakyReLU{T, N}) where {T, N} = LeakyReLU{T, N + 1}(b.α)

# (N=0) Univariate case
function (b::LeakyReLU{<:Any, 0})(x::Real)
    mask = x < zero(x)
    return mask * b.α * x + !mask * x
end
(b::LeakyReLU{<:Any, 0})(x::AbstractVector{<:Real}) = map(b, x)

function Base.inv(b::LeakyReLU{<:Any,N}) where N
    invα = inv.(b.α)
    return LeakyReLU{typeof(invα),N}(invα)
end

function logabsdetjac(b::LeakyReLU{<:Any, 0}, x::Real)
    mask = x < zero(x)
    J = mask * b.α + (1 - mask) * one(x)
    return log(abs(J))
end
logabsdetjac(b::LeakyReLU{<:Real, 0}, x::AbstractVector{<:Real}) = map(x -> logabsdetjac(b, x), x)

# We implement `forward` by hand since we can re-use the computation of
# the Jacobian of the transformation. This will lead to faster sampling
# when using `rand` on a `TransformedDistribution` making use of `LeakyReLU`.
function forward(b::LeakyReLU{<:Any, 0}, x::Real)
    mask = x < zero(x)
    J = mask * b.α + !mask * one(x)
    return (rv=J * x, logabsdetjac=log(abs(J)))
end

# Batched version
function forward(b::LeakyReLU{<:Any, 0}, x::AbstractVector)
    J = let T = eltype(x), z = zero(T), o = one(T)
        @. (x < z) * b.α + (x > z) * o
    end
    return (rv=J .* x, logabsdetjac=log.(abs.(J)))
end

# (N=1) Multivariate case
function (b::LeakyReLU{<:Any, 1})(x::AbstractVecOrMat)
    return let z = zero(eltype(x))
        @. (x < z) * b.α * x + (x > z) * x
    end
end

function logabsdetjac(b::LeakyReLU{<:Any, 1}, x::AbstractVecOrMat)
    # Is really diagonal of jacobian
    J = let T = eltype(x), z = zero(T), o = one(T)
        @. (x < z) * b.α + (x > z) * o
    end

    if x isa AbstractVector
        return sum(log.(abs.(J)))
    elseif x isa AbstractMatrix
        return vec(sum(log.(abs.(J)); dims = 1))  # sum along column
    end
end

# We implement `forward` by hand since we can re-use the computation of
# the Jacobian of the transformation. This will lead to faster sampling
# when using `rand` on a `TransformedDistribution` making use of `LeakyReLU`.
function forward(b::LeakyReLU{<:Any, 1}, x::AbstractVecOrMat)
    # Is really diagonal of jacobian
    J = let T = eltype(x), z = zero(T), o = one(T)
        @. (x < z) * b.α + (x > z) * o
    end

    if x isa AbstractVector
        logjac = sum(log.(abs.(J)))
    elseif x isa AbstractMatrix
        logjac = vec(sum(log.(abs.(J)); dims = 1))  # sum along column
    end

    y = J .* x
    return (rv=y, logabsdetjac=logjac)
end


text:
const normalfam = let
    loghval = log(1/sqrt(2π)) :: AbstractFloat
    logh(x::Real) = loghval :: AbstractFloat
    t(x::Real) = (x,x^2)
    function a(η::Tuple{Real, Real})
        η₁,η₂ = η
        -η₁^2 / 4η₂ - 0.5 * log(-2η₂)
    end
    NaturalExponentialFamily{Tuple{Real, Real},Real}(logh,t,a)
end

export normal
normal(μ,σ) = normalfam((μ/σ^2, -1/2σ^2))

ext: jl


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

mutable struct ParameterValues <: SwaggerModel

    function ParameterValues(;)
        o = new()
        o
    end
end # type ParameterValues

const _property_map_ParameterValues = Dict{Symbol,Symbol}()
const _property_types_ParameterValues = Dict{Symbol,String}()
Base.propertynames(::Type{ ParameterValues }) = collect(keys(_property_map_ParameterValues))
Swagger.property_type(::Type{ ParameterValues }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ParameterValues[name]))}
Swagger.field_name(::Type{ ParameterValues }, property_name::Symbol) =  _property_map_ParameterValues[property_name]

function check_required(o::ParameterValues)
    true
end

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


text:
"""
    LeafModels(interception, energy, photosynthesis, stomatal_conductance, status)
    LeafModels(;interception = missing, energy = missing, photosynthesis = missing,
        stomatal_conductance = missing,status...)

LeafModels component, which is a subtype of `AbstractComponentModel` implementing a component with
a photosynthetic activity. It could be a leaf, or a leaflet, or whatever kind of component
that is photosynthetic. The name `LeafModels` was chosen not because it is generic, but because it
is short, simple and self-explanatory.

# Arguments

- `interception <: Union{Missing,AbstractInterceptionModel}`: An interception model.
- `energy <: Union{Missing,AbstractEnergyModel}`: An energy model, e.g. [`Monteith`](@ref).
- `photosynthesis <: Union{Missing,AbstractAModel}`: A photosynthesis model, e.g. [`Fvcb`](@ref)
- `stomatal_conductance <: Union{Missing,AbstractGsModel}`: A stomatal conductance model, e.g. [`Medlyn`](@ref) or
[`ConstantGs`](@ref)
- `status <: MutableNamedTuple`: a mutable named tuple to track the status (*i.e.* the variables) of
the leaf. Values are set to `0.0` if not provided as VarArgs (see examples)

# Details

The status field depends on the input models. You can get the variables needed by a model using
[`variables`](@ref) on the instantiation of a model. Generally the variables are:

## Light interception model  (see [`AbstractInterceptionModel`](@ref))

- `Rₛ` (W m-2): net shortwave radiation (PAR + NIR). Often computed from a light interception model
- `PPFD` (μmol m-2 s-1): absorbed Photosynthetic Photon Flux Density
- `skyFraction` (0-2): view factor between the object and the sky for both faces.

## Energy balance model (see [`energy_balance`](@ref))

- `Tₗ` (°C): temperature of the object. Often computed from an energy balance model.
- `Rₗₗ` (W m-2): net longwave radiation for the object (TIR)
- `Cₛ` (ppm): stomatal CO₂ concentration
- `ψₗ` (kPa): leaf water potential
- `H` (W m-2): sensible heat flux
- `λE` (W m-2): latent heat flux
- `Dₗ` (kPa): vapour pressure difference between the surface and the saturation vapour

## Photosynthesis model (see [`photosynthesis`](@ref))

- `A` (μmol m-2 s-1): carbon assimilation
- `Gbₕ` (m s-1): boundary conductance for heat (free + forced convection)
- `Cᵢ` (ppm): intercellular CO₂ concentration

# Examples

```julia
# A leaf with a width of 0.03 m, that uses the Monteith and Unsworth (2013) model for energy
# balance, The Farquhar et al. (1980) model for photosynthesis, and a constant stomatal
# conductance for CO₂ of 0.0011 with no residual conductance. The status of
# the leaf is not set yet, all are initialised at `0.0`:
LeafModels(energy = Monteith(),
     photosynthesis = Fvcb(),
     stomatal_conductance = ConstantGs(0.0, 0.0011))

# If we need to initialise some variables at different values, we can call the leaf as:

LeafModels(photosynthesis = Fvcb(),Cᵢ = 380.0)

# Or again:
LeafModels(photosynthesis = Fvcb(), energy = Monteith(), Cᵢ = 380.0, Tₗ = 20.0)
```
"""
struct LeafModels{I <: Union{Missing,AbstractInterceptionModel},E <: Union{Missing,AbstractEnergyModel},A <: Union{Missing,AbstractAModel},Gs <: Union{Missing,AbstractGsModel},S <: Union{MutableNamedTuple,Vector{MutableNamedTuple}}} <: AbstractComponentModel
    interception::I
    energy::E
    photosynthesis::A
    stomatal_conductance::Gs
    status::S
end

function LeafModels(;interception = missing, energy = missing,
                photosynthesis = missing, stomatal_conductance = missing,status...)
    status_vals = collect(values(status))

    length_status = [length(i) for i in status_vals]
    if any(length_status .> 1)
        # One of the variable is given as an array, meaning this is actually several time-steps
        # In this case we make an array of leaves
        max_length_st = maximum(length_status)
        for i in 1:length(status_vals)
            # If the ith status has length one, repeat its value to match the others:
            if length_status[i] == 1
                status_vals[i] = repeat([status_vals[i]], max_length_st)
            else
                length_status[i] != max_length_st && @error "$(keys(status)[i]) should be length $max_length_st or 1"
            end
        end

        # Making a status for each ith value in the user status:
        status_array = MutableNamedTuple[]
        for i in 1:max_length_st
            # status = init_variables_manual(interception, energy, photosynthesis,
            # stomatal_conductance;NamedTuple{keys(status)}([j[i] for j in status_vals])...)
            # push!(leaf_array, LeafModels(interception, energy, photosynthesis, stomatal_conductance, status))
            push!(
                status_array,
                init_variables_manual(
                    interception, energy, photosynthesis, stomatal_conductance;
                    NamedTuple{keys(status)}([j[i] for j in status_vals])...)
            )
        end
        status = status_array
    else
        status = init_variables_manual(interception, energy, photosynthesis,
            stomatal_conductance;status...)
    end

    LeafModels(interception, energy, photosynthesis, stomatal_conductance, status)
end

"""
    Base.copy(l::LeafModels)
    Base.copy(l::LeafModels, status)

Copy a [`LeafModels`](@ref), eventually with new values for the status.
"""
function Base.copy(l::T) where T <: LeafModels
    LeafModels(
        l.interception,
        l.energy,
        l.photosynthesis,
        l.stomatal_conductance,
        deepcopy(l.status)
    )
end

function Base.copy(l::T, status) where T <: LeafModels
    LeafModels(
        l.interception,
        l.energy,
        l.photosynthesis,
        l.stomatal_conductance,
        status
    )
end

"""
    Base.copy(l::AbstractArray{<:LeafModels})

Copy an array-alike of [`LeafModels`](@ref)
"""
function Base.copy(l::T) where T <: AbstractArray{<:LeafModels}
    return [copy(i) for i in l]
end

"""
    Base.copy(l::AbstractDict{N,<:LeafModels} where N)

Copy a Dict-alike of [`LeafModels`](@ref)
"""
function Base.copy(l::T) where {T <: AbstractDict{N,<:AbstractComponentModel} where N}
    return  Dict([k => v for (k, v) in l])
end

"""
    getindex(component::LeafModels,i)

Get a LeafModels component at time-step `i`.
"""
function getindex(component::LeafModels,i) where {I,E,A,Gs}
    LeafModels(
            component.interception,
            component.energy,
            component.photosynthesis,
            component.stomatal_conductance,
            component.status[i]
    )
end

# Same but with a status with only one time-step (will return the same each-time)
function getindex(component::LeafModels{I,E,A,Gs,<:MutableNamedTuples.MutableNamedTuple},i) where {I,E,A,Gs}
    component
end

"""
    DataFrame(components <: AbstractArray{<:AbstractComponentModel})
    DataFrame(components <: AbstractDict{N,<:AbstractComponentModel})

Fetch the data from a [`AbstractComponentModel`](@ref) (or an Array/Dict of) status into
a DataFrame.
"""
function DataFrame(components::T) where T <: Union{AbstractComponentModel,AbstractArray{<:AbstractComponentModel}}
    reduce(vcat, [DataFrame(i) for i in components])
end

function DataFrame(components::T) where {T <: AbstractDict{N,<:AbstractComponentModel} where N}
    df = DataFrame[]
    for (k, v) in components
        df_c = DataFrame(v)
        df_c[!, :component] .= k
        push!(df, df_c)
    end
    reduce(vcat, df)
end

# NB: could use dispatch on concrete types but would enforce specific implementation for each...
function DataFrame(components::T) where T <: AbstractComponentModel
    status = get_status(components)
    if typeof(status) == Vector{MutableNamedTuples.MutableNamedTuple}
        DataFrame([(NamedTuple(j)..., timestep = i) for (i, j) in enumerate(status)])
    else
        DataFrame([NamedTuple(status)])
    end
end

"""
    Component(interception, energy, status)
    Component(;interception = missing, energy = missing, status...)

Generic component, which is a subtype of `AbstractComponentModel` implementing a component with
an interception model and an energy balance model. It can be anything such as a trunk, a
solar panel or else.

# Arguments

- `interception <: Union{Missing,AbstractInterceptionModel}`: An interception model.
- `energy <: Union{Missing,AbstractEnergyModel}`: An energy model.
- `status <: MutableNamedTuple`: a mutable named tuple to track the status (*i.e.* the variables) of
the component. Values are set to `0.0` if not provided as VarArgs (see examples)

# Examples

```julia
# An internode in a plant:
Component(energy = Monteith())
```
"""
struct Component{I <: Union{Missing,AbstractInterceptionModel},E <: Union{Missing,AbstractEnergyModel},S <: MutableNamedTuple} <: AbstractComponentModel
    interception::I
    energy::E
    status::S
end

function Component(;interception = missing, energy = missing,status...)
    status = init_variables_manual(interception, energy;status...)
    Component(interception, energy, status)
end


text:
using StateSpace
using Distributions
using Base.Test

println("Testing Kalman Filter...")

function random_cov(n)
	A = rand(n, n)
	A = A + A'
	return A + n * eye(n)
end

srand(0)

F = diagm(ones(4)) * 0.9
V = random_cov(4)
G = [1.0 0.0 0.5 0.1;
	 0.0 0.5 1.0 1.0]
W = random_cov(2)

m = LinearGaussianSSM(F, V, G, W)
x0 = MvNormal(randn(4), diagm(ones(4) * 100.0))

x1 = predict(m, x0)
# println(mean(x1))

y1 = m.G(1) * mean(x1) + randn(2) / 10
u1 = update(m, x1, y1)
# println(mean(u1))

xx, yy = simulate(m, 100, x0)
yy[1, 50] = NaN  # throw in a missing value

# using PyPlot
# plot(xx', "k")
# plot(yy', "r")
# readline()

fs = filter(m, yy, x0)
# print(fs)

# plot(xx', "k")
# plot(mean(fs)', "r")
# readline()

y_new = fs.observations[:, end] + randn(2) / 10
update!(m, fs, y_new)

ss = smooth(m, fs)

@assert loglikelihood(fs) < loglikelihood(ss)

println("Kalman Filter passed.\n")


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export lengths, set_lengths!, angles, set_angles!
export shape, set_shape!, CellShape, volume, matrix, wrap!

__ptr(cell::UnitCell) = __ptr(cell.__handle)
__const_ptr(cell::UnitCell) = __const_ptr(cell.__handle)

"""
The possible shape for an unit cell are:

- `Chemfiles.Orthorhombic` for unit cells with the three angles are 90°.
- `Chemfiles.Triclinic` for unit cells where the three angles may not be 90°.
- `Chemfiles.Infinite` for unit cells without boundaries.
"""
@enum CellShape begin
    Infinite = lib.CHFL_CELL_INFINITE
    Orthorhombic = lib.CHFL_CELL_ORTHORHOMBIC
    Triclinic = lib.CHFL_CELL_TRICLINIC
end

"""
Create an `UnitCell` from the given `lengths` and `angles`.
"""
function UnitCell(lengths::Vector{Float64}, angles::Vector{Float64}=[90.0, 90.0, 90.0])
    @assert length(lengths) == 3
    @assert length(angles) == 3
    ptr = @__check_ptr(lib.chfl_cell(lengths, angles))
    return UnitCell(CxxPointer(ptr, is_const=false))
end

"""
Create an `UnitCell` from the given 3x3 cell `matrix`.
"""
function UnitCell(matrix::Array{Float64,2})
    @assert size(matrix) == (3, 3)
    ptr = @__check_ptr(lib.chfl_cell_from_matrix(pointer(matrix)))
    return UnitCell(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the `UnitCell` of a `frame`.
"""
function UnitCell(frame::Frame)
    ptr = @__check_ptr(lib.chfl_cell_from_frame(__const_ptr(frame)))
    cell = UnitCell(CxxPointer(ptr, is_const=false))
    copy = deepcopy(cell)
    finalize(cell)
    return copy
end

"""
Get the unit `cell` volume.
"""
function volume(cell::UnitCell)
    result = Ref{Float64}(0)
    __check(lib.chfl_cell_volume(__const_ptr(cell), result))
    return result[]
end

"""
Get the three cell lengths in angstroms.
"""
function lengths(cell::UnitCell)
    result = Float64[0, 0, 0]
    __check(lib.chfl_cell_lengths(__const_ptr(cell), result))
    return result
end

"""
Set the cell `lengths` to the given values. The lengths should be in angstroms.
"""
function set_lengths!(cell::UnitCell, lengths::Vector{Float64})
    @assert length(lengths) == 3
    __check(lib.chfl_cell_set_lengths(__ptr(cell), lengths))
    return nothing
end

"""
Get the three cell angles in degrees.
"""
function angles(cell::UnitCell)
    result = Float64[0, 0, 0]
    __check(lib.chfl_cell_angles(__const_ptr(cell), result))
    return result
end

"""
Set the cell `angles` to the given values. The angles should be in degrees.
"""
function set_angles!(cell::UnitCell, angles::Vector{Float64})
    @assert length(angles) == 3
    __check(lib.chfl_cell_set_angles(__ptr(cell), angles))
return nothing
end

"""
Get the `cell` matricial representation, *i.e.* the representation of the
three base vectors as:

        | a_x   b_x   c_x |
        |  0    b_y   c_y |
        |  0     0    c_z |
"""
function matrix(cell::UnitCell)
    result = Array{Float64}(undef, 3, 3)
    __check(lib.chfl_cell_matrix(__const_ptr(cell), pointer(result)))
    return result
end

"""
Get the `cell` shape, as a `CellShape` value.
"""
function shape(cell::UnitCell)
    result = Ref{lib.chfl_cellshape}(0)
    __check(lib.chfl_cell_shape(__const_ptr(cell), result))
    return CellShape(result[])
end

"""
Set the `cell` shape to the given `shape`.
"""
function set_shape!(cell::UnitCell, shape::CellShape)
    __check(lib.chfl_cell_set_shape(__ptr(cell), lib.chfl_cellshape(shape)))
return nothing
end

"""
Wrap a `vector` in the unit `cell`.
"""
function wrap!(cell::UnitCell, vector::Vector{Float64})
    if length(vector) != 3
        throw(ChemfilesError("Can only wrap vectors containing 3 elements in an unit cell"))
    end
    __check(lib.chfl_cell_wrap(__const_ptr(cell), vector))
    return vector
end

"""
Make a deep copy of a `cell`.
"""
function Base.deepcopy(cell::UnitCell)
    ptr = lib.chfl_cell_copy(__const_ptr(cell))
    return UnitCell(CxxPointer(ptr, is_const=false))
end


text:
module RocketNoopOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: noop()" begin

    println("Testing: operator noop()")

    run_proxyshowcheck("Noop", noop())

    source = from(1:5)

    for i in 1:1000
        source = source |> map(Int, d -> d + 1) |> noop()
    end

    run_testset([
        (
            source = source,
            values = @ts([ 1001, 1002, 1003, 1004, 1005, c ])
        ),
        (
            source = completed() |> noop(),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> noop(),
            values = @ts(e(1))
        ),
        (
            source = never() |> noop(),
            values = @ts()
        )
    ], check_timings = false)

end

end


text:
module StaticTaylorSeries

using Requires
using McCormick # TODO: Remove if package

import Base: ==, +, -, *, /, ^

import Base: iterate, size, eachindex, firstindex, lastindex,
    eltype, length, getindex, setindex!, axes, copyto!

import Base: zero, one, zeros, ones, isinf, isnan, iszero,
             convert, promote_rule, promote, show,
             real, imag, conj, adjoint,
             rem, mod, mod2pi, abs, abs2,
             sqrt, exp, log, sin, cos, tan,
             asin, acos, atan, sinh, cosh, tanh,
             power_by_squaring,
             rtoldefault, isfinite, isapprox, rad2deg, deg2rad

export STaylor1

export getcoeff, derivative, integrate, differentiate,
       evaluate, evaluate!, inverse, set_taylor1_varname,
       show_params_TaylorN, show_monomials, displayBigO, use_show_default,
       get_order, get_numvars, set_variables, get_variables,
       get_variable_names, get_variable_symbols,
       taylor_expand, update!, constant_term, linear_polynomial,
       normalize_taylor, evaluate

include("constructors.jl")
include("conversion.jl")
include("auxiliary.jl")
include("arithmetic.jl")
include("power.jl")
include("functions.jl")
include("other_functions.jl")
include("evaluate.jl")
include("printing.jl")

function __init__()
    @require IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" include("intervals.jl")
end

end # module


text:
@testset "Scheduler Tests" begin

    ancestor_1 = [
        ALifeBenchmark.OR_1(),      # 1 -> c
        ALifeBenchmark.NOP_0(),
        ALifeBenchmark.NOP_0(),
        ALifeBenchmark.JMP(),       # raises error
        ALifeBenchmark.DIVIDE()    # raises error
    ]
    ancestor_2 = [
        ALifeBenchmark.OR_1(),      # 1 -> c
        ALifeBenchmark.SH_L(),      # 2 -> c
        ALifeBenchmark.DIVIDE(),    # raises error
        ALifeBenchmark.NOP_0(),
        ALifeBenchmark.NOP_0()
    ]

    model = ALifeBenchmark.TierraModel(ancestor_1)
    ALifeBenchmark.add_organism!(model, ancestor_2)

    ancestor_1 = model.organisms[1]
    ancestor_2 = model.organisms[2]

    # test initialization

    @test ancestor_1.key == 1
    @test ancestor_2.key == 2

    @test model.reaper_queue == [1, 2]
    @test model.organism_keys == [1, 2]
    @test model.slice_index == 1

    # run program:

    # 2: OR_1
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 0
    @test ancestor_2.c == 1
    @test model.slice_index == 2
    @test model.reaper_queue == [1, 2]

    # 1: OR_1
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 1
    @test model.slice_index == 1
    @test model.reaper_queue == [1, 2]

    # 2: SH_L
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 2
    @test model.reaper_queue == [1, 2]

    # 1: NOP_0
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 1
    @test model.reaper_queue == [1, 2]

    # 2: DIVIDE -> Error
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 2
    @test model.reaper_queue == [2, 1]

    # 1: NOP_0
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 1
    @test model.reaper_queue == [2, 1]

    # 2: NOP_0
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 2
    @test model.reaper_queue == [2, 1]

    # 1: JMP -> Error
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 1
    @test model.reaper_queue == [1, 2]

    # 2: NOP_0
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 2
    @test model.reaper_queue == [1, 2]

    # 1: DIVIDE -> Error
    ALifeBenchmark.execute_slice!(model, slice_size = 1)
    @test ancestor_1.c == 1
    @test ancestor_2.c == 2
    @test model.slice_index == 1
    @test model.reaper_queue == [1, 2]
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# Before running this, please make sure to activate and instantiate the
# environment with [this `Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-breastcancer/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-breastcancer/Manifest.toml).
# For instance, copy these files to a folder 'EX-breastcancer', `cd` to it and
#
# ```julia
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Introduction
# This tutorial covers the concepts of iterative model selection on the popular ["Breast Cancer Wisconsin (Diagnostic) Data Set"](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic))
# from the UCI archives. The tutorial also covers basic data preprocessing and usage of MLJ Scientific Types.

# ## Loading the relevant packages
# For a guide to package intsllation in Julia please refer this [link](https://docs.julialang.org/en/v1/stdlib/Pkg/) taken directly from Juliav1 documentation

using UrlDownload
using DataFrames
using PrettyPrinting
using PyPlot
using MLJ

# Inititalizing a global random seed which we'll use throughout the code to maintain consistency in results

RANDOM_SEED = 42;

# ## Downloading and loading the data
# Using the package UrlDownload.jl, we can capture the data from the given link using the below commands

url = "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data";
feature_names = ["ID", "Class", "mean radius", "mean texture", "mean perimeter", "mean area", "mean smoothness", "mean compactness", "mean concavity", "mean concave points", "mean symmetry", "mean fractal dimension", "radius error", "texture error", "perimeter error", "area error", "smoothness error", "compactness error", "concavity error", "concave points error", "symmetry error", "fractal dimension error", "worst radius", "worst texture", "worst perimeter", "worst area", "worst smoothness", "worst compactness", "worst concavity", "worst concave points", "worst symmetry", "worst fractal dimension"]
data = urldownload(url, true, format = :CSV, header = feature_names);

# ## Exploring the obtained data
# ### Inspecting the class variable

figure(figsize=(8, 6))
hist(data.Class)
xlabel("Classes")
ylabel("Number of samples")

# \figalt{Distribution of target classes}{Target_class.svg}

# ### Inspecting the feature set

df = DataFrame(data)[:, 2:end];

# Printing the 1st 10 rows so as to get a visual idea about the type of data we're dealing with

pprint(first(df,10))

# For checking the statistical attributes of each inividual feature, we can use the __decsribe()__ method

pprint(describe(df))

# As we can see the feature set consists of varying features that have different ranges and quantiles. This can cause trouble for the optimization techniques and might cause convergence
# issues. We can use a feature scaling technique like __Standardizer()__ to handle this.

# But first, let's handle the [scientific types](https://alan-turing-institute.github.io/ScientificTypes.jl/dev/) of all the features. We can use the schema() method from MLJ.jl package to do this

pprint(schema(df))

# As the target variable is 'Textual' in nature, we'll have to change it to a more appropriate scientific type. Using the __coerce()__ method, let's change it to an OrderedFactor.

coerce!(df, :Class => OrderedFactor{2});

# ## Unpacking the values
# Now that our data is fully processed, we can separate the target variable 'y' from the feature set 'X' using the __unpack()__ method.

y, X = unpack(df, ==(:Class),name->true, rng = RANDOM_SEED);

# ## Standardizing the "feature set"
# Now that our feature set is separated from the target variable, we can use the __Standardizer()__ worklow to obtain to standadrize our feature set 'X'.

transformer_instance = Standardizer()
transformer_model = machine(transformer_instance, X)
fit!(transformer_model)
X = MLJ.transform(transformer_model, X);

# ## Train-test split
# After feature scaling, our data is ready to put into a Machine Learning model for classification! Using 80% of data for training, we can perform a train-test split using the __partition()__ method.

train, test = partition(eachindex(y), 0.8, shuffle=true, rng=RANDOM_SEED);

# ## Model compatibility
# Now that we have separate training and testing set, let's see the models compatible with our data!

for m in models(matching(X, y))
    println("Model name = ",m.name,", ","Prediction type = ",m.prediction_type,", ","Package name = ",m.package_name);
end

# ## Analyzing the performance of different models
# Thats a lot of models for our data! To narrow it down, lets analyze the performance of "probabilistic classifiers" from the "ScikitLearn" package.

# ### Creating various empty vectors for our analysis
# - __model_names__ captures the names of the models being iterated
# - __loss_acc captures__ the value of the model accuracy on the test set
# - __loss_ce captures__ the values of the Cross-entropy loss on the test set
# - __loss_f1__ captures the values of F1-Score on the test set

model_names=Vector{String}();
loss_acc=[];
loss_ce=[];
loss_f1=[];

# ### Collecting data for analysis

figure(figsize=(8, 6))
for m in models(matching(X, y))
    if m.prediction_type==Symbol("probabilistic") && m.package_name=="ScikitLearn" && m.name!="LogisticCVClassifier"
        #Excluding LogisticCVClassfiier as we can infer similar baseline results from the LogisticClassifier

        #Capturing the model and loading it using the @load utility
        model_name=m.name
        package_name=m.package_name
        eval(:(clf = @load $model_name pkg=$package_name verbosity=1))

        #Fitting the captured model onto the training set
        clf_machine = machine(clf(), X, y)
        fit!(clf_machine, rows=train)

        #Getting the predictions onto the test set
        y_pred = MLJ.predict(clf_machine, rows=test);

        #Plotting the ROC-AUC curve for each model being iterated
        fprs, tprs, thresholds = roc(y_pred, y[test])
        plot(fprs, tprs,label=model_name);

        #Obtaining different evaluation metrics
        ce_loss = mean(cross_entropy(y_pred,y[test]))
        acc = accuracy(mode.(y_pred), y[test])
        f1_score = f1score(mode.(y_pred), y[test])

        #Adding the different obtained values of the evaluation metrics to the respective vectors
        push!(model_names, m.name)
        append!(loss_acc, acc)
        append!(loss_ce, ce_loss)
        append!(loss_f1, f1_score)
    end
end

#Adding labels and legend to the ROC-AUC curve
xlabel("False Positive Rate")
ylabel("True Positive Rate")
legend(loc="best", fontsize="xx-small")
title("ROC curve")

# \figalt{ROC-AUC Curve}{breastcancer_auc_curve.svg}

# ### Analyzing models
# Let's collect the data in form a dataframe for a more precise analysis

model_info=DataFrame(ModelName=model_names,Accuracy=loss_acc,CrossEntropyLoss=loss_ce,F1Score=loss_f1);

# Now, let's sort the data on basis of the Cross-entropy loss

pprint(sort!(model_info,[:CrossEntropyLoss]));

# It seems like a simple LogisticClassifier works really well with this dataset!

# # Conclusion
# This article covered iterative feature selection on the Breast cancer classification dataset. In this tutorial, we only analyzed the __ScikitLearn__
# models so as to keep the flow of the content precise, but the same workflow can be applied to any compatible model in the __MLJ__ family.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using CellListMap
using StaticArrays
import Random

#
# In this example we compute the histogram of distances, expected to follow the
# function f(f) = ρ(4/3)π(r[i+1]^3 - r[i]^3) with ρ being the density of the system.
#
function distance_histogram(;N=100_000,parallel=true,x=nothing)

    # Number of particles, sides and cutoff
    sides = SVector{3,Float64}(250, 250, 250)
    cutoff = 10.
    box = Box(sides, cutoff)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    end

    # Initialize auxiliary linked lists
    cl = CellList(x, box, parallel=parallel)

    # Function to be evalulated for each pair: build distance histogram
    function build_histogram!(d2, hist)
        d = sqrt(d2)
        ibin = floor(Int, d) + 1
        hist[ibin] += 1
        return hist
    end

    # Preallocate and initialize histogram
    hist = zeros(Int, 10)

    # Run pairwise computation
    map_pairwise!(
        (x, y, i, j, d2, hist) -> build_histogram!(d2, hist),
        hist, box, cl,
        parallel=parallel
    )

    # Normalization (convert to float)
    hist = (N / (N * (N - 1) / 2)) * hist

    return hist
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of Tables.jl compatibility.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# These structures must not be defined inside a @testset. Otherwise, the test
# will fail for Julia 1.0.
struct MyColumnTable{T <: AbstractMatrix}
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    matrix::T
end

struct MyRowTable{T <: AbstractMatrix}
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    matrix::T
end

struct MyMatrixRow{T} <: Tables.AbstractRow
    row::Int
    source::MyRowTable{T}
end

@testset "Tables.jl compatibility" begin
    # A NamedTuple is compliant with Tables.jl API.
    table = (
        x = Int64(1):Int64(3),
        y = 'a':'c',
        z = ["String 1";"String 2";"String 3"]
    )

    # Thus, the following 5 calls must provide the same results.
    result_1 = pretty_table(String, table)
    result_2 = pretty_table(String, Tables.rowtable(table))
    result_3 = pretty_table(String, Tables.columntable(table))
    result_4 = pretty_table(String, Tables.columns(table))
    result_5 = pretty_table(String, Tables.rows(table))

    expected = """
┌───────┬──────┬──────────┐
│     x │    y │        z │
│ Int64 │ Char │   String │
├───────┼──────┼──────────┤
│     1 │    a │ String 1 │
│     2 │    b │ String 2 │
│     3 │    c │ String 3 │
└───────┴──────┴──────────┘
"""

    @test result_1 == result_2 == result_3 == result_4 == result_5 == expected

    # If a header is passed, then it must replace the Tables.jl schema.

    result = pretty_table(
        String,
        table;
        header = ["My col. 1", "My col. 2", "My col. 3"]
    )

    expected = """
┌───────────┬───────────┬───────────┐
│ My col. 1 │ My col. 2 │ My col. 3 │
├───────────┼───────────┼───────────┤
│         1 │         a │  String 1 │
│         2 │         b │  String 2 │
│         3 │         c │  String 3 │
└───────────┴───────────┴───────────┘
"""

    @test result == expected
end

@testset "Tables.jl without schema" begin

    expected = """
┌───┬───────┬─────┬───┐
│ a │     b │   c │ d │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    # Column table
    # --------------------------------------------------------------------------

    Tables.istable(::Type{<:MyColumnTable}) = true
    names(m::MyColumnTable) = getfield(m, :names)
    mat(m::MyColumnTable) = getfield(m, :matrix)
    lookup(m::MyColumnTable) = getfield(m, :lookup)

    Tables.columnaccess(::Type{<:MyColumnTable}) = true
    Tables.columns(m::MyColumnTable) = m
    Tables.getcolumn(m::MyColumnTable, ::Type{T}, col::Int, nm::Symbol) where {T} = mat(m)[:, col]
    Tables.getcolumn(m::MyColumnTable, nm::Symbol) = mat(m)[:, lookup(m)[nm]]
    Tables.getcolumn(m::MyColumnTable, i::Int) = mat(m)[:, i]
    Tables.columnnames(m::MyColumnTable) = names(m)

    table = MyColumnTable(
        [:a, :b, :c, :d],
        Dict(:a => 1, :b => 2, :c => 3, :d => 4),
        data
    )

    result = pretty_table(String, table)

    @test Tables.schema(table) == nothing
    @test result == expected

    # Row table
    # --------------------------------------------------------------------------

    # First, we need to create a object that complies with Tables.jl and that
    # does not have a schema. This is based on Tables.jl documentation.

    Tables.istable(::Type{<:MyRowTable}) = true
    names(m::MyRowTable) = getfield(m, :names)
    mat(m::MyRowTable) = getfield(m, :matrix)
    lookup(m::MyRowTable) = getfield(m, :lookup)

    Tables.rowaccess(::Type{<:MyRowTable}) = true
    Tables.rows(m::MyRowTable) = m
    Base.eltype(m::MyRowTable{T}) where {T} = MyMatrixRow{T}
    Base.length(m::MyRowTable) = size(mat(m), 1)
    Base.iterate(m::MyRowTable, st=1) = st > length(m) ? nothing : (MyMatrixRow(st, m), st + 1)

    Tables.getcolumn(m::MyMatrixRow, ::Type, col::Int, nm::Symbol) =
        getfield(getfield(m, :source), :matrix)[getfield(m, :row), col]
    Tables.getcolumn(m::MyMatrixRow, i::Int) =
        getfield(getfield(m, :source), :matrix)[getfield(m, :row), i]
    Tables.getcolumn(m::MyMatrixRow, nm::Symbol) = getfield(
        getfield(m, :source), :matrix)[
            getfield(m, :row),
            getfield(getfield(m, :source), :lookup)[nm]
        ]
    Tables.columnnames(m::MyMatrixRow) = names(getfield(m, :source))

    table = MyRowTable(
        [:a, :b, :c, :d],
        Dict(:a => 1, :b => 2, :c => 3, :d => 4),
        data
    )

    result = pretty_table(String, table)

    @test Tables.schema(table) == nothing
    @test result == expected
end


text:
using Test
using Random: Random, AbstractRNG, gentype, MersenneTwister, rand!, Sampler, randstring
using RandomMonad
using RandomMonad: wrap, isnormalized

const rng = MersenneTwister()

@testset "Uniform/Wrap" begin
    for W = (Uniform, wrap)
        for (dist, dest) in (Float64     => x -> 0 <= x < 1,
                             1:10        => 1:10,
                             Int         => typemin(Int):typemax(Int),
                             "asd"       => ['a', 's', 'd'],
                             Set(1:3)    => 1:3,
                             Dict(1=>2)  => [1=>2],
                             (1, 2, 3)   => 1:3)

            dist isa Tuple && VERSION < v"1.1" && continue
            w = W(dist)
            if dist isa DataType
                @test w isa RandomMonad.UniformType
            else
                @test w isa RandomMonad.UniformWrap
            end
            @test wrap(w) === w
            if dest isa Function
                @test dest(rand(w))
            else
                @test rand(w) ∈ dest
            end
        end
    end
    s = Sampler(MersenneTwister, 1:3)

    @test_throws MethodError Uniform(1.3)
    @test_throws MethodError Uniform(max)
    @test_throws MethodError Uniform(s)

    @test wrap(1.3) isa RandomMonad.Wrap
    @test wrap(s) === s
    d = Bernoulli()
    @test wrap(d) === d
    @test_throws ArgumentError rand(wrap(1.3))
    @test_throws ArgumentError rand(wrap(max))
end

@testset "support/pmf" begin
    vsp =
        (([4, 1, 4, 2], 1:4, [0.5, 0.25, 0.5, 0.25]),
         (Real[4, 1.2, 2//3], Real[2//3, 1.2, 4], [1/3, 1/3, 1/3]),
         (Integer[1, 0x3, 2], 1:3, [1/3, 1/3, 1/3]),
         ([[1], [1, 3], [4], [1]], [[1], [1, 3], [4]], [0.5, 0.25, 0.25]),
         (2:4, 2:4, [1/3, 1/3, 1/3]),
         ([identity, isless, identity], [identity, isless], [2/3, 1/3, 2/3]),
         (1:4, 1:4, [0.25, 0.25, 0.25, 0.25]),
         (3:-1:1, 1:3, [1/3, 1/3, 1/3]))

    for (v0, vs, vp) in vsp
        for v in (v0, wrap(v0), Tuple(v0), wrap(Tuple(v0)),
                  Uniform(v0), Uniform(Tuple(v0)))

            VERSION < v"1.1" &&
                (v isa Tuple || v isa Uniform && v.val isa Tuple) &&
                continue

            @test support(v) == vs
            fpmf = pmf(v)
            @test pmf(fpmf) === fpmf
            @test issubset(support(fpmf), vs)
            for (x, p) in zip(v0, vp)
                @test pmf(v, x) == p
                @test fpmf(x) == p
                @test fpmf[x] == p # not public API (yet)
            end

            T = gentype(fpmf)
            @test T == gentype(v)
            seed = rand(UInt)
            Random.seed!(rng, seed)
            x, xs = rand(rng, fpmf), rand(rng, fpmf, 3)
            @test x isa T
            @test xs isa Vector{T}

            Random.seed!(rng, seed)
            @test x == rand(rng, fpmf)
            @test xs == rand(rng, fpmf, 3)
        end
    end
    # Complex can't be sorted (just test pmf doesn't error)
    @test sum(values(pmf(randn(ComplexF64, 8)))) == 1.0
    # tuple of arrays can be sorted
    @test support(([3], [1])) == [[1], [3]]

    # normalization
    f = pmf(rand(8), normalized=false)
    @test all(==(1.0), values(f))
    @test length(keys(f)) == 8
    g = normalize!(f)
    @test g === f
    @test all(==(0.125), values(f))
    @test length(keys(f)) == 8
    g = denormalize!(f)
    @test g === f
    @test all(==(1.0), values(f))
    @test length(keys(f)) == 8

    # filter!
    f = pmf([1, 2, 3, 1])
    @test isnormalized(f)
    g = filter!(isodd, f)
    @test !isnormalized(f)
    @test f === g
    @test keys(f)  == [1, 3]
    @test collect(values(f)) == [1/2, 1/4]
    normalize!(f)
    @test isnormalized(f)
    @test keys(f)  == [1, 3]
    @test collect(values(f)) == [2/3, 1/3]
    denormalize!(f)
    @test collect(values(f)) == [1/2, 1/4]

    f = pmf([1, 2, 3, 1])
    denormalize!(f)
    @test !isnormalized(f)
    g = filter!(isodd, f)
    @test !isnormalized(f)
    @test f === g
    @test keys(f)  == [1, 3]
    @test collect(values(f)) == [2.0, 1.0]
    normalize!(f)
    @test isnormalized(f)
    @test keys(f)  == [1, 3]
    @test collect(values(f)) == [2/3, 1/3]
    denormalize!(f)
    @test collect(values(f)) == [2.0, 1.0]
end

@testset "Bernoulli" begin
    @test rand(Bernoulli()) ∈ (0, 1)
    @test rand(Bernoulli()) isa Bool
    for T = (Int, UInt, Bool)
        @test rand(Bernoulli(T)) ∈ (0, 1)
    end
    @test rand(Bernoulli(1)) == 1
    @test rand(Bernoulli(0)) == 0
    # TODO: do the math to estimate proba of failure:
    @test 620 < count(rand(Bernoulli(Bool, 0.7), 1000)) < 780
    for T = (Bool, Int, Float64, ComplexF64)
        r = rand(Bernoulli(T))
        @test r isa T
        @test r ∈ (0, 1)
        r = rand(Bernoulli(T, 1))
        @test r == 1
    end

    b = Bernoulli(0.3)
    @test support(b) == false:true
    for bpmf = (x -> pmf(b, x), pmf(b))
        @test bpmf(-1) == pmf(b, 2) == 0
        @test bpmf(0) == 0.7
        @test bpmf(1) == 0.3
    end
end

@testset "Binomial" begin
    b = Binomial(1000, 0.7)
    @test rand(b) isa Int
    @test eltype(b) == Int
    @test 620 < rand(b) < 780
    b = Binomial(1000)
    @test b.p == 0.5
    @test 420 < rand(b) < 580
    b = Binomial()
    @test b.n == 1
    @test rand(b) ∈ 0:1

    for p in rand(10)
        for n in rand(Shuffle(1:30), 5)
            b = Binomial(n, p)
            @test support(b) == 0:n # unlikely that p ∈ (0.0, 1.0)
            @test sum(pmf(b, x) for x in support(b)) ≈ 1.0
            @test sum(values(pmf(b))) ≈ 1.0
        end
    end
end

@testset "Categorical" begin
    n = rand(1:9)
    @test rand(Categorical(n)) ∈ 1:9
    @test all(∈(1:9), rand(Categorical(n), 10))
    @test rand(Categorical(n)) isa Int
    c = Categorical(Float64(n))
    @test rand(c) isa Int
    @test rand(c) ∈ 1:9
    c = Categorical([1, 7, 2])
    # cf. Bernoulli tests
    @test 620 < count(==(2), rand(c, 1000)) < 780
    @test rand(c) isa Int
    @test rand(Categorical(Float64.((1, 2, 3, 4)))) isa Int
    @test rand(Categorical(3, 1.0:3.0)) isa Float64

    @test_throws ArgumentError Categorical(())
    @test_throws ArgumentError Categorical([])
    @test_throws ArgumentError Categorical(x for x in 1:0)
end

@testset "Multinomial" begin
    for m in (Multinomial(10, Categorical(3)),
              Multinomial(10, 3))
        @test rand(m) isa Vector{Int}
        for a in rand(m, 3)
            @test sum(a) == 10
        end
    end

    for m in (Multinomial(10, Categorical([1, 0, 1])),
              Multinomial(10, [1, 0, 1]))
        a = rand(m)
        @test a[2] == 0
        @test a[1] + a[3] == 10

        # rand! & Pack
        b = rand!(a, m, Val(1))
        @test a === b
        @test a[2] == 0
        @test a[1] + a[3] == 10

        a = rand(Pack(m, 5))
        @test size(a) == (3, 5)
        for j=1:5
            @test a[1, j] + a[3, j] == 10
            @test a[2, j] == 0
        end
    end

    # Pack with equal weigths
    for m in (Multinomial(10, Categorical(3)),
              Multinomial(10, 3))
        a = rand(Pack(m, 20))
        @test size(a) == (3, 20)
    end
end

@testset "MixtureModel" begin
    m = MixtureModel([1, 7, 2], [100:300, Normal(), Uniform(500:1000)], )
    x = rand(m)
    testtype(x) = x isa Float64 || x isa Int
    @test testtype(x)
    xs = rand(m, 1000)
    @test all(testtype, xs)
    # cf. Bernoulli tests
    @test 620 < count(x -> x isa Float64, xs) < 780

    m = MixtureModel((1,2), [Normal(), CloseOpen()])
    @test eltype(m.cat.components) == Distribution{Float64} # not crucial

    m = MixtureModel((0.5, 0.5), [Normal(), CloseOpen()])
    @test m.cat.cdf == [0.5, 1.0]
    # test that the second argument doesn't need to have length defined
    if VERSION >= v"1.4" # for takewhile
        m = MixtureModel((1/3, 1/3, 1/3),
                         Iterators.takewhile(_->true,
                                             (Normal(), Exponential(), CloseOpen())))
        @test m.cat.cdf == [1/3, 2/3, 1.0]
    end

    @test_throws ArgumentError MixtureModel([1, 2, 3], [1:3, Normal()])
end

@testset "Normal" begin
    @test rand(Normal()) isa Float64
    @test rand(Normal(0.0, 1.0)) isa Float64
    @test rand(Normal(0, 1)) isa Float64
    @test rand(Normal(0, 1.0)) isa Float64
    @test rand(Normal(Float32)) isa Float32
    @test rand(Normal(ComplexF64)) isa ComplexF64
end

@testset "Exponential" begin
    @test rand(Exponential()) isa Float64
    @test rand(Exponential(1.0)) isa Float64
    @test rand(Exponential(1)) isa Float64
    @test rand(Exponential(Float32)) isa Float32
end

@testset "Poisson" begin
    p = Poisson()
    @test p.λ == 1
    @test rand(p) isa Int
    @test rand(Poisson(2)) isa Int
    @test rand(Poisson(3.0)) isa Int
    @test rand(Poisson(0)) == 0
end

@testset "rand(::AbstractFloat)" begin
    # check that overridden methods still work
    m = MersenneTwister()
    for F in (Float16, Float32, Float64, BigFloat)
        @test rand(F) isa F
        sp = Random.Sampler(MersenneTwister, RandomMonad.CloseOpen01(F))
        @test rand(m, sp) isa F
        @test 0 <= rand(m, sp) < 1
        for (CO, (l, r)) = (CloseOpen  => (<=, <),
                            CloseClose => (<=, <=),
                            OpenOpen   => (<,  <),
                            OpenClose  => (<,  <=))
            f = rand(CO(F))
            @test f isa F
            @test l(0, f) && r(f, 1)
        end
        F ∈ (Float64, BigFloat) || continue # only types implemented in Random
        sp = Random.Sampler(MersenneTwister, RandomMonad.CloseOpen12(F))
        @test rand(m, sp) isa F
        @test 1 <= rand(m, sp) < 2
    end
    @test CloseOpen(1,   2)          === CloseOpen(1.0, 2.0)
    @test CloseOpen(1.0, 2)          === CloseOpen(1.0, 2.0)
    @test CloseOpen(1,   2.0)        === CloseOpen(1.0, 2.0)
    @test CloseOpen(1.0, Float32(2)) === CloseOpen(1.0, 2.0)
    @test CloseOpen(big(1), 2) isa CloseOpen{BigFloat}

    for CO in (CloseOpen, CloseClose, OpenOpen, OpenClose)
        @test_throws ArgumentError CO(1, 1)
        @test_throws ArgumentError CO(2, 1)

        @test CO(Float16(1), 2) isa CO{Float16}
        @test CO(1, Float32(2)) isa CO{Float32}
    end
end

## adapters

@testset "Pure" begin
    @test rand(Pure(1)) === 1
    @test rand(Pure(1:3)) === 1:3
    @test eltype(Pure(1:3)) == UnitRange{Int}
    s = Set([1, 2, 3])
    @test rand(Pure(s)) === s
    a = rand(Pure(s), 9)
    @test eltype(a) == Set{Int}
    @test all(x -> x === s, a)
end

@testset "algebra" begin
    d = Uniform(Float64) + Uniform(1:3)
    @test eltype(d) == Float64
    @test all(x -> 1.0 <= x < 4, rand(d, 100))
    for op = (+, -, *, ^, /)
        # Keep below for y^x needing x >= 0
        for x = (Keep(x -> x >= 0, Uniform(Int32)), Int32(3))
            for d = (op(x, Uniform(Int32(1):Int32(3))),
                     op(Uniform(Int32(1):Int32(3)), x))
                if op == /
                    @test eltype(d) == Float64
                    @test rand(d) isa Float64
                else
                    @test eltype(d) == Int32
                    @test rand(d) isa Int32
                end
            end
        end
    end
    # test reset!
    @test all(==(1:9), sort!.(rand(Fill(Unique(1:9) + Uniform(0:0), 9), 2)))

    # getindex
    g = Pure('a':'z')[Categorical(3)]
    @test rand(g) ∈ 'a':'c'
    g = Pure(Dict(1=>3, 2=>4))[Uniform(1:2)]
    @test rand(g) ∈ 3:4
end

@testset "Bind" begin
    b = Bind(n -> Fill(Bool, n), 1:5)
    @test eltype(b) == Vector{Bool}
    @test rand(b) isa Vector{Bool}
    for v in rand(b, 4)
        @test v isa Vector{Bool}
        @test length(v) in 1:5
    end

    b = Bind(Zip([Bool, Normal(), 1:9], Shuffle(1:4))) do (t, n)
               Fill(t, n)
           end
    vs = rand(b, 4)
    @test allunique(length.(vs))
    @test all(v -> length(v) ∈ 1:4, vs)
    @test all(v -> eltype(v) ∈ [Bool, Float64, Int], vs)

    # constructor accepts samplers
    @test rand(Bind(Pure, Sampler(MersenneTwister, 1:3))) isa Int
    # Bind doesn't store DataType, but a UniformType
    @test fieldtype(typeof(Bind(Pure, Float64)), 2) == RandomMonad.UniformType{Float64}

    # reset!
    v = rand(Fill(Bind(Pure, Iterate(1:99)), 3), 2)
    @test all(==(1:3), v)
    sp = Sampler(MersenneTwister, Bind(Pure, Iterate(1:99)), Val(1))
    @test rand(sp, 3) == 1:3
    @test rand(sp, 3) == 1:3
end

@testset "Join" begin
    j = Join(Uniform([Uniform(Float64) + 10, Uniform(Float64)]))
    @test eltype(j) == Float64
    @test all(rand(j, 100)) do x
        0 <= x < 1 || 10 <= x < 11
    end
end

@testset "Thunk" begin
    t = let i::Int = 0
            Thunk() do
                i += 1
                Uniform(1:i)
            end
        end
    @test eltype(t) == Int # because of the i::Int in let
    @test rand(t) == 1
    @test rand(t) ∈ 1:2
    v = rand(t, 3)
    for i in 1:3
        @test v[i] ∈ 1:(i+2)
    end
end

@testset "Lift" begin
    d = Lift(x -> 2x, 1:3)
    @test rand(d) ∈ 2:2:6
    @test all(x -> x ∈ 2:2:6, rand(d, 100))
    @test eltype(d) == Int
    @test rand(d) isa Int

    d = Lift{Float64}(x -> x > 0, Normal())
    @test rand(d) isa Float64
    @test all(x -> x ∈ (0.0, 1.0), rand(d, 100))

    d = Lift(+, Float64, [10, 20])
    @test all(rand(d, 30)) do x
        10 <= x < 11 || 20 <= x < 21
    end

    # pmf
    d = Lift(x -> x^2, -2:2)
    f = pmf(d)
    @test keys(f) == [0, 1, 4]
    @test collect(values(f)) == [0.2, 0.4, 0.4]
    f = pmf(Lift(sum, Zip(1:3, (-1, 1))))
    @test f[0] == f[1] == f[3] == f[4] == 1/6
    @test f[2] == 1/3
end

@testset "Keep" begin
    d = Keep(x -> x > 0, Normal())
    @test all(x -> x > 0, rand(d, 1000))
    d = Keep(x -> x < 4, Categorical([1, 7, 2, 10]))
    # cf. Bernoulli tests
    @test 620 < count(==(2), rand(d, 1000)) < 780
    d = Keep(x -> x < 9, 1:20)
    @test all(x -> 1 <= x <= 9, rand(d, 1000))

    # reset!
    s = Random.Sampler(MersenneTwister, Keep(x -> x < 5, Shuffle(1:9)))
    for _ = 1:3
        a = rand(s, 4)
        @test allunique(a)
        @test all(x -> x < 5, a)
    end

    # pmf
    f = pmf(Keep(iseven, 1:9))
    @test f == Dict(2=>1/4, 4=>1/4, 6=>1/4, 8=>1/4)
    f = pmf(Keep(isodd, 1:9))
    @test f == Dict((1:2:9) .=> prevfloat(0.2)) # prevfloat b/c of rounding
end

@testset "Map" begin
    d = Map(x -> 2x, Fill(1:3, 3))
    @test all(∈(2:2:6), rand(d))
    @test all(x -> all(∈(2:2:6), x), rand(d, 10))
    @test eltype(d) == Vector{Int}
    @test rand(d) isa Vector{Int}

    d = Map{Vector{Float64}}(x -> x > 0, Fill(Normal(), 10))
    @test rand(d) isa Vector{Float64}
    @test all(x -> all(∈((0.0, 1.0)), x), rand(d, 5))

    d = Map(+, Fill(Float64, 2), Fill([10, 20], 2))
    @test all(rand(d)) do x
        10 <= x < 11 || 20 <= x < 21
    end
end

@testset "Filter" begin
    f = Filter(x -> x != 0, Fill(-1:1, 10))
    @test eltype(f) == Vector{Int}
    @test all(∈([-1, 1]), rand(f))
    f = Filter(x -> x > 0, Fill(Normal(), 10))
    @test eltype(f) == Vector{Float64}
    @test all(x -> x > 0, rand(f))

    @test_throws MethodError rand(Filter(x -> x != 0, Zip(-1:1, -1:1)), 40)

    if VERSION >= v"1.4"
        z1 = rand(Filter{Tuple}(x -> x != 0, Zip(-1:1, -1:1)), 3)
        @test z1 isa Vector{Tuple}
        z2 = rand(Lift(filter, Pure(x -> x != 0), Zip(-1:1, -1:1)), 3)
        @test z2 isa Vector{Tuple{Vararg{Int64}}}
        for z ∈ (z1, z2), t ∈ z
            @test all(∈((-1, 1)), t)
        end
    end

    # test reset!
    if VERSION >= v"1.4"
        a = rand(Fill(Filter(_->true, Zip(Unique(1:4))), 4), 2)
        @test a isa Vector{Vector{Tuple{Int}}}
    end
end

@testset "Reduce" begin
    r = Reduce(+, Fill(1:3, 2))
    @test rand(r) ∈ 2:6
    @test eltype(r) == Int
    if VERSION >= v"1.1"
        @test all(x -> x==9, rand(Reduce(+, Zip((8,), (1,))), 10))
    end
    f = Fill(Reduce(+, Unique(1:4)), 2)
    @test rand(f) isa Vector{Int}
    # check that reset! works properly (otherwise, rand(f, 3) would never terminate)
    v = rand(f, 3)
    @test v isa Vector{Vector{Int}}
    @test all(allunique, v)
end

@testset "Counts" begin
    c = Counts(Fill(0x1:0x3, 100))
    @test eltype(c) == Dict{UInt8,Int}
    d = rand(c)
    @test d isa Dict{UInt8,Int}
    @test keys(d) ⊆ 1:3
    @test 100 <= sum(values(d)) <= 300
    # check reset! works
    vs = rand(Fill(Counts(Unique(1:3)), 3), 3)
    @test vs isa Vector{Vector{Dict{Int,Int}}}
    for v in vs
        ks = []
        for d in v
            append!(ks, keys(d))
            @test collect(values(d)) == 1:1
        end
        @test sort!(ks) == 1:3
    end
end

@testset "Unique" begin
    u = Unique(1:3)
    @test eltype(u) == Int
    @test rand(u) ∈ 1:3
    a = rand(u, 3)::Vector{Int}
    @test allunique(a)

    z = Fill(Zip(u, u), 3)
    if VERSION >= v"1.1"
        for i=1:3
            rz = rand(z)
            for a = (rand(Fill(u, 3)), first.(rz), last.(rz),
                     rand(Fill(Lift(identity, u), 3)),
                     rand(Fill(Lift(+, (0,), u), 3)),
                     rand(Keep(x->true, u)))
                @test all(in(1:3), a)
                @test allunique(a)
            end
        end
    end

    u = Unique(Bool)
    @test allunique(rand(u, 2))
end

@testset "Iterated" begin
    # simulates Geometric distribution
    # tests that one can use samplers, even if quite useless here
    geom(p) = Iterated{Int}(
        begin
            function it(rng)
                sp = Sampler(rng, Bernoulli(p))
                sp, (sp, 0)
             end

            it(v, sp_i) = v ?
                (nothing, sp_i[2]) :
                (sp_i[1], (sp_i[1], sp_i[2] + 1))
        end)
    f = pmf(rand(geom(0.8), 100), normalized=false)
    @test f(0) > 55 # > about 60 in maybe 1 in a million cases

    # result of games in three parts where one must win at least twice
    # your chance of success increases if you won last game
    halting = Iterated{Vector{Bool}}(
        function (v, res=Bool[])
            if v isa AbstractRNG
                return Bernoulli(Bool), res
            else
                push!(res, v)
                (sum(res) == 2 || length(res) == 3) && return nothing, res
                Bernoulli(Bool, v ? 2/3 : 1/3), res
            end
        end)
    @test issubset(rand(halting, 5),
                   [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1],
                    [1, 0, 0], [1, 0, 1], [1, 1]])
    f = pmf(halting)
    @test collect(values(f)) ≈ [4, 2, 1, 2, 2, 1, 6] / 18
    denormalize!(f)
    @test collect(values(f)) ≈ [4, 2, 1, 2, 2, 1, 6]
end

@testset "Iterate" begin
    it = Iterate(1:0)
    @test_throws ArgumentError rand(it)
    @test_throws ArgumentError rand(it, 2)
    a = rand(it, 0)
    @test a isa Vector{Int}
    @test isempty(a)

    it = Iterate(1:3)
    @test rand(it) == 1
    @test rand(it) == 1
    @test rand(it, 3) == 1:3
    @test_throws ArgumentError rand(it, 4)

    a = rand(Fill(it, 3), 2)
    for v in a
        @test v == 1:3
    end

    it = Iterate(Iterators.countfrom(9))
    @test rand(it) == 9
    @test rand(it) == 9
    @test rand(Fill(it, 9)) == 9:17
end

@testset "SubSeq" begin
    s = SubSeq(1:100, .1)
    @test eltype(s) == Vector{Int}
    a = rand(s)
    b = Int[]
    c = rand!(b, s, Val(1))
    @test b === c
    for v = (a, b)
        @test v isa Vector{Int}
        @test issorted(v)
        @test allunique(v)
        @test all(∈(1:100), v)
    end
end

@testset "SubIter" begin
    s = SubIter(Iterators.countfrom(1), .1)
    @test eltype(s) == Int
    for n = (Val(1), Val(Inf))
        sp = Sampler(MersenneTwister, s, n)
        v = [rand(sp) for _=1:10]
        @test issorted(v)
        @test allunique(v)
        @test v isa Vector{Int}
        @test v[end] > 10 # very unlikely to fail
        RandomMonad.reset!(sp)
        @test rand(sp) < v[end] # very unlikely to fail
    end
end

@testset "$ShuffleAlgo" for ShuffleAlgo = (FisherYates, SelfAvoid, Shuffle)
    u = ShuffleAlgo(1:3)
    @test eltype(u) == Int
    @test rand(u) ∈ 1:3
    a = rand(u, 3)::Vector{Int}
    @test allunique(a)

    f = Fill(u, 3)
    z = Fill(Zip(u, u), 3)
    for i=1:3
        rz = rand(z)
        for a = (rand(f), first.(rz), last.(rz))
            @test all(in(1:3), a)
            @test allunique(a)
        end
    end

    u = ShuffleAlgo('a':'z')
    @test rand(u) isa Char
    @test allunique(rand(u, 26))

    # test reset!
    a = rand(Fill(Fill(ShuffleAlgo(1:9), 5), 2))
    @test a isa Vector{Vector{Int}}
    @test !allunique(vcat(a...))

    # test reset!(sp, 1) works
    @test all(in(1:4), vcat(rand(Fill(ShuffleAlgo(1:4), 1), 10)...))
    # test reset!(sp, 2) works
    @test all(rand(Fill(ShuffleAlgo(1:4), 2), 100)) do x
        x[1] != x[2]
    end

    # test Val(1)-sampler works
    s = Random.Sampler(MersenneTwister, ShuffleAlgo(1:9), Val(1))
    @test allunique(rand(s, 9))

    # old bug with FisherYates
    @test counts(Lift(Fill(Fill(ShuffleAlgo(1:9), 1), 4)) do xs
                      for i = 1:length(xs[1])
                          all(x -> x[i] == xs[1][i], xs) && return true
                      end
                      false
                 end,
                 1000)[false] > 990 # would fail probably less than once in 30000

    # compatibility with rand!
    s = Random.Sampler(MersenneTwister, ShuffleAlgo(1:4))
    # when rand(s, 4) is called in a loop, we test that `reset!(s, 4)`
    # is called each time; otherwise, `s` would get exhausted after
    # the first call and an exception would be thrown
    for _ = 1:3
        a = rand(s, 4)
        @test all(in(1:4), a)
        @test allunique(a)
    end

    # multivariate
    if ShuffleAlgo == Shuffle
        a = Char[]

        s = Shuffle(['a'], 0)
        @test rand(s) == Char[]
        @test rand!(a, s, Val(1)) == Char[]

        s = Shuffle(['a'], 1)
        @test rand(s) == ['a']
        @test rand!(a, s, Val(1)) == ['a']

        s = Shuffle(['a', 'b'], 2)
        @test sort!(rand(s)) == ['a', 'b']
        @test sort!(rand!(a, s, Val(1))) == ['a', 'b']

        s = Shuffle(['a', 'b', 'c'], 3)
        @test sort!(rand(s)) == ['a', 'b', 'c']
        @test sort!(rand!(a, s, Val(1))) == ['a', 'b', 'c']

        s = Shuffle(collect(randstring(100)), 50) # FisherYates
        v = rand(s)
        @test v isa Vector{Char}
        @test length(v) == 50
        @test rand!(a, s, Val(1)) === a
        @test length(a) == 50

        s = Shuffle(collect(randstring(1000)), 2) # SelfAvoid
        v = rand(s)
        @test v isa Vector{Char}
        @test length(v) == 2
        @test rand!(a, s, Val(1)) === a
        @test length(a) == 2
    end
end

## containers

@testset "Zip" begin
    z = Zip(Normal(), 1:3)
    @test eltype(z) == Tuple{Float64,Int}
    @test rand(z) isa Tuple{Float64,Int}
    @test all(x -> x ∈ 1:3, last.(rand(z, 100)))

    z = Zip()
    @test rand(z) == ()
    z = Zip(1:3)
    @test rand(z) isa Tuple{Int}
    @test rand(z)[1] in 1:3
    z = Zip(1:2, Normal(), Int8)
    @test rand(z) isa Tuple{Int,Float64,Int8}
    @test all(x -> x isa Tuple{Int,Float64,Int8}, rand(z, 9))

    # Zip(::Sampler...)
    s = Sampler(MersenneTwister, 1:3)
    z = Zip(s, s)
    @test rand(z) isa Tuple{Int,Int}
    @test all(in(1:3), rand(z))
    for t in rand(z, 3)
        @test all(in(1:3), rand(z))
    end

    # recursive rand!
    z1 = Zip(Fill(1:3, 2), Fill(1:4, 3))
    z2 = Zip(Fill(4:6, 2), Fill(5:8, 3))
    t = rand(z1)
    a, b = t[1][1], t[2][1]
    rand!(t, z2, Val(2))
    @test t[1][1] != a
    @test t[2][1] != b
    @test all(∈(5:8), t[2])
    rand!([t], Fill(z1, 1), Val(3))
    @test all(∈(1:4), t[2])
    @test_throws ArgumentError rand!(t, Zip(1:2, 1:3), Val(1)) # can't mutate tuple
    @test_throws ArgumentError rand!((t..., t...), z2, Val(2)) # not same size

    # check that sub-samplers are initialized in Sampler(rng, ::Zip, Val(1))
    for n = (1, Inf)
        sp = Sampler(MersenneTwister, Zip(Iterate(1:9)), Val(n))
        for i = 1:9
            @test rand(sp) == (i,)
        end
        sp = Sampler(MersenneTwister, Zip(Shuffle(1:9)), Val(n))
        v = [rand(sp) for _=1:9]
        @test allunique(v)
    end

    # support & pmf
    z = Zip(1:2, [1, 3, 1])
    @test vec(collect(support(z))) ==
        [(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)]
    for zpmf in (x -> pmf(z, x), pmf(z))
        @test zpmf((1, 1)) == 1/3
        @test zpmf((2, 1)) == 1/3
        @test zpmf((1, 2)) == 0
        @test zpmf((2, 2)) == 0
        @test zpmf((1, 3)) == 1/6
        @test zpmf((2, 3)) == 1/6
        @test zpmf((1,))   == 0.0
        @test zpmf((1, 2, 3)) == 0.0
    end
end

@testset "Fill" begin
    rng = MersenneTwister()
    for f in (Fill(1:9, (2, 3)),
              Fill(Uniform(1:9), 2, 3))
        @test f isa Distribution{Array{Int,2}}
        a = rand(f)
        b = similar(a)
        rand!(rng, b, f, Val(1))
        c = similar(a)
        rand!(c, f, Val(1))
        for A = (a, b, c)
            @test A isa Array{Int,2}
            @test size(A) == (2, 3)
            @test all(x -> x ∈ 1:9, A)
        end
        @test_throws ArgumentError rand!(similar(a, 3, 2), f, Val(1))
        @test_throws ArgumentError rand!(similar(a, 6), f, Val(1))
    end
    f = Fill(Float64)
    @test rand(f) isa Array{Float64,0}
    f = Fill(Float64, 0x2, Int128(3))
    a = rand(f)
    @test a isa Array{Float64,2}
    @test size(a) == (2, 3)
    f = Fill(1:9, 9)
    a = rand!(Int8[], f, Val(1))
    @test a isa Vector{Int8}
    @test length(a) == 9
    @test all(in(1:9), a)

    # Fill(::Sampler, ...)
    s = Sampler(MersenneTwister, 1:3)
    f = Fill(s, 2)
    @test all(in(1:3), rand(f))

    # recursive rand!
    inner = Fill(1:3, 3)
    a = rand(Fill(inner, 4))
    b = copy(a)
    rand!(a, inner)
    @test a[1] !== b[1]
    if VERSION >= v"1.1" # copy!
        copy!(b, a)
        rand!(a, Fill(inner, 4), Val(1))
        @test a[1] !== b[1]
        copy!(b, a)
        rand!(a, Fill(inner, 4), Val(2))
        @test a[1] === b[1]
    end

    # pmf
    f = pmf(Fill(1:2, 3))
    @test keys(f) ==
        collect.(sort(vec(collect(Iterators.product(1:2, 1:2, 1:2)))))
    @test all(==(1/8), values(f))
    @test isnormalized(f)
    denormalize!(f)
    @test all(==(1.0), values(f))
    f = pmf(Fill([1, 1, 2], 2), normalized=false)
    @test Dict(f) ==
        Dict([2, 2] => 1.0,[1, 1] => 4.0,[1, 2] => 2.0,[2, 1] => 2.0)
end

@testset "variate_size" begin
    @test variate_size(Normal()) == ()
    @test variate_size(Fill(Normal(), 2, 3)) == (2, 3)
    @test_throws MethodError variate_size(SubSeq(1:9, 0.3))
end

@testset "Pack" begin
    d = Fill(1:3)
    p = Pack(d, 2)
    @test eltype(p) == Vector{Int}
    v = rand(p)
    w = rand!(copy(v), p, Val(1))
    for u = (v, w)
        @test u isa Vector{Int}
        @test size(u) == (2,)
        @test all(∈(1:3), u)
    end

    d = Shuffle(1:9, 9)
    p = Pack(d)
    @test eltype(p) == Vector{Int}
    @test size(rand(p)) == (9,)

    d = Shuffle([[1], [1, 2]])
    @test_throws ArgumentError Pack(d)

    # Fill-like behavior for univariate distributions
    d = Shuffle(1:6)
    p = Pack(d, 2, 3)
    v = rand(p)
    @test allunique(v)
    @test all(∈(1:6), v)
    @test size(v) == (2, 3)

    d = Fill(1:9, 9)
    p = Pack(d, 0)
    @test size(rand(p)) == (9, 0)

    d = Fill(1:9) # 0-d
    p = Pack(d)
    @test rand(p) isa Array{Int,0}
    p = Pack(d, 3)
    @test rand(p) isa Vector{Int}

    d = Fill(Shuffle(1:9), 3, 3)
    for p = (Pack(d, 2, 0x4),
             Pack(d, (2, 4)))
        @test eltype(p) == Array{Int,4}
        v = rand(p)
        @test v isa Array{Int,4}
        @test size(v) == (3,3,2,4)
        for idx in CartesianIndices((2, 4))
            sub = view(v, CartesianIndices((3, 3)), idx)
            @test all(∈(1:9), sub)
            @test allunique(sub)
        end
    end
end

@testset "discretize!" begin # TODO: more tests
    f = pmf(discretize!(randn(1000), 10))
    @test length(keys(f)) ∈ (8, 9, 10) # 8 in about 1/1000 cases
end


text:
function install(m::Module; kw...)
    options = read_options(m; kw...)
    return install(m, options)
end

function print_install_help(io::IO)
    println(io, "Comonicon - Installation CLI.")
    println(io)
    println(io, "Install the CLI script to `.julia/bin` if not specified with subcommands.")
    println(io)
    printstyled(io, "USAGE\n\n"; bold = true)
    printstyled(io, " "^4, "julia --project deps/build.jl [command]\n\n"; color = :cyan)
    printstyled(io, "COMMAND\n\n"; bold = true)

    printstyled(io, " "^4, "app"; color = :light_blue, bold = true)
    printstyled(io, " [tarball]"; color = :blue)
    println(io, " "^15, "build the application, optionally make a tarball.\n")

    printstyled(io, " "^4, "sysimg"; color = :light_blue, bold = true)
    printstyled(io, " [tarball]"; color = :blue)
    println(io, " "^12, "build the system image, optionally make a tarball.\n")

    printstyled(io, " "^4, "tarball"; color = :light_blue, bold = true)
    println(io, " "^21, "build application and system image then make tarballs")
    println(io, " "^32, "for them.\n")

    printstyled(io, "EXAMPLE\n\n"; bold = true)
    printstyled(io, " "^4, "julia --project deps/build.jl sysimg\n\n"; color = :cyan)
    println(
        io,
        " "^4,
        "build the system image in the path defined by Comonicon.toml or in deps by default.\n\n",
    )
    printstyled(io, " "^4, "julia --project deps/build.jl sysimg tarball\n\n"; color = :cyan)
    println(io, " "^4, "build the system image then make a tarball on this system image.\n\n")
    printstyled(io, " "^4, "julia --project deps/build.jl app tarball\n\n"; color = :cyan)
    println(
        io,
        " "^4,
        "build the application based on Comonicon.toml and make a tarball from it.\n\n",
    )
end

function install(m::Module, options::ComoniconOptions.Comonicon)
    if isempty(ARGS)
        if options.install.quiet
            logger = NullLogger()
        else
            logger = ConsoleLogger()
        end

        with_logger(logger) do
            install_script(m, options)
        end
        return
    elseif "-h" in ARGS || "--help" in ARGS || "help" in ARGS
        return print_install_help(stdout)
    elseif first(ARGS) == "sysimg" && !isnothing(options.sysimg)
        if length(ARGS) == 1
            return build_sysimg(m, options)
        elseif length(ARGS) == 2 && ARGS[2] == "tarball"
            return build_tarball_sysimg(m, options)
        end
    elseif first(ARGS) == "app" && !isnothing(options.application)
        if length(ARGS) == 1
            return build_application(m, options)
        elseif length(ARGS) == 2 && ARGS[2] == "tarball"
            return build_tarball_app(m, options)
        end
    elseif first(ARGS) == "tarball" && (!isnothing(options.sysimg) || !isnothing(options.application))
        if length(ARGS) == 1
            return build_tarball(m, options)
        end
    end
end

"install a script as the CLI"
function install_script(m::Module, options::ComoniconOptions.Comonicon; project=pkgdir(m))
    bin = expanduser(joinpath(options.install.path, "bin"))
    shadow = joinpath(bin, options.name * ".jl")

    if isnothing(options.sysimg)
        sysimg_path = nothing
    else
        download_sysimg(m, options)
        sysimg_path = pkgdir(m, options.sysimg.path, "lib", sysimg(options.name))
    end

    shell_script = cmd_script(
        m,
        shadow;
        project,
        sysimg = sysimg_path,
        compile = options.install.compile,
        optimize = options.install.optimize,
    )
    file = joinpath(bin, options.name)

    # start writing
    if !ispath(bin)
        @info "cannot find Julia bin folder creating $bin"
        mkpath(bin)
    end

    # generate contents
    @info "generating $shadow"
    open(shadow, "w+") do f
        println(
            f,
            "#= generated by Comonicon for $(options.name) =# using $m; exit($m.command_main())",
        )
    end

    @info "generating $file"
    open(file, "w+") do f
        println(f, shell_script)
    end

    if options.install.completion
        install_completion(m, options)
    end

    chmod(file, 0o777)
    return
end

function install_completion(m::Module, options::ComoniconOptions.Comonicon, shell=detect_shell())
    completions_dir = expanduser(joinpath(options.install.path, "completions"))
    isnothing(shell) && return

    @info "generating auto-completion script for $shell"
    script = completion_script(shell, m)
    script === nothing && return

    if !ispath(completions_dir)
        mkpath(completions_dir)
    end

    completion_file = joinpath(completions_dir, "_" * options.name)
    @info "writing to $completion_file"
    write(completion_file, script)
    return
end


text:
using SymbolicControlSystems
using Test
using ControlSystems
using SymbolicControlSystems: s, z

macro test_both(G,sym)
    quote
        if isdiscrete($(esc(G)))
            @test Sym($(esc(G))) == $(esc(sym))
            @test $(esc(G)) == tf($(esc(sym)), $(esc(G)).Ts)
        else
            @test Sym($(esc(G))) == $(esc(sym))
            @test $(esc(G)) == tf($(esc(sym)))
            @test tf(ss($(esc(G)))) == tf($(esc(sym)))
        end
    end
end

@testset "SymbolicControlSystems.jl" begin
    @testset "To symbolic" begin
        @info "Testing To symbolic"

        pol = 2s^2+3s+4
        @test SymbolicControlSystems.expand_coeffs(pol, s) == [2,3,4]

        pol = 2s^3+3s^2+4*s
        @test SymbolicControlSystems.expand_coeffs(pol, s) == [2,3,4,0]

        pol = 2s^3+4*s
        @test SymbolicControlSystems.expand_coeffs(pol, s) == [2,0,4,0]

        pol = 2z^2+3z+4
        @test SymbolicControlSystems.expand_coeffs(pol, z) == [2,3,4]

        pol = 2z^3+3z^2+4*z
        @test SymbolicControlSystems.expand_coeffs(pol, z) == [2,3,4,0]

        pol = 2z^3+4*z
        @test SymbolicControlSystems.expand_coeffs(pol, z) == [2,0,4,0]

        @vars ω0
        sys = ControlSystems.DemoSystems.resonant(;ω0)
        ssys = Sym(sys)
        n, d = sp.fraction(simplify(ssys))
        @test SymbolicControlSystems.expand_coeffs(n, s) == [1.0*ω0^3]
        @test sum(SymbolicControlSystems.expand_coeffs(d, s) - [1, 0.5, ω0^2+0.25^2]) == 0
    end

    @testset "Sym -> tf and vice versa" begin
        @info "Testing Sym -> tf and vice versa"

        @vars a b c
        @test_both tf([a], [b, c])  a/(b*s + c)
        @test_both tf([a], [b, c], 0.1)  a/(b*z + c)

        @test_both tf([a], [b, c, 1])  a/(b*s^2 + c*s + 1)
        @test_both tf([a], [b, c, 1], 0.1)  a/(b*z^2 + c*z + 1)

        # Test that denominator is still monic
        @test string(Sym(tf([a], [1, c, b]))) == string(a/(s^2 + c*s + b))
        @test string(Sym(tf([a], [1, c, b], 0.1))) == string(a/(z^2 + c*z + b))

        # Test that denominator is still monic
        @test string(Sym(tf([b, a], [1, c, 1]))) == string((b*s + a)/(s^2 + c*s + 1))
        @test string(Sym(tf([b, a], [1, c, 1], 0.1))) == string((b*z + a)/(z^2 + c*z + 1))

        # Test that denominator is still monic
        @test_both tf([a], [1, c, b])  a/(s^2 + c*s + b)
        @test_both tf([a], [1, c, b], 0.1)  a/(z^2 + c*z + b)

        # Test that denominator is still monic
        @test_both tf([b, a], [1, c, 1])  (b*s + a)/(s^2 + c*s + 1)
        @test_both tf([b, a], [1, c, 1], 0.1)  (b*z + a)/(z^2 + c*z + 1)

    end

    @testset "sym2num" begin
        @info "Testing sym2num"
        @vars a b c
        G = tf([a], [b, c])
        @test sym2num(G, a=>1, b=>1, c=>1) == tf(1,[1,1])
        @test sym2num(G, a=>1, b=>2, c=>1) == tf(1,[2,1])
        @test sym2num(G, a=>1, b=>1, c=>2) == tf(1,[1,2])

        G = tf([b, a], [1, c, 1])
        @test sym2num(G, a=>1, b=>1, c=>1) == tf([1,1],[1,1,1])
        @test sym2num(G, a=>1, b=>2, c=>1) == tf([2,1],[1,1,1])
        @test sym2num(G, a=>1, b=>1, c=>2) == tf([1,1],[1,2,1])

        G = tf([a], [b, c], 0.1)
        @test_throws ErrorException sym2num(G, a=>1, b=>1, c=>1) == tf(1,[1,1], 0.1)

        @test sym2num(G, 0.1, a=>1, b=>1, c=>1) == tf(1,[1,1], 0.1)
        @test sym2num(G, 0.1, a=>1, b=>2, c=>1) == tf(1,[2,1], 0.1)
        @test sym2num(G, 0.1, a=>1, b=>1, c=>2) == tf(1,[1,2], 0.1)

        G = tf([b, a], [1, c, 1], 0.1)
        @test sym2num(G, 0.1, a=>1, b=>1, c=>1) == tf([1,1],[1,1,1], 0.1)
        @test sym2num(G, 0.1, a=>1, b=>2, c=>1) == tf([2,1],[1,1,1], 0.1)
        @test sym2num(G, 0.1, a=>1, b=>1, c=>2) == tf([1,1],[1,2,1], 0.1)

    end

    @testset "Tustin and C-code" begin
        @info "Testing Tustin and C-code"

        @vars J c
        G = tf(1.,[J^2,c,1])
        Gs = Sym(G)
        Gn = tf(1.,[1,1,1])
        Gt = tustin(G, 0.1)
        @test denvec(Gt)[1][1] == 1 # test that it's monic
        @test_throws ErrorException sym2num(Gt, J=>1, c=>1)
        Gtn = sym2num(Gt, 0.1, J=>1, c=>1)
        @test minreal(Gtn) ≈ c2d(Gn, 0.1, :tustin)
        code = SymbolicControlSystems.ccode(Gt)
        @test occursin("static double u[3] = {0};", code)
        @test occursin("static double y[3] = {0};", code)
        @test occursin("double ui, double J, double c", code)

        code = SymbolicControlSystems.ccode(ss(Gt), cse=false)
        @test occursin("double u, double J, double c", code)

        code = SymbolicControlSystems.ccode(ss(Gt), cse=true)
        @test occursin("double u, double J, double c", code)

        @vars w T d # Define symbolic variables
        h = 0.01
        G = tf([w^2], [1, 2*d*w, w^2]) * tf(1, [T, 1])
        Gd = tustin(G, h) # Discretize
        code = SymbolicControlSystems.ccode(Gd, cse=true)
        path = mktempdir()
        filename = joinpath(path, "code.c")
        outname = joinpath(path, "test.so")
        write(joinpath(path, filename), code)
        run(`gcc $filename -lm -shared -o $outname`)

        function c_lsim(u, T, d, w)
            Libc.Libdl.dlopen(outname) do lib
                fn = Libc.Libdl.dlsym(lib, :transfer_function)
                map(u) do u
                    @ccall $(fn)(u::Float64, T::Float64, d::Float64, w::Float64)::Float64
                end
            end
        end

        function c_lsim_ss(u, T, d, w)
            y = zeros(1)
            Libc.Libdl.dlopen(outname) do lib
                fn = Libc.Libdl.dlsym(lib, :transfer_function)
                map(u) do u
                    @ccall $(fn)(y::Ref{Cdouble}, u::Float64, T::Float64, d::Float64, w::Float64)::Cvoid
                    y[]
                end
            end
        end

        function c_lsim_ss(u)
            Y = Libc.Libdl.dlopen(outname) do lib
                fn = Libc.Libdl.dlsym(lib, :transfer_function)
                map(u) do u
                    y = zeros(2)
                    @ccall $(fn)(y::Ref{Cdouble}, u::Float64)::Cvoid
                    y
                end
            end
            reduce(hcat, Y)'
        end

        u = randn(1,1000); # Random input signal
        T_, d_, w_ = 0.03, 0.2, 2.0 # Define system parameters
        y = c_lsim( u,  T_,  d_,  w_); # Filter u through the C-function filter
        Gd_ = sym2num(Gd, h, Pair.((T, d, w), (T_, d_, w_))...) # Replace symbols with numeric constants
        y_,_ = lsim(ss(Gd_), u); # Filter using Julia
        @test norm(y-y_)/norm(y_) < 1e-10

        code = SymbolicControlSystems.ccode(ss(Gd), cse=true)
        write(joinpath(path, filename), code)
        run(`gcc $filename -lm -shared -o $outname`)
        y = c_lsim_ss( u,  T_,  d_,  w_); # Filter u through the
        y_,_ = lsim(ss(Gd_), u);
        @test norm(y-y_)/norm(y_) < 1e-10 # TODO: figure out why this is more sensitive

        G_ = sym2num(G, Pair.((T, d, w), (T_, d_, w_))...)
        y_,_ = lsim(c2d(ss(G_), h, :tustin), u);
        @test norm(y-y_)/norm(y_) < 1e-10 # TODO: figure out why this is more sensitive

        # test without symbols

        J = 2
        c = 1
        G = tf(1.,[J^2,c,1])
        Gn = tf(1.,[1,1,1])
        Gt = tustin(G, 0.1)
        code = SymbolicControlSystems.ccode(Gt)
        @test occursin("static double u[3] = {0};", code)
        @test occursin("static double y[3] = {0};", code)
        @test occursin("double ui", code)

        code = SymbolicControlSystems.ccode(ss(Gt), cse=false)
        @test occursin("double u", code)

        code = SymbolicControlSystems.ccode(ss(Gt), cse=true)
        @test occursin("double u", code)

        ## test multi output
        Gt2 = c2d([ss(G); ss(Gn)], 0.1, :tustin)
        code = SymbolicControlSystems.ccode(Gt2, cse=true)
        write(joinpath(path, filename), code)
        run(`gcc $filename -lm -shared -o $outname`)
        @test occursin("double *y, double u", code)

        u = randn(1,1000); # Random input signal
        y_,_ = lsim(Gt2, u); # Filter using Julia
        y = c_lsim_ss(u);
        @test norm(y'-y_)/norm(y_) < 1e-10

    end

    # @testset "symbolics2sym" begin
    #     @info "Testing symbolics2sym"

    #     @variables s m k c b d a
    #     Ld1 = ((-exp(-0.001b) - (2exp(-0.001sqrt(k))))*((0.001 - (5.0e-7c*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (1.0 - (5.0e-7k*(m^-1)))^2) + (exp(-0.001sqrt(k))^2 + 2exp(-0.001b)*exp(-0.001sqrt(k)))*(1.0 - (5.0e-7k*(m^-1))) + ((0.001 - (5.0e-7c*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (1.0 - (5.0e-7k*(m^-1)))^2)*(1.0 - (5.0e-7k*(m^-1))) + ((1.0 - (5.0e-7k*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (5.0e-7c*k*(m^-2) - (0.001k*(m^-1)))*(1.0 + 5.0e-7(c^2)*(m^-2) - (0.001c*(m^-1)) - (5.0e-7k*(m^-1))))*(0.001 - (5.0e-7c*(m^-1))) - (5.0e-10((-exp(-0.001b) - (2exp(-0.001sqrt(k))))*((1.0 - (5.0e-7k*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (5.0e-7c*k*(m^-2) - (0.001k*(m^-1)))*(1.0 + 5.0e-7(c^2)*(m^-2) - (0.001c*(m^-1)) - (5.0e-7k*(m^-1)))) + (exp(-0.001sqrt(k))^2 + 2exp(-0.001b)*exp(-0.001sqrt(k)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + ((0.001 - (5.0e-7c*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (1.0 - (5.0e-7k*(m^-1)))^2)*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + ((1.0 - (5.0e-7k*(m^-1)))*(5.0e-7c*k*(m^-2) - (0.001k*(m^-1))) + (5.0e-7c*k*(m^-2) - (0.001k*(m^-1)))*(1.0 + 5.0e-7(c^2)*(m^-2) - (0.001c*(m^-1)) - (5.0e-7k*(m^-1))))*(1.0 + 5.0e-7(c^2)*(m^-2) - (0.001c*(m^-1)) - (5.0e-7k*(m^-1))))*((1.0e-6 - (5.0e-10b) - (5.0e-10c*(m^-1)))^-1)) - (exp(-0.001b)*(exp(-0.001sqrt(k))^2)))*((1.0e-9 + 2.5e-16(b^2) + 0.001(0.001 - (5.0e-7c*(m^-1)))*(0.001 - (5.0e-7b) - (5.0e-7c*(m^-1))) - (5.0e-13b) - (2.5e-16k*(m^-1)) - (5.0e-10(0.001(0.001 - (5.0e-7b) - (5.0e-7c*(m^-1)))*(1 + 5.0e-7(b^2) - (0.001b)) + 0.001(0.001 - (5.0e-7b) - (5.0e-7c*(m^-1)))*(1.0 + 5.0e-7(c^2)*(m^-2) - (0.001c*(m^-1)) - (5.0e-7k*(m^-1))) + 2.5e-16c*k*(m^-2) - (5.0e-13k*(m^-1)))*((1.0e-6 - (5.0e-10b) - (5.0e-10c*(m^-1)))^-1)))^-1)
    #     Ldsp = symbolics2sym(Ld1)
    #     symbolicsval = Symbolics.substitute(Ld1, Dict(b=>1.0, k=>3.0, c=>5.0, m=>7.0))

    #     @vars s m k c b d a
    #     for (sym,val) in [(b,1.0), (k,3.0), (c,5.0), (m,7.0)]
    #         Ldsp = SymPy.subs(Ldsp, (sym, val))
    #     end
    #     @test float(Ldsp) ≈ symbolicsval.val rtol=1e-6
    # end
end


text:
using URIParser, MbedTLS, WebSockets

"""
# Engine
"""
type Engine
    router::Router
end

# right now, I'll just use the default Mux
Engine()=Engine(Pivot.TRouter())

handle!(fn::Handler,
        e::Engine,
        method::Verb,
        path::AbstractString) = handle!(fn, e.router, method, parseurl(path))

handle!(e::Engine,
        fn::Handler,
        method::Verb,
        path::AbstractString) = handle!(fn, e.router, method, parseurl(path))

use!(e::Engine, fn::Function) = use!(e.router, fn)
use!(fn::Function, e::Engine) = use!(e.router, fn)
use!(e::Engine, fns::Function...) = map((fn) -> use!(e, fn), fns)

"""
# buildhandler

returns a function that is considered the app of mux middleware.

This function fetched the proper endpoint, along with the proper handler, then passes the context through the middleware
"""
function buildhandler(router::Router)
    function (ctx)
        ep, params = fetch(router, ctx.data[:request][:path])
        method = STI[ctx.data[:request][:method]]
        ctx.data[:params] = params
        ctx |> proper_method(method, ep.handlermap)
    end
end

"""
# run

runs the server
"""
function run(e::Engine, port::Number=8080; keyspath="", reference=Dict())
    Pivot.finalize!(e.router)
    reference[:http] = HttpHandler() do req::Request, res::Response
        # need to parse the path
        rq = req |> Mux.todict |> Pivot.splitquery
        ctx = Context(req)
        ctx.data[:router] = e.router
        ctx.data[:request] = rq
        ctx.data[:reference] = reference
        str_resp = mux(e.router.middleware, buildhandler(e.router))(ctx)
        ctx.response.data = convert(Vector{UInt8}, str_resp)
        Response(ctx)
    end

    reference[:socket] = WebSocketHandler() do req, client
        while true
            msg = read(client) |> utf8
            write(client, msg)
        end
    end

    reference[:server] = Server(reference[:http],
                                reference[:socket])

    # ssl for https
    # TODO: test remotely
    if keyspath != ""
        rel = create_rp(keyspath)
        generate_ssl(rel(""))

        cert = MbedTLS.crt_parse_file(rel("keys/server.crt"))
        key = MbedTLS.parse_keyfile(rel("keys/server.key"))
        HttpServer.run(reference[:server], port=port, ssl=(cert, key))
        return
    end

    HttpServer.run(reference[:server], port)
end


text:
# module TestRclBase

using Rcl
using Test

@testset "init and shutdown" begin
    # uninitialized
    @test rcl.ok() == false
    @test rcl.instance_id() == 0

    # initialize
    rcl.init()
    @test rcl.ok() == true
    @test rcl.instance_id() != 0

    # shut down
    rcl.shutdown()
    @test rcl.ok() == false
    @test rcl.instance_id() == 0

    # init again
    rcl.init()
    @test rcl.ok() == true
    @test rcl.instance_id() != 0

    # shut down again
    rcl.shutdown()
    @test rcl.ok() == false
    @test rcl.instance_id() == 0

    @testset "double init/shutdown exceptions" begin
        # uninitialized
        @test !rcl.ok()

        # Initialize
        rcl.init()
        @test rcl.ok()

        # Try to initialize again
        @test_throws RclError rcl.init()

        # Still initialized
        @test rcl.ok()

        # shutdown
        rcl.shutdown()
        @test !rcl.ok()

        # Try to shut down again
        @test_throws RclError rcl.shutdown()
        @test !rcl.ok()
    end
end

# end # module TestRclBase
nothing


text:
using HPAT
using HPAT.API.LinearRegression
using MPI

@acc hpat function calcLinearRegression(file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points", file_name)
    responses = DataSource(Vector{Float64},HDF5,"/responses", file_name)
    coeffs = LinearRegression(points, responses)
    return coeffs
end

c2 = calcLinearRegression(HPAT.getDefaultDataPath()*"/linear_regression_train_large.hdf5")

if MPI.Comm_rank(MPI.COMM_WORLD)==0	println("linear regression Coeffs: ", c2) end


text:
const CINT_SIZE = sizeof(Cint)

function writeint(mem::Memory, int::Cint)
    buf = reinterpret(UInt8, [int])

    checkbounds(mem, CINT_SIZE)
    unsafe_copyto!(mem.ptr, pointer(buf), CINT_SIZE)
    return CINT_SIZE
end

function readint(mem::Memory)
    buf = Vector{UInt8}(undef, CINT_SIZE)

    checkbounds(mem, CINT_SIZE)
    unsafe_copyto!(pointer(buf), mem.ptr, CINT_SIZE)
    return reinterpret(Cint, buf)[1]
end

mutable struct SimpleDoubleBuffer
    buffer::Array{UInt8, 2}
    next::Bool
end

SimpleDoubleBuffer(buf_size::Integer) = SimpleDoubleBuffer(Array{UInt8}(undef, 2, buf_size), false)
function get_buffer!(db::SimpleDoubleBuffer)
    out_buffer = db.buffer[db.next+1, :]
    db.next = !db.next  # Update index
    return out_buffer
end

function max_compressed_size(in_size::Integer, block_size::Integer)
    num_blocks = ceil(Int, in_size / block_size)
    compressed_size = LZ4_compressBound(block_size) + CINT_SIZE

    return compressed_size * num_blocks
end

mutable struct LZ4FastCompressor <: TranscodingStreams.Codec
    streamptr::Ptr{LZ4_stream_t}
    acceleration::Cint
    block_size::Int
    buffer::SimpleDoubleBuffer
end

"""
    LZ4FastCompressor(; kwargs...)

Creates an LZ4 compression codec.

# Keywords
- `acceleration::Integer=0`: acceleration factor
- `block_size::Integer=1024`: The size in bytes to encrypt into each block.
"""
function LZ4FastCompressor(; acceleration::Integer=0, block_size::Integer=1024)
    if block_size > LZ4_MAX_INPUT_SIZE
        throw(ArgumentError("`block_size` larger than $LZ4_MAX_INPUT_SIZE."))
    end

    return LZ4FastCompressor(C_NULL, acceleration, block_size, SimpleDoubleBuffer(block_size))
end

const LZ4FastCompressorStream{S} = TranscodingStream{LZ4FastCompressor,S} where S<:IO

"""
    LZ4FastCompressor(stream::IO; kwargs...)

Creates an LZ4 compression stream. See `LZ4FastCompressor()` and `TranscodingStream()` for arguments.
"""
function LZ4FastCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:acceleration, :block_size))
    return TranscodingStream(LZ4FastCompressor(; x...), stream; y...)
end

"""
    TranscodingStreams.expectedsize(codec::Union{LZ4FastCompressor, LZ4HCCompressor}, input::Memory)

Returns the expected size of the transcoded data.
"""
function TranscodingStreams.expectedsize(codec::LZ4FastCompressor, input::Memory)::Int
    max_compressed_size(length(input), codec.block_size)
end

"""
   TranscodingStreams.minoutsize(codec::Union{LZ4FastCompressor, LZ4HCCompressor}, input::Memory)

Returns the minimum output size of `process`.
"""
function TranscodingStreams.minoutsize(codec::LZ4FastCompressor, input::Memory)::Int
    LZ4_compressBound(input.size) + CINT_SIZE
end

"""
   TranscodingStreams.initialize(codec::LZ4FastCompressor)

Initializes the LZ4 Compression Codec.
"""
function TranscodingStreams.initialize(codec::LZ4FastCompressor)::Nothing
    codec.streamptr = LZ4_createStream()
    nothing
end

"""
    TranscodingStreams.finalize(codec::LZ4FastCompressor)

Finalizes the LZ4 Compression Codec.
"""
function TranscodingStreams.finalize(codec::LZ4FastCompressor)::Nothing
    LZ4_freeStream(codec.streamptr)
    codec.streamptr = C_NULL
    nothing
end

"""
    TranscodingStreams.startproc(codec::LZ4FastCompressor, mode::Symbol, error::Error)

Starts processing with the codec
"""
function TranscodingStreams.startproc(
    codec::LZ4FastCompressor,
    mode::Symbol,
    error::Error
)::Symbol

    try
        LZ4_resetStream(codec.streamptr)
        :ok
    catch err
        error[] = err
        :error
    end
end

"""
    TranscodingStreams.process(codec::LZ4FastCompressor, input::Memory, output::Memory, error::Error)

Compresses the data from `input` and writes to `output`.
The process follows the compression example from
https://github.com/lz4/lz4/blob/dev/examples/blockStreaming_doubleBuffer.c
wherein each block of encoded data is prefixed by the number of bytes contained in that block.
Decoding can be done through the LZ4SafeDecompressor.
"""
function TranscodingStreams.process(
    codec::LZ4FastCompressor,
    input::Memory,
    output::Memory,
    error::Error
)::Tuple{Int,Int,Symbol}

    input.size == 0 && return (0, 0, :end)
    try
        in_buffer = pointer(get_buffer!(codec.buffer))

        data_size = min(input.size, codec.block_size)
        out_buffer = Vector{UInt8}(undef, LZ4_compressBound(data_size))
        unsafe_copyto!(in_buffer, input.ptr, data_size)

        compressed_size = LZ4_compress_fast_continue(
            codec.streamptr,
            in_buffer,
            pointer(out_buffer),
            data_size,
            length(out_buffer),
            codec.acceleration,
        )

        checkbounds(output, compressed_size + CINT_SIZE)
        writeint(output, compressed_size)
        unsafe_copyto!(output.ptr + CINT_SIZE, pointer(out_buffer), compressed_size)

        return (data_size, compressed_size + CINT_SIZE, :ok)
    catch err
        error[] = err
        return (0, 0, :error)
    end
end

mutable struct LZ4SafeDecompressor <: TranscodingStreams.Codec
    streamptr::Ptr{LZ4_streamDecode_t}
    block_size::Int
    buffer::SimpleDoubleBuffer
end

"""
    LZ4SafeDecompressor(; kwargs...)

Creates an LZ4 compression codec.

# Keywords
- `block_size::Integer=1024`: The size in bytes of unecrypted data contained in each block.
    Must match or exceed the compression `block_size` or decompression will fail.
"""
function LZ4SafeDecompressor(; block_size::Integer=1024)
    return LZ4SafeDecompressor(C_NULL, block_size, SimpleDoubleBuffer(block_size))
end

const LZ4SafeDecompressorStream{S} = TranscodingStream{LZ4SafeDecompressor,S} where S<:IO

"""
    LZ4SafeDecompressorStream(stream::IO; kwargs...)

Creates an LZ4 compression stream. See `LZ4SafeDecompressor()` and `TranscodingStream()` for arguments.
"""
function LZ4SafeDecompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:block_size,))
    return TranscodingStream(LZ4SafeDecompressor(; x...), stream; y...)
end

"""
    TranscodingStreams.expectedsize(codec::LZ4SafeDecompressor, input::Memory)

Returns the expected size of the transcoded data.
"""
function TranscodingStreams.expectedsize(codec::LZ4SafeDecompressor, input::Memory)::Int
    max(input.size * 2, codec.block_size)
end

"""
   TranscodingStreams.minoutsize(codec::LZ4SafeDecompressor, input::Memory)

Returns the minimum output size of `process`.
"""
function TranscodingStreams.minoutsize(codec::LZ4SafeDecompressor, input::Memory)::Int
    max(input.size * 2, codec.block_size)
end

"""
   TranscodingStreams.initialize(codec::LZ4SafeDecompressor)

Initializes the LZ4 Compression Codec.
"""
function TranscodingStreams.initialize(codec::LZ4SafeDecompressor)::Nothing
    codec.streamptr = LZ4_createStreamDecode()
    nothing
end

"""
    TranscodingStreams.finalize(codec::LZ4SafeDecompressor)

Finalizes the LZ4SafeDecompression Codec.
"""
function TranscodingStreams.finalize(codec::LZ4SafeDecompressor)::Nothing
    LZ4_freeStreamDecode(codec.streamptr)
    codec.streamptr = C_NULL
    nothing
end

"""
    TranscodingStreams.startproc(codec::LZ4SafeDecompressor, mode::Symbol, error::Error)

Starts processing with the codec
"""
function TranscodingStreams.startproc(
    codec::LZ4SafeDecompressor,
    mode::Symbol,
    error::Error
)::Symbol

    try
        LZ4_setStreamDecode(codec.streamptr)
        :ok
    catch err
        error[] = err
        :error
    end
end

"""
    TranscodingStreams.process(codec::LZ4SafeDecompressor, input::Memory, output::Memory, error::Error)

Compresses the data from `input` and writes to `output`.
The process follows the decompression example from
https://github.com/lz4/lz4/blob/dev/examples/blockStreaming_doubleBuffer.c
and requires each block of encoded data to be prefixed by the number of bytes contained in that block.
"""
function TranscodingStreams.process(
    codec::LZ4SafeDecompressor,
    input::Memory,
    output::Memory,
    error::Error
)::Tuple{Int,Int,Symbol}

    input.size == 0 && return (0, 0, :end)
    try
        out_buffer = get_buffer!(codec.buffer)
        data_size = readint(input)

        decompressed_size = LZ4_decompress_safe_continue(
            codec.streamptr,
            input.ptr+CINT_SIZE,
            pointer(out_buffer),
            data_size,
            length(out_buffer)
        )

        checkbounds(output, decompressed_size)
        unsafe_copyto!(output.ptr, pointer(out_buffer), decompressed_size)

        return (data_size + CINT_SIZE, decompressed_size, :ok)
    catch err
        error[] = err
        return (0, 0, :error)
    end
end