text:
"""
    @autodef struct_def

Sets up methods to automatically evaluate the Jacobian of the function. By default, it defines
these methods for a new type `NewFun`:

    jacobian!(::InPlace, ::ForwardAD, fun::NewFun, J, y, z)
    jacobian!(::InPlace, ::FiniteDifference, fun::NewFun, J, y, z)
    jacobian!(::StaticReturn, ::ForwardAD, fun::NewFun, J, y, z)
    jacobian!(::StaticReturn, ::FiniteDifference, fun::NewFun, J, y, z)

These methods are all optimized to be non-allocating.

The method will modify the type definition, adding fields and type parameters to set them
up for efficient evaluation of the jacobians. The changes are usually transparent to the user.

If an inner constructor is provided, the signature will be unchanged, but will be modified
to initialize the new fields and provide the new type parameters.

# Limitations
* `RobotDynamics` must be defined the local module (cannot be hidden by an alias)
* If the type has type parameters and no inner constructor, the type parameters must be
  explicitly provided when constructing the object. No automatic inference will be done.
  It's best to explicitly provide an inner constructor with the desired behavior.

# Examples

```julia
@autodef struct MyFun <: RobotDynamics.AbstractFunction end
```

will define
```julia
struct MyFun{CH} <: RobotDynamics.AbstractFunction
    cfg::ForwardDiff.JacobianConfig{Nothing, Float64, CH, Tuple{Vector{ForwardDiff.Dual{Nothing, Float64, CH}}, Vector{(ForwardDiff).Dual{Nothing, Float64, CH}}}}
    cache::FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}
    function MyFun()
        model = new{0}()
        _n = input_dim(model)
        _m = output_dim(model)
        cfg = ForwardDiff.JacobianConfig(nothing, zeros(Float64, _m), zeros(Float64, _n))
        model = new{length(cfg).seeds}(cfg)
        _n = input_dim(model)
        _m = output_dim(model)
        cache = FiniteDiff.JacobianCache(zeros(Float64, _n), zeros(Float64, _m))
        model = new{length(cfg).seeds}(cfg, cache)
    end
end

function RobotDynamics.jacobian!(::InPlace, ::ForwardAD, fun::MyFun, J, y, z)
    f_aug!(y, z_) = begin
        RobotDynamics.evaluate!(fun, y, RobotDynamics.setinputs!(z, z_))
        ForwardDiff.jacobian!(J, f_aug!, y, z, fun.cfg)
    end
end
...  # other Jacobian methods
```

"""
macro autodiff(type::Expr)
    _autodiff(__module__, type, [InPlace(), StaticReturn()], [ForwardAD(), FiniteDifference()])
end

macro autodiff(type::Expr, args...)
    sigs = FunctionSignature[]
    methods = DiffMethod[]
    for arg in args
        if arg isa Symbol
            earg = eval(arg)
            if earg <: FunctionSignature
                push!(sigs, earg())
            elseif earg <: DiffMethod
                push!(methods, earg())
            end
        end
    end
    _autodiff(__module__, type, sigs, methods)
end

function _autodiff(mod, struct_expr::Expr, sigs::Vector{FunctionSignature}, methods::Vector{DiffMethod})
    jac_defs = Expr[]

    # Get signature for jacobian methods
    type_expr = struct_expr.args[2].args[1]

    # Check if it's a subtype of DiscreteDynamics
    parent_expr = get_struct_parent(struct_expr)
    parent = mod.eval(parent_expr)
    is_discrete_dynamics = parent <: RobotDynamics.DiscreteDynamics

    # Check if it's a subtype of ScalarFunction
    is_scalar_fun = parent <: RobotDynamics.ScalarFunction

    for method in methods
        for sig in sigs
            if !is_scalar_fun
                push!(jac_defs, gen_jacobian(sig, method, type_expr, mod))
                if is_discrete_dynamics
                    push!(jac_defs, gen_dynamics_jacobian(sig, method, type_expr, mod))
                end
            end
        end
        struct_expr = modify_struct_def(method, struct_expr, mod, is_scalar_fun)

        if is_scalar_fun
            gradfun, hessfun = gen_hessgrad(method, type_expr, mod)
            push!(jac_defs, gradfun)
            push!(jac_defs, hessfun)
        end
    end
    return quote
        $struct_expr
        $(Expr(:block, jac_defs...))
    end
end

function get_call_signature(sig, diff, fname, fargs, type_expr, mod)
    # Get the name of the type
    local type_name
    if type_expr isa Expr
        if type_expr.head == :curly
            type_name = type_expr.args[1]
        else
            display(type_expr)
            error("Got an unexpected type expression. Can't generate the Jacobian.")
        end
    elseif type_expr isa Symbol
        type_name = type_expr
    end

    # Convert type to a GlobalRef
    (type_name isa Symbol) && (type_name = GlobalRef(mod, type_name))

    # Create the call signature
    if isnothing(sig)
        Expr(:call, GlobalRef(@__MODULE__, fname), :(::$(typeof(diff))), :(fun::$type_name), fargs...)
    else
        Expr(:call, GlobalRef(@__MODULE__, fname), :(::$(typeof(sig))), :(::$(typeof(diff))), :(fun::$type_name), fargs...)
    end
end

"""
Takes an inner constructor of the form:
```
function MyType(args...)
    ...
    new(args...)
end
```
and adds type params:
```
function MyType{params...}(args...) where {params...}
    ...
    new{params...}(args...)
end
```

"""
function add_params_to_innercon(fun, params)
    @assert fun.head == :function
    sig = fun.args[1]
    newsig = :($(copy(sig)) where {$(params...)})
    typename = newsig.args[1].args[1]
    newsig.args[1].args[1] = :($typename{$(params...)})

    newfun = copy(fun)
    newfun.args[1] = newsig
    newcall = newfun.args[2].args[end]
    newcall.args[1] = :(new{$(params...)})
    newfun
end

function get_struct_parent(struct_expr::Expr)
    @assert struct_expr.head == :struct
    typedef = struct_expr.args[2]
    if typedef.head != :(<:)
        error("Struct has no parent!")
    else
        parent_name = typedef.args[2]
        if parent_name isa Expr && parent_name.head == :curly
            parent_params = parent_name.args[2:end]
            parent_name = Expr(:where, parent_name, parent_params...)
        end
        return parent_name
    end
end

function get_parent_name(parent, mod)
    loc = nothing
    if parent isa Symbol
        return (parent, loc)
    end
    name_w_params = parent
    if parent isa Expr && parent.head == :where
        name_w_params = parent.args[1]
        loc = parent
    end

    name = name_w_params
    if name_w_params isa Expr && name_w_params.head == :curly
        name = name_w_params.args[1]
        loc = name_w_params
    end

    if name isa GlobalRef || name isa Symbol
        return (name, loc)
    end

    if name isa Expr && name.head == :.
        return (name, loc)
        # return (GlobalRef(mod.eval(name.args[1]), name.args[2].value), loc)
    end
    error("Couldn't get parent name")
end

function add_field_to_struct(struct_expr0::Expr, newfield::Vector{Expr}, init_field, pname::Union{Symbol,Nothing}, init_param, mod)

    struct_expr = copy(struct_expr0)
    @assert struct_expr.head == :struct
    typedef = struct_expr.args[2]
    body = struct_expr.args[3]

    # Check for valid sub-typing
    parent_expr = get_struct_parent(struct_expr)
    parent = mod.eval(parent_expr)
    if !(parent <: RobotDynamics.AbstractFunction)
        error("Type must be a sub-type of RobotDynamics.AbstractFunction")
    end
    type_param = Symbol(inputtype(parent))

    # Resolve the parent name in the original scope
    parent_name, loc = get_parent_name(parent_expr, mod)
    if !(parent_name isa GlobalRef)
        local gref
        if parent_name isa Symbol
            gref = GlobalRef(mod, parent_name)
        else
            typename = parent_name.args[2].value
            parent_module_symb = parent_name.args[1]
            parent_module = mod.eval(parent_module_symb)
            gref = GlobalRef(parent_module, parent_name.args[2].value)
        end
        if isnothing(loc)
            typedef.args[2] = gref
        else
            loc.args[1] = gref
        end
    end
    # @show parent_name.mod
    # @show parent_name

    # Add type parameter
    has_params = false
    if !isnothing(pname)
        if typedef.args[1] isa Expr && typedef.args[1].head == :curly
            has_params = true
            typeparams = typedef.args[1].args

            # Check if the parameter name alraedy exists
            for param in typeparams
                if (param isa Symbol && param == pname) || (param isa GlobalRef && param.name == pname)
                    error("Cannot add parameter, the type parameter $pname is already being used.")
                end
            end

            # Add to type params
            push!(typeparams, pname)
        elseif typedef.args[1] isa Symbol  # add new param
            typename = typedef.args[1]
            typedef.args[1] = :($typename{$pname})
        end
    end

    # Add field
    names = map(newfield) do field
        @assert field isa Expr && field.head == :(::)
        name = field.args[1]

        # Check if the field name already exists
        for field in body.args
            if field isa Expr && field.head == :(::) && field.args[1] == name
                error("Cannot add field, the field $name already exists")
            end
        end

        # Add the field to the end of the struct
        push!(body.args, field)

        name
    end

    # Modify the constructor
    constructor_found = false
    for field in body.args
        if field isa Expr && field.head == :function
            constructor_found |= add_config_to_constructor(field, type_param, names, init_field, init_param)
        end
    end

    # Add new constructor, if needed
    if !constructor_found
        inner_con = new_default_constructor(struct_expr0)
        push!(body.args, inner_con)
        add_config_to_constructor(inner_con, type_param, names, init_field, init_param)
    end
    return struct_expr
end

function new_default_constructor(struct_expr::Expr)
    @assert struct_expr.head == :struct
    struct_expr = copy(struct_expr)
    typedef = struct_expr.args[2]
    body = struct_expr.args[3]

    # Get the type name
    if typedef.head == :(<:)
        typename = typedef.args[1]
    else
        typename = typedef
    end

    # Get parameters
    params = []
    if typename isa Expr && typename.head == :curly
        params = typename.args[2:end]
    end

    # Get the field names
    args = map(filter(body.args) do arg
        (arg isa Expr && arg.head == :(::)) || arg isa Symbol
    end) do arg
        arg isa Expr ? arg.args[1] : arg
    end

    # Call to "new"
    inner_con = quote
        function $typename($(args...))
            new($(args...))
        end
    end
    inner_con = inner_con.args[2]
    @assert inner_con.head == :function

    # Add type params
    if !isempty(params)
        # Add "where" clause
        call = inner_con.args[1]
        inner_con.args[1] = Expr(:where, call, params...)

        # Add to "new"
        newcall = inner_con.args[2].args[end]
        newcall.args[1] = :(new{$(params...)})
    end

    return inner_con
end

function add_config_to_constructor(con, type_param, fieldname, init_field, init_param)
    @assert con.head == :function
    is_constructor = false
    constructor_body = con.args[2]
    for (j, line) in enumerate(reverse(constructor_body.args))
        i = length(constructor_body.args) - j + 1  # line number

        isnewcall(call) = call.args[1] == :new || (call.args[1] isa Expr && call.args[1].head == :curly && call.args[1].args[1] == :new)
        hasnewcall = line isa Expr && line.head == :call && isnewcall(line)
        hasnewassign = line isa Expr && line.head == :(=) && line.args[2] isa Expr && line.args[2].head == :call && isnewcall(line.args[2])
        if hasnewcall
            callfun = line
        elseif hasnewassign
            callfun = line.args[2]
        else
            continue
        end

        newparam = init_param(fieldname)
        add_newparam = !isnothing(newparam)

        # Add new param, so that the "new" call looks like
        # new{...,new_param}(...)
        if add_newparam
            @assert callfun.head == :call
            hasparams = callfun.args[1] isa Expr
            new_param0 = :(0)
            if hasparams  # append type param
                @assert callfun.args[1].head == :curly
                push!(callfun.args[1].args, new_param0)
            else  # add new type param
                @assert callfun.args[1] == :new
                callfun.args[1] = :(new{$new_param0})
            end
        end

        # Append assignment to make the line look like
        # model = new{..., new_param}(...)
        outname = :__model__
        if hasnewcall
            line = :($outname = $line)
            constructor_body.args[i] = line
        else
            # if a model is output, create an alias with the given name
            @assert line.head == :(=)
            if outname != line.args[1]
                insert!(constructor_body.args, i+1, :($outname = $(line.args[1])))
            end
        end

        # add the new argument to "new"
        callfun_new = copy(callfun)
        if add_newparam
            callfun_new.args[1].args[end] = newparam
        end
        append!(callfun_new.args, fieldname)  # add to the arguments to "new"

        # add lines after the call to "new" that initialize the new variables
        init_expr = init_field(outname, fieldname, type_param, callfun_new)
        for arg in reverse(init_expr.args)
            insert!(constructor_body.args, i+1, arg)
        end
        is_constructor = true
        break
    end
    return is_constructor
end

#########################
# ForwardDiff
#########################

function gen_jacobian(sig::StaticReturn, diff::ForwardAD, type_expr, mod)
    fname = :jacobian!
    fargs = (:J, :y, :z)
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :evaluate)
    fun_body = quote
        f(_z) = $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        J .= ForwardDiff.jacobian(f, getdata(z))
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function gen_jacobian(sig::InPlace, diff::ForwardAD, type_expr, mod)
    fname = :jacobian!
    fargs = (:J, :y, :z)
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :(evaluate!))
    fun_body = quote
        f!(_y,_z) = $eval(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
        ForwardDiff.jacobian!(J, f!, y, getdata(z), fun.cfg)
        return nothing
    end
    jacfun = Expr(:function, callsig, fun_body)
end

function gen_hessgrad(diff::ForwardAD, type_expr, mod)
    fname = :gradient!
    fargs = (:grad, :z)
    callsig = get_call_signature(nothing, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :evaluate)
    fun_body = quote
        f(_z) = $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        ForwardDiff.gradient!(grad, f, getdata(z), fun.gradcfg)
        return nothing
    end
    gradfun = Expr(:function, callsig, fun_body)

    fname = :hessian!
    fargs = (:hess, :z)
    callsig = get_call_signature(nothing, diff, fname, fargs, type_expr, mod)
    fun_body = quote
        f(_z) = $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        ForwardDiff.hessian!(hess, f, getdata(z), fun.hesscfg)
        return nothing
    end
    hessfun = Expr(:function, callsig, fun_body)
    return gradfun, hessfun
end

function gen_dynamics_jacobian(sig::StaticReturn, diff::ForwardAD, type_expr, mod)
    fname = :dynamics_error_jacobian!
    fargs = (:J2, :J1, :y2, :y1, :(z2::AbstractKnotPoint), :(z1::AbstractKnotPoint))
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :dynamics_error)
    fun_body = quote
        f1(_z) = $eval(fun, z2, RobotDynamics.StaticKnotPoint(z1, _z))
        J1 .= ForwardDiff.jacobian(f1, RobotDynamics.getdata(z1))

        f2(_z) = $eval(fun, RobotDynamics.StaticKnotPoint(z2, _z), z1)
        J2 .= ForwardDiff.jacobian(f2, RobotDynamics.getdata(z2))
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function gen_dynamics_jacobian(sig::InPlace, diff::ForwardAD, type_expr, mod)
    fname = :dynamics_error_jacobian!
    fargs = (:J2, :J1, :y2, :y1, :(z2::AbstractKnotPoint), :(z1::AbstractKnotPoint))
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :dynamics_error!)
    fun_body = quote
        f1!(_y, _z) = $eval(fun, _y, y1, z2, RobotDynamics.StaticKnotPoint(z1, _z))
        ForwardDiff.jacobian!(J1, f1!, y2, RobotDynamics.getdata(z1), fun.cfg)

        f2!(_y, _z) = $eval(fun, _y, y1, RobotDynamics.StaticKnotPoint(z2, _z), z1)
        ForwardDiff.jacobian!(J2, f2!, y2, RobotDynamics.getdata(z2), fun.cfg)
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function init_cfg(outname, fieldname, type_param, callfun)
    if length(fieldname) == 1
        init = [
            :($(fieldname[1]) = ForwardDiff.JacobianConfig(nothing, zeros($type_param, _m), zeros($type_param, _n)))
        ]
    else
        init = [
            :($(fieldname[1]) = ForwardDiff.GradientConfig(nothing, zeros($type_param, _n)))
            :($(fieldname[2]) = ForwardDiff.HessianConfig(nothing, zeros($type_param, _n)))
        ]
    end
    quote
        _n = RobotDynamics.state_dim($outname) + RobotDynamics.control_dim($outname)
        _m = RobotDynamics.output_dim($outname)
        $(init...)
        $outname = $callfun
    end
end

function init_cfg_param(fieldname)
    if length(fieldname) == 1
        :(length($(fieldname[1]).seeds))
    else
        :(length($(fieldname[1]).seeds))
    end
end

function modify_struct_def(::ForwardAD, struct_expr::Expr, mod, is_scalar_fun)
    pname = :JCH
    parent_name = get_struct_parent(struct_expr)
    parent = mod.eval(parent_name)
    type_param = Symbol(inputtype(parent))
    if is_scalar_fun
        newfield = [
            :(gradcfg::ForwardDiff.GradientConfig{Nothing, Float64, JCH, Vector{ForwardDiff.Dual{Nothing, Float64, JCH}}})
            :(hesscfg::ForwardDiff.HessianConfig{Nothing, Float64, JCH, Vector{ForwardDiff.Dual{Nothing, ForwardDiff.Dual{Nothing, Float64, JCH}, JCH}}, Vector{ForwardDiff.Dual{Nothing, Float64, JCH}}})
        ]
    else
        newfield = [
            :(cfg::ForwardDiff.JacobianConfig{Nothing, $type_param, $pname, Tuple{Vector{ForwardDiff.Dual{Nothing, $type_param, $pname}}, Vector{ForwardDiff.Dual{Nothing, $type_param, $pname}}}})
        ]
    end

    struct_expr = copy(struct_expr)
    @assert struct_expr.head == :struct

    return add_field_to_struct(struct_expr, newfield, init_cfg, pname, init_cfg_param, mod)
end

#########################
# FiniteDiff
#########################

function gen_jacobian(sig::StaticReturn, diff::FiniteDifference, type_expr, mod)
    fname = :jacobian!
    fargs = (:J, :y, :z)
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :evaluate)
    fun_body = quote
        f!(_y,_z) = _y .= $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function gen_jacobian(sig::InPlace, diff::FiniteDifference, type_expr, mod)
    fname = :jacobian!
    fargs = (:J, :y, :z)
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :evaluate!)
    fun_body = quote
        f!(_y,_z) = $eval(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
        FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function gen_hessgrad(diff::FiniteDifference, type_expr, mod)
    fname = :gradient!
    fargs = (:grad, :z)
    callsig = get_call_signature(nothing, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :evaluate)
    fun_body = quote
        f(_z) = $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        FiniteDiff.finite_difference_gradient!(grad, f, getdata(z), fun.gradcache)
        return nothing
    end
    gradfun = Expr(:function, callsig, fun_body)

    fname = :hessian!
    fargs = (:hess, :z)
    callsig = get_call_signature(nothing, diff, fname, fargs, type_expr, mod)
    fun_body = quote
        cache = fun.hesscache
        cache.xmm .= z; cache.xmp .= z; cache.xpm .= z; cache.xpp .= z
        f(_z) = $eval(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
        FiniteDiff.finite_difference_hessian!(hess, f, getdata(z), cache)
        return nothing
    end
    hessfun = Expr(:function, callsig, fun_body)
    return gradfun, hessfun
end

function gen_dynamics_jacobian(sig::StaticReturn, diff::FiniteDifference, type_expr, mod)
    fname = :dynamics_error_jacobian!
    fargs = (:J2, :J1, :y2, :y1, :(z2::AbstractKnotPoint), :(z1::AbstractKnotPoint))
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :dynamics_error)
    fun_body = quote
        f1!(_y, _z) = _y .= $eval(fun, z2, RobotDynamics.StaticKnotPoint(z1, _z))
        FiniteDiff.finite_difference_jacobian!(J1, f1!, RobotDynamics.getdata(z1), fun.cache)

        f2!(_y, _z) = _y .= $eval(fun, RobotDynamics.StaticKnotPoint(z2, _z), z1)
        FiniteDiff.finite_difference_jacobian!(J2, f2!, RobotDynamics.getdata(z2), fun.cache)
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function gen_dynamics_jacobian(sig::InPlace, diff::FiniteDifference, type_expr, mod)
    fname = :dynamics_error_jacobian!
    fargs = (:J2, :J1, :y2, :y1, :(z2::AbstractKnotPoint), :(z1::AbstractKnotPoint))
    callsig = get_call_signature(sig, diff, fname, fargs, type_expr, mod)

    eval = GlobalRef(@__MODULE__, :dynamics_error!)
    fun_body = quote
        f1!(_y, _z) = $eval(fun, _y, y1, z2, RobotDynamics.StaticKnotPoint(z1, _z))
        FiniteDiff.finite_difference_jacobian!(J1, f1!, RobotDynamics.getdata(z1), fun.cache)

        f2!(_y, _z) = $eval(fun, _y, y1, RobotDynamics.StaticKnotPoint(z2, _z), z1)
        FiniteDiff.finite_difference_jacobian!(J2, f2!, RobotDynamics.getdata(z2), fun.cache)
        return nothing
    end
    Expr(:function, callsig, fun_body)
end

function init_cache(outname, fieldname, type_param, callfun)
    if length(fieldname) == 1
        init = [
            :($(fieldname[1]) = FiniteDiff.JacobianCache(zeros($type_param, _n), zeros($type_param, _m)))
        ]
    else
        init = [
            :($(fieldname[1]) = FiniteDiff.GradientCache(zeros(_n), zeros(_n), Val(:forward)))
            :($(fieldname[2]) = FiniteDiff.HessianCache(zeros(_n), zeros(_n), zeros(_n), zeros(_n), Val(:hcentral), Val(true)))
        ]
    end
    quote
        _n = RobotDynamics.state_dim($outname) + RobotDynamics.control_dim($outname)
        _m = RobotDynamics.output_dim($outname)
        $(init...)
        $outname = $callfun
    end
end

init_cache_param(fieldname) = nothing

function modify_struct_def(::FiniteDifference, struct_expr::Expr, mod, is_scalar_fun)
    pname = nothing
    parent_name = get_struct_parent(struct_expr)
    parent = mod.eval(parent_name)
    type_param = Symbol(inputtype(parent))
    if is_scalar_fun
        newfield = [
            :(gradcache::FiniteDiff.GradientCache{Nothing, Nothing, Nothing, Vector{Float64}, Val{:forward}(), Float64, Val{true}()})
            :(hesscache::FiniteDiff.HessianCache{Vector{Float64}, Val{:hcentral}(), Val{true}()})
        ]
    else
        newfield = [
            :(cache::FiniteDiff.JacobianCache{Vector{$type_param}, Vector{$type_param}, Vector{$type_param}, UnitRange{Int64}, Nothing, Val{:forward}(), $type_param})
        ]
    end

    struct_expr = copy(struct_expr)
    @assert struct_expr.head == :struct

    return add_field_to_struct(struct_expr, newfield, init_cache, pname, init_cache_param, mod)
end


text:
using BinaryBuilder

name = "p7zip"
version = v"16.02"

# Collection of sources required to build p7zip
sources = [
    ArchiveSource("https://downloads.sourceforge.net/project/p7zip/p7zip/16.02/p7zip_16.02_src_all.tar.bz2",
                  "5eb20ac0e2944f6cb9c2d51dd6c4518941c185347d4089ea89087ffdd6e2341f"),
    FileSource("https://downloads.sourceforge.net/project/sevenzip/7-Zip/19.00/7z1900.exe",
               "759aa04d5b03ebeee13ba01df554e8c962ca339c74f56627c8bed6984bb7ef80"),
    FileSource("https://downloads.sourceforge.net/project/sevenzip/7-Zip/19.00/7z1900-x64.exe",
               "0f5d4dbbe5e55b7aa31b91e5925ed901fdf46a367491d81381846f05ad54c45e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/p7zip_*/

if [[ ${target} == *mingw* ]]; then
    # It's incredibly frustrating to build p7zip on mingw, so instead we just redistribute 7z
    apk add p7zip

    mkdir ${prefix}/bin
    cd ${prefix}/bin
    if [[ ${target} == i686* ]]; then
        7z x -y ${WORKSPACE}/srcdir/7z1900.exe 7z.exe 7z.dll
    else
        7z x -y ${WORKSPACE}/srcdir/7z1900-x64.exe 7z.exe 7z.dll
    fi
    chmod +x 7z.exe 7z.dll
else
    # Build requirements
    apk add nasm yasm

    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/12-CVE-2016-9296.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/13-CVE-2017-17969.patch
    atomic_patch -p4 ${WORKSPACE}/srcdir/patches/15-Enhanced-encryption-strength.patch

    # Convert from target to makefile:
    target_makefile()
    {
        case "${target}" in
            x86_64-linux*)       echo makefile.linux_amd64_asm;;
            i686-linux*)         echo makefile.linux_x86_asm_gcc_4.X;;
            powerpc64le*linux*)  echo makefile.linux_cross_ppc64le;;
            aarch64*linux*)      echo makefile.linux_cross_aarch64;;
            arm-*linux*)         echo makefile.linux_cross_arm;;
            x86_64-*freebsd*)    echo makefile.freebsd6+;;
            x86_64-*darwin*)     echo makefile.macosx_llvm_64bits;;
            aarch64-*darwin*)    echo makefile.macosx_llvm_64bits;;
        esac
    }
    cp $(target_makefile) makefile.machine

    # clang doesn't like this c++11 narrowing, so we disable the error
    if [[ "${target}" == *darwin* ]] || [[ "${target}" == *freebsd* ]]; then
        CXXFLAGS="${CXXFLAGS} -Wno-c++11-narrowing"
    fi

    install_license DOC/License.txt

    make -j${nproc} 7za CC="${CC} ${CFLAGS}" CXX="${CXX} ${CXXFLAGS}"
    mkdir -p ${prefix}/bin
    cp -a bin/7za ${prefix}/bin/7z
fi
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    ExecutableProduct("7z", :p7zip),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
version = v"16.02.1"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

ext: jl

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# ## ロジスティック回帰
# ### 例29
# シグモイド関数を定義しておきます。
# ```julia
# sigmoid(x::Matrix, β::Vector) = 1 / (1 + exp(-x*β))
# ```
using Joe, Plots
using Plots.PlotMeasures # hide
Plots.reset_defaults() # hide
default( # hide
    titlefont  = font("JuliaMono", default(:titlefontsize),  ),  # hide
    guidefont  = font("JuliaMono",  default(:guidefontsize),  ), # hide
    tickfont   = font("JuliaMono", default(:tickfontsize),   ), # hide
    legendfont = font("JuliaMono",  default(:legendfontsize), ), # hide
    left_margin = 30px, # hide
    bottom_margin = 30px # hide
) # hide
p29 = plot(xlims=(-10,10), xlabel="x",ylabel="P(Y=1|x)",
            legend=:topleft, title="ロジスティック曲線")
x29 = insert_ones(-10:0.1:10)
for i in [0, 0.2, 0.5, 1, 2, 10]
    plot!(p29, x29[:,2], Joe.sigmoid(x29,[0,i]),label="$i")
end

savefig(p29,joinpath(@OUTPUT,"fig2-1.svg")) # hide
# \fig{fig2-1}

# ## ニュートン・ラフソン法
# ### 例31
# $x^2-1=0$を4を初期値として数値計算で解く問題。10回でどの程度近づくか。
# せっかくなので、 微分をあらわに書き下すのではなくて、自動微分パッケージを使ってみることにします。
# 見た目が$x \leftarrow x - \dfrac{f(x)}{f'(x)}$という公式のまんまですね。
using Zygote
f(x) = x^2 -1
let # hide
x=4
for i in 1:10
    x -= f(x)/f'(x)
    println(x)
end
end # hide

# ### 例31
# 2変数関数が2つある場合。 Zygote.jlにjacobianがなかったので、ForwardDiff.jlを使ってみます。

using ForwardDiff

f(x,y) = [x^2+y^2-1, x+y]
let # hide
z = [3,4]
for i in 1:10
    z -= ForwardDiff.jacobian(x->f(x[1], x[2]),z) \ f(z[1],z[2])
    println(z)
end
end # hide

# ### 例33
# 最尤推定によりロジスティック回帰を行う問題。
# $\nabla l(\beta_0, \beta) = 0$を求める問題を解いてみる。
using Joe, Random, LinearAlgebra
using Plots
N=1000; p=2; Random.seed!(1)
X = insert_ones(randn(N,p))
β = randn(p+1)
prob = @. 1/(1 + exp($*(X,β)))
threshold=0.5
y = ifelse.(rand(N) .> prob,1,-1) # ここまでデータ生成
@show y β;
# せっかくなのでデータを可視化してみます。
data1 = X[y .== 1, 2:3]
data2 = X[y .== -1, 2:3]
p33=scatter(data1[:,1],data1[:,2],ylims=(-5,5),marker=:auto, label="y=1")
scatter!(p33,data2[:,1],data2[:,2],marker=:auto,label="y=-1")
plot!(p33,x-> -β[1]/β[3] - β[2]/β[3]*x, label="β")
savefig(p33,joinpath(@OUTPUT,"fig2-2.svg")) # hide
# \fig{fig2-2}
γ = randn(2+1) #初期値
γ2 = copy(γ) #別解用にコピー
@show γ
W(v::Vector) = @.(v/(1+ v)^2) |> diagm
t = true
for i in 1:10
    global γ
    s=X*γ
    v = @. exp(-y*s)
    u = @. y*v/(1+v)
    γ += ((X'*W(v)*X) \ X') * u
    #δ'*δ < 0.001 && (t = false)
    @show γ
end
plot!(x-> -γ[1]/γ[3] - γ[2]/γ[3]*x, label="γ")
# 今度は横着をして、ヘシアンを使って計算してみます。
using Zygote, LinearAlgebra
l(γ,X=X,y=y) = sum(@. log( 1 /(1+exp(*($*(X,γ),-y))))) #対数尤度関数

for i in 1:10
    global γ2
    δ =  Zygote.hessian(l,γ2) \ l'(γ2)
    γ2 -= δ
    @show γ2
end

plot!(p33,x-> -γ2[1]/γ2[3] - γ2[2]/γ2[3]*x, label="γ2")
savefig(p33,joinpath(@OUTPUT,"fig2-2-2.svg")) # hide
# \fig{fig2-2-2}
# \lineskip
# だいたい同じような境界線が得られることが分かりました。
# ### 例34
# データの半分を使ってロジスティック回帰で学習し、残りのデータで検証する練習。
# まずデータを生成します。
using Joe, Random, Distributions, Plots, LinearAlgebra
Random.seed!(123)
n=100
x34 = vcat(randn(n).+1,randn(n).-1) |> insert_ones
y34 = vcat(ones(n),-ones(n));
# 訓練データとテストデータを分けて可視化します。
index = sample(1:2n,2n,replace=false); #ランダムなインデックスを作り
train = index[begin:n];# 学習用とテスト用にインデックスを分ける。
test = index[n+1:end];
X_train = x34[train,:];y_train = y34[train]
X_test = x34[test,:];y_test = y34[test]
p34 = scatter(X_train[:,2],y_train,label="train")
scatter!(p34,X_test[:,2],y_test,label="test")
savefig(p34,joinpath(@OUTPUT,"fig2-2-3.svg")) # hide
# \fig{fig2-2-3}
# \lineskip
# \note{
    # **訓練データとテストデータを分ける方法**
    #
    # 原著のpythonコードでは以下のようになっています。
    # ```python
    # import numpy as np
    # train = np.random.choice(2*n, int(n), replace=false)
    # test = list(set(range(n))-set(train))
    # ```
    # setを使った書き方は便利ですが、juliaでは重複を許さずにランダムに並べたインデックスを二つに分ける
    # 方法を使いました。
    # ```julia
    # using StatsBase
    # index = sample(1:n,n,replace=false);
    # train = index[begin:Int(n/2)];
    # test = index[Int(n/2)+1:end];
    # ```
    # でもよく考えるとsetdiffを使えば十分でした。
    # ```julia
    # train = sample(1:2n,n,replace=false); test = setdiff(1:2n,train)
    # ```
# }
# ループは原著の方法に従いました。
β = [0,0] #初期値
γ = randn(2)
while sum(β-γ)^2 > 0.001
    global β, γ
    local W
    β = γ
    s = X_train*β
    v = @. exp(-s*y_train)
    u = @. y_train*v/(1+v)
    w = @. v/(1+v)^2
    W = diagm(w)
    z = @. s + u/w
    γ = (X_train'*W*X_train)\(X_train'*W*z)
    @show γ;
end
# 次にこのγを使って、テストデータを予測してみます。 正誤をカウントする関数を次のように定義します。
# ```julia
# function table_count(test,pred)
#     @assert length(test) == length(pred)
#     m = unique(test) |> length
#     count = zeros(Int,m,m)
#     for i in 1:length(test)
#         count[test[i],pred[i]] += 1
#     end
#     return count
# end
# ```
using Joe:table_count
y_pred = X_test*γ .|> x -> ifelse(x>0,2,1)
y_answer = y_test .|> x -> ifelse(x>0,2,1)
table = table_count(y_answer,y_pred)
@show table;
正答率 = sum(diag(table)) / sum(table)
# ## 2.3 線形判別と二次判別
# ### 例35
# まずは与えられた平均と分散共分散行列から二次元正規分布のデータを2種類生成します。
using Distributions, Random,Plots, LinearAlgebra, Parameters
μ₁=[2,2];  Σ₁ = [2 0; 0 2]
μ₂=[-3,-3]; Σ₂ = [1 -0.8; -0.8 1]

N = 100;Random.seed!(123)
data1 = rand(MvNormal(μ₁,Σ₁),100) |> transpose
data2 = rand(MvNormal(μ₂,Σ₂),100) |> transpose
p35 = scatter(data1[:,1],data1[:,2])
scatter!(p35, data2[:,1],data2[:,2])

savefig(p35,joinpath(@OUTPUT,"fig2-3.svg")) # hide
# \fig{fig2-3}
# \lineskip
# ここまでデータの生成。
# 生成したデータから平均と分散共分散行列を計算する。
μ̂₁ = mean(data1,dims=1)'; Σ̂₁ = cov(data1,dims=1);
μ̂₂ = mean(data2,dims=1)' ; Σ̂₂ = cov(data2,dims=1);
# \warning{mean関数で行列の列方向の平均をとる時、得られるデータは二次元の横ベクトルになっている。
# 縦ベクトルとして扱いたいので、転置(正確には複素共役）をしている。}
# それぞれの分布の情報を複合型にまとめておきます。
# ```julia
# @with_kw struct QDA
#     μ::Array # 平均
#     Σ::Matrix # 分散共分散行列
#     invΣ::Array = inv(Σ) # 分散共分散行列の逆行列
#     detΣ = det(Σ) # 分散共分散行列の行列式
# end
# ```
# \note{しれっとパラメータを初期化した複合型を使っていますが、ここではParameters.jlとい
# うサードパーティライブラリを使っています。}
# また、多変量正規分布の対数尤度を取得する関数を以下のように定義します。
# ここでは、juliaのfunction-like objectの機能を利用しています。
# また、可変長引数(x...)を利用して、データが3次元以上の場合にも拡張しています。
# キーワード付き引数に事前確率を加えています。
# \note{function-like objectの機能を使うと、複合型のインスタンスを関数のように
# 使うことができ、再利用しやすくなります。}
# ```julia
# function (qda::QDA)(x...,prior=1.0)
#     data = collect(x)
#     @unpack μ,invΣ,detΣ = qda
#     a = -0.5*(data-μ)' * invΣ * (data-μ)
#     a[1] - log(detΣ) - log(prior)
# end
# ```
using Joe: QDA
param1=QDA(μ= μ̂₁,Σ = Σ̂₁);param2=QDA(μ = μ̂₂,Σ = Σ̂₂)

hanbetsu(x,y) = param1(x,y) - param2(x,y)
x35=-5:0.1:5;y35=-5:0.1:5;
scatter(data1[:,1],data1[:,2])
scatter!( data2[:,1],data2[:,2])
p35_2=contour!(x35,y35, hanbetsu.(x35,y35'), title="QDA")
savefig(p35_2,joinpath(@OUTPUT,"fig2-4.svg")) # hide
# \fig{fig2-4}
#
# 線形判別を行う場合は、分散共分散行列が等しいことを仮定します。
# data1とdata2を、それぞれ中心化後に統合したデータセットについて、
# 新たな分散共分散行列を求めることにします。
Σ_L =  vcat(data1 .- μ̂₁' , data2 .- μ̂₂') |> cov
param1_L=QDA(μ= μ̂₁,Σ =Σ_L );param2_L=QDA(μ = μ̂₂,Σ = Σ_L)
hanbetsu_L(x,y) = param1_L(x,y) - param2_L(x,y)
p35_3=contour!(p35,x35,y35, hanbetsu_L.(x35,y35'), title="LDA")
savefig(p35_3,joinpath(@OUTPUT,"fig2-5.svg")) # hide
# \fig{fig2-5}

# ### 例36 (Fisherのあやめ)
# まずRDatasetsからirisのデータセットを読み込みます。
using ScikitLearn, StatsBase, Random
using Joe:QDA
@sk_import datasets: load_iris
iris = load_iris()
x = iris["data"]
y = iris["target"]
# 訓練データとテストデータを分けます。
Random.seed!(123)
n = length(y)
index = sample(1:n,n,replace=false); #ランダムなインデックスを作り
train = index[begin:Int(n/2)];# 学習用とテスト用にインデックスを分ける。
test = index[Int(n/2)+1:end];
X_train= x[train,:]; X_test = x[test,:];
y_train = y[train]; y_test = y[test];
# それぞれのあやめの訓練データの平均と分散共分散行列を求めて保存する。
Params = QDA[]
for i in 1:3
    μ,Σ =  mean_and_cov(X_train[y_train .==i,:])
    push!(Params, QDA(μ=μ', Σ=Σ)) #μは行ベクトルであることに注意
end

# テストデータで検証する。対数尤度を最大化するラベルを選択します。
y_pred = similar(y_test)
for i in 1:length(y_test)
    y_pred[i] = argmax([param(X_test[i,:]...) for param in Params])
end

using Joe:table_count
table_count(y_test,y_pred)

# ### 問29 事前確率が分かっている場合
# あやめの事前確率が[0.5,0.25,0.25]だったときは、対数をとって足せば良いでしょう。
# キーワード付き引数で事前確率を指定できるようにしたので、
# 少し書き換えるだけで対応可能です。今回の場合、予測結果に対して影響は少ないようです。
priors = [0.5,0.25,0.25]
y_pred2 = similar(y_test)
for i in 1:length(y_test)
    y_pred2[i] = argmax([Params[j](X_test[i,:]...;prior=priors[j]) for j in 1:3])
end
table_count(y_test,y_pred2)
# ## K近傍法
# まず、K近傍法による関数を定義します。原著の記法をほぼ踏襲しています。
# タイブレーキングが少し重複した書き方になっているのが残念。
# ```julia
# using LinearAlgebra, StatsBase
# function knn(X_train::Matrix,y_train::Vector,X_test::Vector, k)
#     n = size(X_train)[1]
#     distance = [norm(X_train[i,:]-X_test) for i in 1:n]
#     S = sortperm(distance)[1:k]
#     u = counts(y_train[S], 1:k)
#     u_max = maximum(u)
#     m = findall(c->c==u_max,u)
#     while length(m)!==1 #タイブレーキング
#         k -=1
#         S = S[1:k]
#         u = counts(S, 1:k)
#         u_max = maximum(u)
#         m = findall(c->c==u_max,u)
#     end
#     return m[1]
# end
# ```
# テストデータが複数の場合は多重ディスパッチで対応します。
# ```julia
# function knn(X_train::Matrix,y_train::Vector,X_test::Matrix, k)
#     l = size(X_test)[1]
#     w = Array{Int}(undef,l)
#     for i in 1:l
#         w[i] = knn(X_train,y_train,X_test[i,:], k)
#     end
#     return w
# end
# ```
# ### 例37 Fisherのあやめ returns
# 例36と同じデータを使ってK近傍法でやってみます。
using Joe:knn,table_count
y_pred_knn = knn(X_train,y_train,X_test,3)
table_count(y_test,y_pred_knn)
# この乱数の種の場合だと、二次判別とそれほど精度は変わりません。
# ## ROC曲線
# ### 例38
# x, y という表記はpositive, negativeに
# 添え字を陽性をp(positive), 陰性をn(negative)に変えています。
# また、varはσに変更しています。
using Distributions, Plots
μₚ = 1; μₙ = -1
σₚ=1; σₙ =1
Nₚ = 1000; Nₙ = 10000
positive = rand(Normal(μₚ,σₚ),Nₚ)
negative = rand(Normal(μₙ,σₙ),Nₙ)
θ = exp.(-10:0.1:100);
U = Vector{Float64}(undef,length(θ))
V = Vector{Float64}(undef,length(θ))
for i in 1:length(θ)
    global U,V # hide
    U[i] = sum(@. pdf(Normal(μₚ,σₚ),negative) / pdf(Normal(μₙ,σₙ),negative) > θ[i]) / Nₙ
    V[i] = sum(@. pdf(Normal(μₚ,σₚ),positive) / pdf(Normal(μₙ,σₙ),positive) > θ[i]) / Nₚ
end

AUC = 0
for i in 1:length(θ)-1
    global AUC # hide
    AUC += abs(U[i+1]-U[i])*V[i]
end
p38 = plot(U,V,xlabel="False Positive",ylabel="False Negative",
            title="ROC curve",legend=false,
            ann = (0.5,0.5,"AUC = $(round(AUC,digits=2))"))

savefig(p38,joinpath(@OUTPUT,"fig2-6.svg")) # hide
# \fig{fig2-6}


text:
#TODO: Update internal expand
function _expand_subgraphs(graph::OptiGraph,overlap::Int64)
    subproblem_graphs = []
    boundary_linkedges_list = []
    hypergraph,hyper_map = gethypergraph(graph)

    #NOTE: This could all be calculated simultaneously
    for subgraph in getsubgraphs(mg)

        println("Performing neighborhood expansion...")
        subnodes = all_nodes(subgraph)
        hypernodes = [hyper_map[node] for node in subnodes]

        #Return new hypernodes and hyperedges covered through the expansion.
        overlap_hnodes = Plasmo.neighborhood(hypergraph,hypernodes,overlap)
        overlap_hedges = Plasmo.induced_edges(hypergraph,overlap_hnodes)
        boundary_hedges = Plasmo.incident_edges(hypergraph,overlap_hnodes)

        overlap_nodes = [hyper_map[node] for node in overlap_hnodes]
        overlap_edges = [hyper_map[edge] for edge in overlap_hedges]
        boundary_edges = [hyper_map[edge] for edge in boundary_hedges]

        #Setup subproblem graphs
        subproblem_graph = OptiGraph()
        subproblem_graph.modelnodes = overlap_nodes
        subproblem_graph.linkedges = overlap_edges

        push!(subproblem_graphs,subproblem_graph)
        push!(boundary_linkedges_list,boundary_edges)
    end

    return subproblem_graphs,boundary_linkedges_list
end

#Links can be formulated as either constraints or penalties
function _assign_links(subgraphs,subgraph_boundary_edges,input_primal_links,input_dual_links)
    subgraph_primal_links = []
    subgraph_dual_links = []

    for (i,edge_set) in enumerate(subgraph_boundary_edges)
        primal_links = LinkConstraintRef[]
        dual_links = LinkConstraintRef[]

        for edge in edge_set
            linkrefs = edge.linkrefs
            for linkref in linkrefs
                if linkref in input_primal_links
                    push!(primal_links,linkref)
                elseif linkref in input_dual_links
                    push!(dual_links,linkref)
                else
                    #use attached node to assign link
                    target_node = constraint_object(linkref).attached_node
                    #target_node = edge.nodes[end]
                    if !(target_node in all_nodes(subgraphs[i]))
                        push!(dual_links,linkref)  #send primal info to target node, receive dual info
                    else
                        push!(primal_links,linkref)   #receive primal info at target node, send back dual info
                    end
                end
            end
        end
        push!(subgraph_primal_links,primal_links)
        push!(subgraph_dual_links,dual_links)
    end
    return subgraph_primal_links,subgraph_dual_links
end

#Find boundary edges of expanded subgraphs
function _find_boundaries(optigraph::OptiGraph,subgraphs::Vector{OptiGraph})
    boundary_linkedges_list = []
    for subgraph in subgraphs
        subnodes = all_nodes(subgraph)
        boundary_edges = Plasmo.incident_edges(optigraph,subnodes)
        push!(boundary_linkedges_list,boundary_edges)
    end
    return boundary_linkedges_list
end


text:
# BigBed Record
# =============

# type Record is defined in reader.jl.

"""
    BigWig.Record()

Create an unfilled bigWig record.
"""
function Record()
    return Record(0, 0, NaN32)
end

function Base.convert(::Type{GenomicFeatures.Interval}, record::Record)
    return GenomicFeatures.Interval(
        chrom(record), chromstart(record), chromend(record),
        GenomicFeatures.STRAND_BOTH, value(record))
end

function Base.convert(::Type{GenomicFeatures.Interval{Record}}, record::Record)
    return convert(GenomicFeatures.Interval, record)
end

function BioCore.isfilled(record::Record)
    return isdefined(record, :reader)
end

function Base.copy(record::Record)
    copy = Record(record.chromstart, record.chromend, record.value)
    copy.header = record.header
    if isdefined(record, :reader)
        copy.reader = record.reader
    end
    return copy
end

function Base.show(io::IO, record::Record)
    if get(io, :compact, false)
        print(io, chrom(record), ':', chromstart(record), '-', chromend(record), "  ", value(record))
    else
        print(io, summary(record), ':')
        if isfilled(record)
            println(io)
            println(io, "  chromosome: ", chrom(record))
            println(io, "       start: ", chromstart(record))
            println(io, "         end: ", chromend(record))
              print(io, "       value: ", value(record))
        else
            print(io, " <not filled>")
        end
    end
end

# Accessor functions
# ------------------

"""
    chromid(record::Record)::UInt32

Get the chromosome ID of `record`.
"""
function chromid(record::Record)::UInt32
    checkfilled(record)
    return record.header.chromid
end

"""
    chrom(record::Record)::String

Get the chromosome name of `record`.
"""
function chrom(record::Record)::String
    checkfilled(record)
    return record.reader.chrom_names[chromid(record)]
end

function haschrom(record::Record)
    return isfilled(record)
end

function BioCore.seqname(record::Record)
    return chrom(record)
end

function BioCore.hasseqname(record::Record)
    return haschrom(record)
end

"""
    chromstart(record::Record)::Int

Get the start position of `record`.
"""
function chromstart(record::Record)::Int
    checkfilled(record)
    return record.chromstart + 1
end

function haschromstart(record::Record)
    return isfilled(record)
end

function BioCore.leftposition(record::Record)
    return chromstart(record)
end

function BioCore.hasleftposition(record::Record)
    return haschromstart(record)
end

"""
    chromend(record::Record)::Int

Get the end position of `record`.
"""
function chromend(record::Record)::Int
    checkfilled(record)
    return record.chromend % Int
end

function haschromend(record::Record)
    return isfilled(record)
end

function BioCore.rightposition(record::Record)
    return chromend(record)
end

function BioCore.hasrightposition(record::Record)
    return haschromend(record)
end

"""
    value(record::Record)::Float32

Get the value of `record`.
"""
function value(record::Record)::Float32
    checkfilled(record)
    return record.value
end

function hasvalue(record::Record)
    return isfilled(record)
end

function checkfilled(record::Record)
    if !isfilled(record)
        throw(ArgumentError("bigWig record is not filled"))
    end
end


text:
# Html components utilities
function input_field(id, label; value="")
    html_div([dcc_input(id=id, type="text", value=value),html_span("  "*label),html_break(1)])
end

function html_break(n)
    html_div([html_br() for i = 1:n])
end

function row_div(row_data)
    html_div(className="row",[html_div(className="six columns", data) for data in row_data])
end

function load_markdown(filename)
    s = open(filename) do file
        read(file, String)
    end
    html_div(dcc_markdown(s))
end

export input_field, html_break, row_div, load_markdown


text:
using NBodyIPs, JuLIP, Test, IPFitting, DataFrames
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB
using NBodyIPs: BondLengthDesc, BondAngleDesc
using NBodyIPs.Polys: Chebyshev
const Lsq = IPFitting.Lsq
const Err = IPFitting.Errors

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.2*r0, 30, calc)

TRANSFORM = "exp( - 2*(r/$r0 - 1))"

# exp( - 2*(0.9*r0/r0-1))

rcut2 = cutoff(calc)*1.6
CUTOFF2 = "(:cos, $(rcut2-1), $rcut2)"
D2 = BondLengthDesc(TRANSFORM, CUTOFF2)

rcut3 = cutoff(calc)*1.6
CUTOFF3 = "(:cos, $(rcut3-1), $rcut3)"
D3 = BondLengthDesc(TRANSFORM, CUTOFF3)

##
err_erms = Float64[]
err_frms = Float64[]
degrees = [4, 6, 8, 10] # , 10, 12]
rr = range(0.9*r0, stop=cutoff(calc), length=200)
for deg3 in degrees
   # B = [B1; gen_basis(2, D2, deg2); gen_basis(3, D3, deg3)]
   # B = [nbpolys(2, D2, 8, Chebyshev()); nbpolys(3, D3, deg3, Chebyshev())]
   B = nbpolys(3, D3, deg3, Chebyshev())
   @show length(B)
   db = LsqDB("", B, data)
   IP, errs = Lsq.lsqfit( db,
                         E0 = 0.0,
                         configweights = Dict("rand" => 1.0),
                         dataweights   = Dict("E" => 100.0, "F" => 1.0) )
   # push!(err_erms, Err.relrmse(errs, "rand", "E"))
   # push!(err_frms, Err.relrmse(errs, "rand", "F"))
end

##
df = DataFrame( :degrees => degrees,
                :relrms_E => err_erms,
                :relrms_F => err_frms )
display(df)

(@test minimum(err_erms) < 0.0001) |> println
(@test minimum(err_frms) < 0.01) |> println

# BOND LENGTH
# │ Row │ degrees │ rms_E       │ rms_F      │
# ├─────┼─────────┼─────────────┼────────────┤
# │ 1   │ 4       │ 0.00134475  │ 0.118625   │
# │ 2   │ 6       │ 0.000244722 │ 0.0362398  │
# │ 3   │ 8       │ 0.000115287 │ 0.0179257  │
# │ 4   │ 10      │ 5.28124e-5  │ 0.00888845 │

# BOND ANGLE
# │ Row │ degrees │ relrms_E    │ relrms_F   │
# ├─────┼─────────┼─────────────┼────────────┤
# │ 1   │ 4       │ 0.000523011 │ 0.070237   │
# │ 2   │ 6       │ 0.000445928 │ 0.0477666  │
# │ 3   │ 8       │ 0.000100544 │ 0.016773   │
# │ 4   │ 10      │ 3.7269e-5   │ 0.00568809 │
# │ 5   │ 12      │ 2.00043e-5  │ 0.00419453 │


text:
module TestDuplicates

using Test, DataFrames
const ≅ = isequal

@testset "nonunique" begin
    df = DataFrame(a = [1, 2, 3, 3, 4])
    udf = DataFrame(a = [1, 2, 3, 4])
    @test nonunique(df) == [false, false, false, true, false]
    @test udf == unique(df)
    unique!(df)
    @test df == udf

    pdf = DataFrame(a = CategoricalArray(["a", "a", missing, missing, "b", missing, "a", missing]),
                    b = CategoricalArray(["a", "b", missing, missing, "b", "a", "a", "a"]))
    updf = DataFrame(a = CategoricalArray(["a", "a", missing, "b", missing]),
                    b = CategoricalArray(["a", "b", missing, "b", "a"]))
    @test nonunique(pdf) == [false, false, false, true, false, false, true, true]
    @test nonunique(updf) == falses(5)
    @test updf ≅ unique(pdf)
    unique!(pdf)
    @test pdf ≅ updf
end

end # module


text:
using W2vUtils
using MultivariateStats
using Gadfly

wv = load(ARGS[1], W2vData)

(words, dists) = distance(wv, ARGS[2]; n=15)

vecs = W2vUtils.projection(wv, words)

model = fit(PCA, vecs'; maxoutdim=2)

transvecs = transform(model, vecs')

pca_plot = plot(x=transvecs[1, :], y=transvecs[2, :], label=words, Geom.point, Geom.label)
draw(PDF(ARGS[3], 4inch, 3inch), pca_plot)


text:
using HttpServer
using WebSockets, Requests
import HttpServer.mimetypes

#global Dict to store open connections in
global connections = Dict{Int,WebSocket}()
global usernames   = Dict{Int,String}()

function decodeMessage( msg )
    String(copy(msg))
end

wsh = WebSocketHandler() do req, client
    global connections
    @show connections[client.id] = client
    while true
        msg = read(client)
        msg = decodeMessage(msg)
        if startswith(msg, "setusername:")
            println("SETTING USERNAME: $msg")
            usernames[client.id] = msg[13:end]
        end
        if startswith(msg, "say:")
            println("EMITTING MESSAGE: $msg")
            for (k,v) in connections
                if k != client.id
                    write(v, (usernames[client.id] * ": " * msg[5:end]))
                end
            end
        end
    end
end

dir = Pkg.dir("Blocks","src") # ,"client.html"
println(dir)
# onepage = readstring(Pkg.dir("CodeServer","src","client.html"))
httph = HttpHandler() do req::Request, res::Response

  page = ""
  if ismatch(r"/?username=Hello/",req.resource)
     Response(readstring(dir*"\\html\\client.html"))
  else
    uri = URI(req.resource)
    key = match(r"(?:\.(\w+$))", uri.path)[1]
    mime = mimetypes[ key ]
    # println(dir, "     ", replace(uri.path, "/", "\\"), "     Key: ", mime)
    # for windows
    Response(200, Dict{AbstractString,AbstractString}([("Content-Type",mime)]),readstring(dir* replace(uri.path, "/", "\\") ))
    #Response(200, Dict{AbstractString,AbstractString}([("Content-Type",mime)]),readstring(dir* uri.path ))
  end

end

server = Server(httph, wsh)
println("Server listening on 8000...")
run(server,8000)


text:
function intrinsic_value(contract::VLCallOptionContract, currentPriceValue::Float64)::VLResult

    # initialize -
    payoffValue = 0.0
    profitLossValue = 0.0

    # get data from the contract object -
    sense = contract.sense
    strikePrice = contract.strikePrice
    premiumValue = contract.premimumValue
    number_of_contracts = contract.numberOfContracts
    contract_multiplier = contract.contractMultiplier

    # compute the iv -
    if (sense == :buy)

        # compute the P/L -
        payoffValue = max(0.0, (currentPriceValue - strikePrice))
        profitLossValue = (payoffValue - premiumValue)

    elseif (sense == :sell)

        # compute the P/L -
        payoffValue = min(0.0, -1.0 * (currentPriceValue - strikePrice))
        profitLossValue = (payoffValue + premiumValue)
    end

    # compute the intrinsic value -
    iv = (contract_multiplier) * (number_of_contracts) * payoffValue
    pl = (contract_multiplier) * (number_of_contracts) * profitLossValue

    # make a named tuple -
    named_tuple = (intrinsic_value = iv, pl_value = pl)

    # return -
    return VLResult(named_tuple)
end

function intrinsic_value(contract::VLPutOptionContract, currentPriceValue::Float64)::VLResult

    # initialize -
    payoffValue = 0.0
    profitLossValue = 0.0

    # get data from the contract object -
    sense = contract.sense
    strikePrice = contract.strikePrice
    premiumValue = contract.premimumValue
    number_of_contracts = contract.numberOfContracts
    contract_multiplier = contract.contractMultiplier

    # compute the iv -
    if (sense == :buy)

        # compute the P/L -
        payoffValue = max(0.0, (strikePrice - currentPriceValue))
        profitLossValue = (payoffValue - premiumValue)

    elseif (sense == :sell)

        # compute the P/L -
        payoffValue = min(0.0, -1.0 * (strikePrice - currentPriceValue))
        profitLossValue = (payoffValue + premiumValue)
    end

    # compute the intrinsic value -
    pl = (contract_multiplier) * (number_of_contracts) * profitLossValue
    iv = (contract_multiplier) * (number_of_contracts) * payoffValue

    # make a named tuple -
    named_tuple = (intrinsic_value = iv, pl_value = pl)

    # return -
    return VLResult(named_tuple)
end

function intrinsic_value(equityObject::VLEquityAsset, currentPriceValue::Float64)::VLResult

    # initialize -
    iv = 0.0

    # get data from equityObject -
    purchasePricePerShare = equityObject.purchasePricePerShare
    numberOfShares = equityObject.numberOfShares

    # compute the intrinsic value -
    iv = numberOfShares * (currentPriceValue - purchasePricePerShare)

    # make a named tuple -
    named_tuple = (intrinsic_value = iv, pl_value = iv)

    # return -
    return VLResult(named_tuple)
end

function intrinsic_value(contractSet::Set{VLAbstractAsset}, underlyingPriceValue::Float64)::VLResult

    # initialize -
    tmp_iv_array = Float64[]

    # go through each contract, compute the iv, store -
    for (_, contract) in enumerate(contractSet)

        # compute -
        result = intrinsic_value(contract, underlyingPriceValue)
        if (isa(result.value, Exception) == true)
            return result
        end
        iv_value = result.value.intrinsic_value

        # store -
        push!(tmp_iv_array, iv_value)
    end

    # sum -
    total_iv_value = sum(tmp_iv_array)

    # return -
    return VLResult{Float64}(total_iv_value)
end

function compute_profit_loss_value(contractSet::Set{VLAbstractAsset}, underlyingPriceValue::Float64)::VLResult

    # initialize -
    tmp_pl_array = Float64[]

    # go through each contract, compute the iv, store -
    for (index, contract) in enumerate(contractSet)

        # compute -
        result = intrinsic_value(contract, underlyingPriceValue)
        if (isa(result.value, Exception) == true)
            return result
        end
        pl_value = result.value.pl_value

        # store -
        push!(tmp_pl_array, pl_value)
    end

    # sum -
    total_pl_value = sum(tmp_pl_array)

    # return -
    return VLResult{Float64}(total_pl_value)
end


text:
## -----------  Test setup

# Make a univalued history for testing.
# Key spacing governed by `keySpacing`.
function make_test_history(n :: Integer, dimV;
    offSet :: Float64 = 0.0, keySpacing :: Int = 1)

    h = History(Array{Float64});
    nEl = prod(dimV);
    idx = 1;
    for j = 1 : n
        v = reshape(1 : nEl, dimV...) .+ offSet .+ j;
        push!(h, idx, v);
        idx += keySpacing;
    end
    return h
end

function make_test_mvhistory(n :: Integer)
    h = MVHistory();
    nameV = Symbol.(collect('d' : ('d' + (n-1))));
    offSetV = 1.0 : n;
    for j = 1 : n
        sHist = make_test_history(4, [3,1]; offSet = offSetV[j], keySpacing = j);
        ValueHistoriesLH.add_series!(h, nameV[j], sHist);
    end
    return h
end

test_dir() = joinpath(@__DIR__, "test_files");

mkpath(test_dir());

# ---------------


text:
using SafeTestsets

include("util.jl")

## SLM/
include("SLM/default_slm_kwargs.jl")
include("SLM/util.jl")
include("SLM/property.jl")
include("SLM/solve_slm_system.jl")
include("SLM/slm.jl")
include("SLM/eval.jl")

## solve/
include("solve/init_derivatives.jl")
include("solve/solve_nojump.jl")
include("solve/solve.jl")

## models/

# Li2020
include("models/li2020.jl")
include("models/li2020/nojump_eqm.jl")
include("models/li2020/eqcond.jl")
include("models/li2020/augment_variables.jl")
include("models/li2020/solve.jl")

# const GROUP = get(ENV, "GROUP", "All")
# const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
# const is_TRAVIS = haskey(ENV,"TRAVIS")

# Start Test Script

# @time begin
#     if GROUP == "All" || GROUP == "Interface"
#         # @time @safetestset "Basic SDO Examples" begin include("BasicSDOExamples.jl") end
#     end
# end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/CompilerSupportLibraries_jll.jl

module CompilerSupportLibraries_jll

using Libdl, Base.BinaryPlatforms

const PATH_list = String[]
const LIBPATH_list = String[]

export libgfortran, libstdcxx, libgomp

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libgfortran_handle = C_NULL
libgfortran_path = ""
libstdcxx_handle = C_NULL
libstdcxx_path = ""
libgomp_handle = C_NULL
libgomp_path = ""

if Sys.iswindows()
    if arch(HostPlatform()) == "x86_64"
        const libgcc_s = "libgcc_s_seh-1.dll"
    else
        const libgcc_s = "libgcc_s_sjlj-1.dll"
    end
    const libgfortran = string("libgfortran-", libgfortran_version(HostPlatform()).major, ".dll")
    const libstdcxx = "libstdc++-6.dll"
    const libgomp = "libgomp-1.dll"
elseif Sys.isapple()
    if arch(HostPlatform()) == "aarch64"
        const libgcc_s = "@rpath/libgcc_s.2.dylib"
    else
        const libgcc_s = "@rpath/libgcc_s.1.dylib"
    end
    const libgfortran = string("@rpath/", "libgfortran.", libgfortran_version(HostPlatform()).major, ".dylib")
    const libstdcxx = "@rpath/libstdc++.6.dylib"
    const libgomp = "@rpath/libgomp.1.dylib"
else
    const libgcc_s = "libgcc_s.so.1"
    const libgfortran = string("libgfortran.so.", libgfortran_version(HostPlatform()).major)
    const libstdcxx = "libstdc++.so.6"
    const libgomp = "libgomp.so.1"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libgcc_s_handle = dlopen(libgcc_s)
    global libgcc_s_path = dlpath(libgcc_s_handle)
    global libgfortran_handle = dlopen(libgfortran)
    global libgfortran_path = dlpath(libgfortran_handle)
    global libstdcxx_handle = dlopen(libstdcxx)
    global libstdcxx_path = dlpath(libstdcxx_handle)
    global libgomp_handle = dlopen(libgomp)
    global libgomp_path = dlpath(libgomp_handle)
end

is_available() = true

end  # module CompilerSupportLibraries_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
function apply_analytic(f, grid, t, laser, x₀, y₀, z₀, f_unit)
    mapgrid(grid) do (x,y,z)
        r = SVector{3}(x-x₀, y-y₀, z-z₀)

        f(r, t, laser) .|> f_unit
    end
end

function apply_analytic(f, grid, t, laser;
    f_unit,
    downsample_grid=true,
    x₀=zero(recursive_bottom_eltype(grid)),
    y₀=zero(recursive_bottom_eltype(grid)),
    z₀=zero(recursive_bottom_eltype(grid)))

    r₀ = SVector{3}(x₀, y₀, z₀)
    # TODO: Can this be done more elegantly?
    ret_type = Base._return_type(f, (typeof(r₀),typeof(t),typeof(laser)))
    @debug "Return type of supplied function: $ret_type"
    if ret_type isa Number
        _scalarness = ScalarQuantity()
    else
        _scalarness = VectorQuantity()
    end
    sz = approx_target_size(domain_discretization(typeof(grid)), _scalarness)
    grid = downsample_grid ? downsample(grid, approx_size=sz) : grid
    data = apply_analytic(f, grid, t, laser, x₀, y₀, z₀, f_unit)

    if _scalarness === ScalarQuantity()
        ScalarField(data, grid)
    else
        VectorField(data, grid)
    end
end

function analytic_slice(f, grid, t, laser, dir, slice_location;
    downsample_grid=false,
    full = false,
    f_unit,
    x₀=zero(recursive_bottom_eltype(grid)),
    y₀=zero(recursive_bottom_eltype(grid)),
    z₀=zero(recursive_bottom_eltype(grid))
    )

    grid_slice = selectdim(grid, dir, slice_location)
    slice = apply_analytic(f, grid_slice, t, laser;
        downsample_grid, f_unit, x₀, y₀, z₀)
    if full
        result = getindex(slice; NamedTuple{(dir,)}(slice_location)...)
    else
        result = dropgriddims(dropdims(slice, dims=dir), dims=dir)
    end
    downsample_grid ? downsample(result) : result
end


text:
# This file should be executed manually and locally to update the examples

push!(LOAD_PATH,joinpath(@__DIR__, "../src/"))
using Documenter, Literate, StableApproxEPW

# Generating some documentation files with `Literate.jl`
ex_file = joinpath(@__DIR__, "../example/example.jl")
Literate.markdown(ex_file, dirname(ex_file)*"/../docs/src/"; flavor=Literate.DocumenterFlavor())
Literate.notebook(ex_file, dirname(ex_file))

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
function read_data(file_path::String)
	file_stream::IOStream = open(file_path, "r")
	file_string::String = readstring(file_stream)
	close(file_stream)

	ax::Float64 = read_key(file_string, "ax", Float64)
	ay::Float64 = read_key(file_string, "ay", Float64)
	number_of_frames::Array{Int64, 1} = read_key(file_string, "number_of_frames", Array{Int64, 1})
	deltat::Float64 = read_key(file_string, "deltat", Float64)
	K::Array{Int64, 1} = read_key(file_string, "K", Array{Int64, 1})
	DE::Array{Float64, 1} = read_key(file_string, "DE", Array{Float64, 1})

	return (
		ax,
		ay,
		number_of_frames,
		deltat,
		K,
		DE)
end


text:
using Test
using BenchmarkTools
using Pkg
Pkg.activate(".")

using Revise
using OptFlowUtils


text:
"""
    savegif(filename::String, o::Output; kw...)

Write the output array to a gif.

# Keywords

- `fps`: `Real` frames persecond. Defaults to the `fps` of the output, or `25`.
- `minval`: Minimum value in the grid(s) to normalise for conversion to an RGB pixel.
    `Number` or `Tuple` for multiple grids.
- `maxval`: Maximum value in the grid(s) to normalise for conversion to an RGB pixel.
    `Number` or `Tuple` for multiple grids.
- `font`: `String` name of font to search for. A default will be guessed.
- `text`: `TextCongif()` or `nothing` for no text. Default is `TextCongif(; font=font)`.
- `scheme`: ColorSchemes.jl scheme, `ObjectScheme()` or `Greyscale()`
- `imagegen`: `ImageGenerator` like `Image` or `Layout`. Will be detected automatically
"""
function savegif(filename::String, o::Output, ruleset=Ruleset();
    minval=minval(o), maxval=maxval(o),
    scheme=ObjectScheme(), imagegen=autoimagegen(init(o), scheme),
    font=autofont(), text=TextCongif(), textconfig=text, kw...
)
    im_o = NoDisplayImageOutput(o;
        imageconfig=ImageConfig(init(o);
            minval=minval, maxval=maxval, imagegen=imagegen, textconfig=textconfig
        )
    )
    savegif(filename, im_o, ruleset; kw...)
end
function savegif(filename::String, o::ImageOutput, ruleset=Ruleset();
    imagegen=imagegen(o), fps=fps(o), kw...) length(o) == 1 && @warn "The output has length 1: the saved gif will be a single image"
    ext = extent(o)
    simdata = SimData(ext, ruleset)
    println(tspan(ext))
    images = map(collect(firstindex(o):lastindex(o))) do f
        @set! simdata.currentframe = f
        grid_to_image!(o, simdata)
    end
    array = cat(images..., dims=3)
    FileIO.save(filename, array; fps=fps, kw...)
end

"""
    GifOutput <: ImageOutput

    GifOutput(init; filename, tspan, kw...)

Output that stores the simulation as images and saves a Gif file on completion.

# Arguments:
- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords

- `filename`: File path to save the gif file to.
- `tspan`: `AbstractRange` timespan for the simulation
- `aux`: NamedTuple of arbitrary input data. Use `get(data, Aux(:key), I...)`
    to access from a `Rule` in a type-stable way.
- `mask`: `BitArray` for defining cells that will/will not be run.
- `padval`: padding value for grids with neighborhood rules. The default is `zero(eltype(init))`.
- `font`: `String` font name, used in default `TextConfig`. A default will be guessed.
- `text`: [`TextConfig`](@ref) object or `nothing` for no text.
- `scheme`: ColorSchemes.jl scheme, or `Greyscale()`
- `imagegen`: [`ImageGenerator`](@ref)
- `minval`: minimum value(s) to set colour maximum
- `maxval`: maximum values(s) to set colour minimum
"""
mutable struct GifOutput{T,F<:AbstractVector{T},E,GC,IC,G,N} <: ImageOutput{T,F}
    frames::F
    running::Bool
    extent::E
    graphicconfig::GC
    imageconfig::IC
    gif::G
    filename::N
end
function GifOutput(; frames, running, extent, graphicconfig, imageconfig, filename, kw...)
    gif = _allocgif(imageconfig, extent)
    GifOutput(frames, running, extent, graphicconfig, imageconfig, gif, filename)
end

filename(o::GifOutput) = o.filename
gif(o::GifOutput) = o.gif

showimage(image, o::GifOutput, data::SimData) = gif(o)[:, :, currentframe(data)] .= image

finalisegraphics(o::GifOutput, data::AbstractSimData) = savegif(o)

savegif(o::GifOutput) = savegif(filename(o), o)
function savegif(filename::String, o::GifOutput, fps=fps(o); kw...)
    FileIO.save(filename, gif(o); fps=fps, kw...)
end

_allocgif(i::ImageConfig, e::Extent) = _allocgif(imagegen(i), i::ImageConfig, e::Extent)
function _allocgif(::ImageGenerator, i::ImageConfig, e::Extent)
    zeros(ARGB32, gridsize(e)..., length(tspan(e)))
end
function _allocgif(p::Layout, i::ImageConfig, e::Extent)
    zeros(ARGB32, (gridsize(e) .* size(p.layout))..., length(tspan(e)))
end


text:
function libsizenormalization(x::Array{Float64,2})

   # Library size normalization (see Splatter Package Bioconductor [https://github.com/Oshlack/splatter/blob/382e22e6b22f32f2bb86a4152af369eda19aaba3/R/splat-estimate.R])
   libsize = mapslices(sum, x, dims = 2)
   libmedian = median(libsize)
   normalizedcounts = (x ./ libsize) .* libmedian

   normalizedcounts
end

function counthistograms(x::Array{Float64,2}, particles::Array{Array{Float64,2},1},ncells::Int)
   tmp1 = DataFrame(x);
   tmp1 = DataFrames.melt(tmp1);
   test = repeat(["Ground Truth"], size(tmp1,1));
   tmp1 = hcat(test,tmp1);

   tmp2 = DataFrame(particles[1]);
   tmp2 = DataFrames.melt(tmp2);
   test = repeat(["Generated Data"], size(tmp2,1));
   tmp2 = hcat(test,tmp2);

   data = vcat(tmp1,tmp2)
   data = data[:,[1,3]]

   rename!(data, :x1 => :Setting, :value => :ExpressionCount)

   p1 =  plot(data,   x = "ExpressionCount",
         xgroup="Setting",
         Geom.subplot_grid(Geom.histogram,Coord.cartesian(xmin=0.0, ymin=0.0)),
         Guide.title("Histogram of Expression Counts for n=$ncells Highly Expressed Cells"),
         Guide.xlabel("Expression Count"))
end

function cellwise_meanvarianceplot(x::Array{Float64,2}, particles::Array{Float64,2},ncells::Int;logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_rbm = log.(mapslices(var,particles,dims=2) .+ 1.0)
      means_rbm = log.(mapslices(mean,particles,dims=2) .+ 1.0)

      variance_true = log.(mapslices(var,x,dims=2) .+ 1.0)
      means_true = log.(mapslices(mean,x,dims=2) .+ 1.0)
   else
      # Take means and variances of each row
      variance_rbm = mapslices(var,particles,dims=2)
      means_rbm = mapslices(mean,particles,dims=2)

      variance_true = mapslices(var,x,dims=2)
      means_true = mapslices(mean,x,dims=2)
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["Generated Data"], size(means_rbm,1))
   mean_var_rbm = DataFrame(hcat(test,means_rbm,variance_rbm))

   data = vcat(mean_var_true,mean_var_rbm)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Cell-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log(Mean)"),
      Guide.ylabel("log(Variance)"))
end

function scvi_cellwise_meanvarianceplot(x::Array{Float64,2}, particles::Array{Float64,2},ncells::Int;logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_rbm = log.(mapslices(var,particles,dims=2) .+ 1.0)
      means_rbm = log.(mapslices(mean,particles,dims=2) .+ 1.0)

      variance_true = log.(mapslices(var,x,dims=2) .+ 1.0)
      means_true = log.(mapslices(mean,x,dims=2) .+ 1.0)
   else
      # Take means and variances of each row
      variance_rbm = mapslices(var,particles,dims=2)
      means_rbm = mapslices(mean,particles,dims=2)

      variance_true = mapslices(var,x,dims=2)
      means_true = mapslices(mean,x,dims=2)
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["scVI Generated Data"], size(means_rbm,1))
   mean_var_rbm = DataFrame(hcat(test,means_rbm,variance_rbm))

   data = vcat(mean_var_true,mean_var_rbm)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Cell-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log(Mean)"),
      Guide.ylabel("log(Variance)"))
end

function compare_cellwise_meanvarianceplot(x::Array{Float64,2}, dbmgen::Array{Float64,2}, scvigen::Array{Float64,2};logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_dbm = log.(mapslices(var,dbmgen,dims=2) .+ 1.0)
      means_dbm = log.(mapslices(mean,dbmgen,dims=2) .+ 1.0)

      variance_scvi = log.(mapslices(var,scvigen,dims=2) .+ 1.0)
      means_scvi = log.(mapslices(mean,scvigen,dims=2) .+ 1.0)

      variance_true = log.(mapslices(var,x,dims=2) .+ 1.0)
      means_true = log.(mapslices(mean,x,dims=2) .+ 1.0)
   else
      # Take means and variances of each row
      variance_dbm = mapslices(var,dbmgen,dims=2)
      means_dbm = mapslices(mean,dbmgen,dims=2)

      variance_scvi = mapslices(var,scvigen,dims=2)
      means_scvi = mapslices(mean,scvigen,dims=2)

      variance_true = mapslices(var,x,dims=2)
      means_true = mapslices(mean,x,dims=2)
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["NB-DBM Generated Data"], size(means_dbm,1))
   mean_var_dbm = DataFrame(hcat(test,means_dbm,variance_dbm))

   test = repeat(["scVI Generated Data"], size(means_scvi,1))
   mean_var_scvi = DataFrame(hcat(test,means_scvi,variance_scvi))

   data = vcat(mean_var_true,mean_var_dbm, mean_var_scvi)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Cell-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log(Mean)"),
      Guide.ylabel("log(Variance)"))
end

function genewise_meanvarianceplot(x::Array{Float64,2}, particles::Array{Float64,2};logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_rbm = log10.(mapslices(var,particles,dims=1)' .+ 1.0)
      means_rbm = log10.(mapslices(mean,particles,dims=1)' .+ 1.0)

      variance_true = log10.(mapslices(var,x,dims=1)' .+ 1.0)
      means_true = log10.(mapslices(mean,x,dims=1)' .+ 1.0)
   else
      # Take means and variances of each row
      variance_rbm = mapslices(var,particles,dims=1)'
      means_rbm = mapslices(mean,particles,dims=1)'

      variance_true = mapslices(var,x,dims=1)'
      means_true = mapslices(mean,x,dims=1)'
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["Generated Data"], size(means_rbm,1))
   mean_var_rbm = DataFrame(hcat(test,means_rbm,variance_rbm))

   data = vcat(mean_var_true,mean_var_rbm)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   ticks = [0, 0.5, 1, 2]
   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Gene-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log10(Mean)"),
      Guide.ylabel("log10(Variance)"))
end

function scvi_genewise_meanvarianceplot(x::Array{Float64,2}, particles::Array{Float64,2};logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_rbm = log10.(mapslices(var,particles,dims=1)' .+ 1.0)
      means_rbm = log10.(mapslices(mean,particles,dims=1)' .+ 1.0)

      variance_true = log10.(mapslices(var,x,dims=1)' .+ 1.0)
      means_true = log10.(mapslices(mean,x,dims=1)' .+ 1.0)
   else
      # Take means and variances of each row
      variance_rbm = mapslices(var,particles,dims=1)'
      means_rbm = mapslices(mean,particles,dims=1)'

      variance_true = mapslices(var,x,dims=1)'
      means_true = mapslices(mean,x,dims=1)'
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["scVI Generated Data"], size(means_rbm,1))
   mean_var_rbm = DataFrame(hcat(test,means_rbm,variance_rbm))

   data = vcat(mean_var_true,mean_var_rbm)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   ticks = [0, 0.5, 1, 2]
   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Gene-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log10(Mean)"),
      Guide.ylabel("log10(Variance)"))
end

function compare_genewise_meanvarianceplot(x::Array{Float64,2}, dbmgen::Array{Float64,2}, scvigen::Array{Float64,2}; logarithm::Bool)

   if logarithm
      # Take means and variances of each row
      variance_dbm = log10.(mapslices(var,dbmgen,dims=1)' .+ 1.0)
      means_dbm = log10.(mapslices(mean,dbmgen,dims=1)' .+ 1.0)

      variance_scvi = log10.(mapslices(var,scvigen,dims=1)' .+ 1.0)
      means_scvi = log10.(mapslices(mean,scvigen,dims=1)' .+ 1.0)

      variance_true = log10.(mapslices(var,x,dims=1)' .+ 1.0)
      means_true = log10.(mapslices(mean,x,dims=1)' .+ 1.0)
   else
      # Take means and variances of each row
      variance_dbm = mapslices(var,dbmgen,dims=1)'
      means_dbm = mapslices(mean,dbmgen,dims=1)'

      variance_scvi = mapslices(var,scvigen,dims=1)'
      means_scvi = mapslices(mean,scvigen,dims=1)'

      variance_true = mapslices(var,x,dims=1)'
      means_true = mapslices(mean,x,dims=1)'
   end

   test = repeat(["Ground Truth"], size(means_true,1))
   mean_var_true = DataFrame(hcat(test,means_true,variance_true))

   test = repeat(["NB-DBM Generated Data"], size(means_dbm,1))
   mean_var_dbm = DataFrame(hcat(test,means_dbm,variance_dbm))

   test = repeat(["scVI Generated Data"], size(means_scvi,1))
   mean_var_scvi = DataFrame(hcat(test,means_scvi,variance_scvi))

   data = vcat(mean_var_true,mean_var_dbm,mean_var_scvi)
   rename!(data, :x1 => :Setting, :x2 => :log_Mean , :x3 => :log_Var)

   ticks = [0, 0.5, 1, 2]
   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Mean", y = "log_Var",
      Guide.title("Gene-Wise Mean-Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log10(Mean)"),
      Guide.ylabel("log10(Variance)"))
end

function compare_tsne(tsne_original::Array{Float64,2}, tsne_dbm::Array{Float64,2}, tsne_scvi::Array{Float64,2})

   test = repeat(["Ground Truth"], size(tsne_original,1))
   tsne_original = DataFrame(hcat(test,tsne_original))

   test = repeat(["NB-DBM Generated Data"], size(tsne_dbm,1))
   tsne_dbm = DataFrame(hcat(test,tsne_dbm))

   test = repeat(["scVI Generated Data"], size(tsne_scvi,1))
   tsne_scvi = DataFrame(hcat(test,tsne_scvi))

   data = vcat(tsne_original,tsne_dbm,tsne_scvi)
   rename!(data, :x1 => :Setting, :x2 => :TSNE1 , :x3 => :TSNE2)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "TSNE1", y = "TSNE2",
      Guide.title("t-SNE Plot"),
      Geom.subplot_grid(layer(Geom.point)),
      Guide.xlabel("t-SNE 1"),
      Guide.ylabel("t-SNE 2"))
end

function compare_pca(pca_original::Array{Float64,2}, pca_dbm::Array{Float64,2}, pca_scvi::Array{Float64,2})

   test = repeat(["Ground Truth"], size(pca_original,1))
   pca_original = DataFrame(hcat(test,pca_original))

   test = repeat(["NB-DBM Generated Data"], size(pca_dbm,1))
   pca_dbm = DataFrame(hcat(test,pca_dbm))

   test = repeat(["scVI Generated Data"], size(pca_scvi,1))
   pca_scvi = DataFrame(hcat(test,pca_scvi))

   data = vcat(pca_original,pca_dbm,pca_scvi)
   rename!(data, :x1 => :Setting, :x2 => :PC1 , :x3 => :PC2)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "PC1", y = "PC2",
      Guide.title("Principal Component Analysis"),
      Geom.subplot_grid(layer(Geom.point)),
      Guide.xlabel("PC 1"),
      Guide.ylabel("PC 2"))
end

function meandispersionplot(x::Array{Float64,2}, rbm::NegativeBinomialBernoulliRBM; libsizenormalization::Bool)

   dispersion = 1 ./ rbm.inversedispersion
   if libsizenormalization
      normalizedcounts = libsizenormalization(x)
   else
      normalizedcounts = x
   end

   mean_normalizedcounts = mapslices(mean, normalizedcounts, dims=1)'

   data = DataFrame(hcat(dispersion,mean_normalizedcounts))
   rename!(data, :x1 => :Dispersion, :x2 => :MeanNormalizedCounts)

   p = plot(data, x = "MeanNormalizedCounts", y = "Dispersion", Geom.point,
      Guide.title("Mean-Dispersion Plot"),
      Guide.xlabel("Mean Normalized Counts"),
      Guide.ylabel("Dispersion Estimate"),
      Theme(major_label_font_size=25pt, minor_label_font_size=20pt))
   p
end

function genewise_violinplot(x::Array{Float64,2}, particles::Array{Array{Float64,2},1},gene::Int)
   violinplot = DataFrame(hcat(x[:,gene],particles[1][:,gene]));
   rename!(violinplot, :x1 => :Ground_Truth, :x2 => :Generated_Data);
   violinplot = melt(violinplot);
   p_violin = plot(violinplot, x="variable", y="value", Geom.histogram)
end

function cellwise_violinplot(x::Array{Float64,2}, particles::Array{Float64,2},cell::Int)
   violinplot = DataFrame(hcat(x[cell,:],particles[cell,:]));
   rename!(violinplot, :x1 => :Ground_Truth, :x2 => :Generated_Data);
   violinplot = melt(violinplot);
   p_violin = plot(violinplot, x="variable", y="value", Geom.violin)
end

function genewise_varplot(x::Array{Float64,2}, particles::Array{Array{Float64,2},1})
   truevariance = DataFrame(log.(mapslices(var,x,dims=1)'))
   generatedvariance = DataFrame(log.(mapslices(var,particles[1],dims=1)'))

   test = repeat(["Ground Truth"], size(truevariance,1))
   var_true = DataFrame(hcat(test,truevariance,makeunique=true))

   test = repeat(["Generated Data"], size(generatedvariance,1))
   var_generated = DataFrame(hcat(test,generatedvariance,makeunique=true))

   data = vcat(var_true,var_generated)
   rename!(data, :x1 => :Setting, :x1_1 => :log_Var)

   abline = Geom.abline(color="red", style=:dash)
   p2 = plot(data, xgroup="Setting", x = "log_Var",
      Guide.title("Gene-Wise Variance Plot"),
      Geom.subplot_grid(layer(Geom.point), layer(intercept=[0.0], slope=[1.0], abline)),
      Guide.xlabel("log(Mean)"),
      Guide.ylabel("log(Variance)"),
      Theme(major_label_font_size=25pt, minor_label_font_size=20pt))
end

# https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html
function coefficientofvariation(particles::Array{Array{Float64,2},1}, rbm::NegativeBinomialBernoulliRBM)
   tmp = DataFrame(mapslices(std,particles[1],dims=1)' ./ mu(rbm))
   rename!(tmp, :x1 => :CV)
   tmp
end

function biologicalcoefficientofvariation(particles::Array{Array{Float64,2},1}, rbm::NegativeBinomialBernoulliRBM)
   tmp = DataFrame(mapslices(std,particles[1],dims=1)' ./ mu(rbm))
   rename!(tmp, :x1 => :CV)
   tmp
end

function zeroinflationplot_compare(x::Array{Float64,2}, particles::Array{Array{Float64,2},1}, particles_nozeroinflation::Array{Array{Float64,2},1})

   # genewise mean-expression vs proportion of zero counts for ground truth
   means_true = log.(mapslices(mean,x,dims=1)' .+ 1.0)
   zeroproportion = Array{Float64,2}((mapslices(sum, x .== 0, dims = 1) / size(x,1))')

   plotting_true = DataFrame(hcat(means_true,zeroproportion))
   rename!(plotting_true, :x1 => :Means, :x2 => :Zeroproportions)

   # genewise mean-expression vs proportion of zero counts for generated data
   means_generated = log.(mapslices(mean,particles[1],dims=1)' .+ 1.0)
   zeroproportion_generated = Array{Float64,2}((mapslices(sum, particles[1] .== 0, dims = 1) / size(particles[1],1))')

   plotting_generated = DataFrame(hcat(means_generated,zeroproportion_generated))
   rename!(plotting_generated, :x1 => :Means, :x2 => :Zeroproportions)

   # genewise mean-expression vs proportion of zero counts for generated data
   means_generated_noZI = log.(mapslices(mean,particles_nozeroinflation[1],dims=1)' .+ 1.0)
   zeroproportion_generated_noZI = Array{Float64,2}((mapslices(sum, particles_nozeroinflation[1] .== 0, dims = 1) / size(particles_nozeroinflation[1],1))')

   plotting_generated_noZI = DataFrame(hcat(means_generated_noZI,zeroproportion_generated_noZI))
   rename!(plotting_generated_noZI, :x1 => :Means, :x2 => :Zeroproportions)

   data_final = vcat(plotting_true, plotting_generated, plotting_generated_noZI)

   test_true = repeat(["Ground Truth"], size(plotting_true,1));
   test_generated = repeat(["ZINB-DBM"], size(plotting_generated,1));
   test_generated_noZI = repeat(["NB-DBM"], size(plotting_generated_noZI,1));

   tmp_final = vcat(test_true,test_generated,test_generated_noZI)

   data_final = hcat(tmp_final,data_final)
   rename!(data_final, :x1 => :Dataset)

   p2 = plot(data_final, x = "Means", y = "Zeroproportions",
      color = "Dataset", Geom.point, Stat.smooth(smoothing=0.75),
      Guide.title("Mean-Zero Relationship"),
      Guide.xlabel("log(Mean)"),
      Guide.ylabel("Proportion of Zeros"),
      Theme(point_size=1.5mm, major_label_font_size=25pt,
      minor_label_font_size=20pt,
      point_label_font_size=20pt,
      key_title_font_size=20pt,
      key_label_font_size=15pt))
   p2
end

function zeroinflationplot(x::Array{Float64,2}, particles::Array{Float64,2})

   # genewise mean-expression vs proportion of zero counts for ground truth
   means_true = log.(mapslices(mean,x,dims=1)' .+ 1.0)
   zeroproportion = Array{Float64,2}((mapslices(sum, x .== 0, dims = 1) / size(x,1))')

   plotting_true = DataFrame(hcat(means_true,zeroproportion))
   rename!(plotting_true, :x1 => :Means, :x2 => :Zeroproportions)

   # genewise mean-expression vs proportion of zero counts for generated data
   means_generated = log.(mapslices(mean,particles,dims=1)' .+ 1.0)
   zeroproportion_generated = Array{Float64,2}((mapslices(sum, particles .== 0, dims = 1) / size(particles,1))')

   plotting_generated = DataFrame(hcat(means_generated,zeroproportion_generated))
   rename!(plotting_generated, :x1 => :Means, :x2 => :Zeroproportions)

   data_final = vcat(plotting_true, plotting_generated)

   test_true = repeat(["Ground Truth"], size(plotting_true,1));
   test_generated = repeat(["NB-DBM"], size(plotting_generated,1));

   tmp_final = vcat(test_true,test_generated)

   data_final = hcat(tmp_final,data_final)
   rename!(data_final, :x1 => :Dataset)

   p2 = plot(data_final, x = "Means", y = "Zeroproportions",
      color = "Dataset", Geom.point,
      Guide.title("Mean-Zero Relationship"),
      Guide.xlabel("log10(Mean)"),
      Guide.ylabel("Proportion of Zeros"),
      Theme(major_label_font_size=25pt, minor_label_font_size=20pt, key_title_font_size=20pt, key_label_font_size=15pt))

   p2
end

function scvi_zeroinflationplot(x::Array{Float64,2}, particles::Array{Float64,2})

   # genewise mean-expression vs proportion of zero counts for ground truth
   means_true = log.(mapslices(mean,x,dims=1)' .+ 1.0)
   zeroproportion = Array{Float64,2}((mapslices(sum, x .== 0, dims = 1) / size(x,1))')

   plotting_true = DataFrame(hcat(means_true,zeroproportion))
   rename!(plotting_true, :x1 => :Means, :x2 => :Zeroproportions)

   # genewise mean-expression vs proportion of zero counts for generated data
   means_generated = log.(mapslices(mean,particles,dims=1)' .+ 1.0)
   zeroproportion_generated = Array{Float64,2}((mapslices(sum, particles .== 0, dims = 1) / size(particles,1))')

   plotting_generated = DataFrame(hcat(means_generated,zeroproportion_generated))
   rename!(plotting_generated, :x1 => :Means, :x2 => :Zeroproportions)

   data_final = vcat(plotting_true, plotting_generated)

   test_true = repeat(["Ground Truth"], size(plotting_true,1));
   test_generated = repeat(["scVI"], size(plotting_generated,1));

   tmp_final = vcat(test_true,test_generated)

   data_final = hcat(tmp_final,data_final)
   rename!(data_final, :x1 => :Dataset)

   p2 = plot(data_final, x = "Means", y = "Zeroproportions",
      color = "Dataset", Geom.point,
      Guide.title("Mean-Zero Relationship"),
      Guide.xlabel("log10(Mean)"),
      Guide.ylabel("Proportion of Zeros"),
      Theme(major_label_font_size=25pt, minor_label_font_size=20pt, key_title_font_size=20pt, key_label_font_size=15pt))

   p2
end

function zerospergene(x::Array{Float64,2}, particles::Array{Array{Float64,2},1})

   Ground_Truth = Array{Float64,2}((mapslices(sum, x .== 0, dims = 1) / size(x,1))')
   Generated = Array{Float64,2}((mapslices(sum, particles[1] .== 0, dims = 1) / size(particles[1],1))')

   combined_zeroproportion = DataFrame(hcat(Ground_Truth, Generated))
   rename!(combined_zeroproportion, :x1 => :Ground_Truth, :x2 => :Generated)

   tmp = DataFrames.melt(combined_zeroproportion)

   p = plot(tmp, x="variable", y="value" , Geom.boxplot,
         Guide.title("Distribution of Zeros per Gene"),
         Guide.xlabel(""),
         Guide.ylabel("Proportion of Zeros per Gene"))
   p
end

function bicluster_heatmap(M::Array{Float64,2}; k = 2, y_names = Array[])

      # cluster and sort cols
      c = Clustering.kmeans(M,k)
      idy =  sortperm(assignments(c))
      M = M[:,idy]

      # cluster and sort rows
      M_t = convert(Array{Float64,2}, M')
      c2 = Clustering.kmeans(M_t,k)
      idx =  sortperm(assignments(c2))
      M_t = M_t[:,idx]

      # generate names of y if non were supplied
      if isempty(y_names)
          y_names = [string("y", i) for i = 1:1:size(M)[2]]
      end

      Plots.heatmap(1:size(M_t)[2], convert(Array{String,1}, y_names[idy]), M_t)
end

function compare_umap(umap_original::DataFrame, umap_dbm::DataFrame, umap_scvi::DataFrame)

   # Generate color palette
   palette = Scale.color_discrete().f(5)

   # Set up data frame for plotting
   groups = vcat(repeat(["Original"], size(umap_original,1)), repeat(["NB-DBM"], size(umap_dbm,1)), repeat(["scVI"], size(umap_scvi,1)))
   plotting_data = hcat(vcat(umap_original, umap_dbm, umap_scvi), groups)
   x = ["UMAP 1", "UMAP 2", "Cluster", "Groups"]
   names!(plotting_data, Symbol.(x))

   # Plot
   p = plot(plotting_data, xgroup="Groups", x="UMAP 1", y="UMAP 2", color="Cluster", Scale.color_discrete_manual(palette..., levels=["1","2","3","4","5"]),
      Geom.subplot_grid(Geom.point), Guide.title(""), Guide.xlabel("UMAP 1")
   )

   p
end


text:
import Base.valtype
import Distributions.logpdf

"""
    ODESolverType

Types inheriting from abstract type `ODESolverType` declare the numerical
schemes used for finding solutions to backward ODEs
"""
abstract type ODESolverType end

"""
    ODEChangePt

Types inheriting from abstract type `ODEChangePt` decide upon which ODE solvers
are to be used and when to ultimately compute the triplet H,Hν,c.
"""
abstract type ODEChangePt end

"""
    NoChangePt <: ODEChangePt

Struct indicating that only solvers for H,Hν,c are to be employed. The field `λ`
indicates the amount of space (in the units of the number of elements of vector)
that needs to be nevertheless reserved for L,M⁺,μ (even though the latter are
not used).
"""
struct NoChangePt <: ODEChangePt
    λ::Int64
    NoChangePt(λ=0) = new(λ)
end

"""
    SimpleChangePt <: ODEChangePt

Struct indicating that both types of solvers are to be used; for H,Hν,c as well
as for L,M⁺,μ. The ODE solvers for L,M⁺,μ are used first on the terminal part
of the interval, where `λ` gives the lenght of this terminal interval. The ODE
solvers for H,Hν,c are used from then on.
"""
struct SimpleChangePt <: ODEChangePt
    λ::Int64
    SimpleChangePt(λ=0) = new(λ)
end

"""
    get_change_pt(changePt::ODEChangePt)

Return the length of terminal interval over which ODE solvers for L,M⁺,μ are
used
"""
get_change_pt(changePt::ODEChangePt) = changePt.λ

"""
    AbstractObsScheme

Types inheriting from abstract type `AbstractObsScheme` define the scheme
according to which a stochastic process has been observed
"""
abstract type AbstractObsScheme end

"""
    PartObs <: AbstractObsScheme

Type acting as a flag for partially observed diffusions
"""
struct PartObs <: AbstractObsScheme end

"""
    FPT <: AbstractObsScheme

Observation scheme in which only first passage times are observed
"""
struct FPT <: AbstractObsScheme end

"""
    ParamUpdateType

Types inheriting from abstract type `ParamUpdateType` define the way in which
parameters are to be updated by the MCMC sampler
"""
abstract type ParamUpdateType end

"""
    ConjugateUpdt <: ParamUpdateType

Type acting as a flag for update from full conditional (conjugate to a prior)
"""
struct ConjugateUpdt <: ParamUpdateType end

"""
    MetropolisHastingsUpdt <: ParamUpdateType

Flag for performing update according to Metropolis Hastings step
"""
struct MetropolisHastingsUpdt <: ParamUpdateType end

"""
    LangevinUpdt <: ParamUpdateType
Flag for performing Metropolis adjusted Langevin updates
"""
struct LangevinUpdt <: ParamUpdateType end

"""
    BlockingSchedule

Types inheriting from abstract type `BlockingSchedule` define the schedule
according to which blocking is done
"""
abstract type BlockingSchedule end

"""
    struct NoBlocking <: BlockingSchedule end

Regular updates with no blocking
"""
struct NoBlocking <: BlockingSchedule end

abstract type ActionType end

struct Verbose <: ActionType end
struct SavePath <: ActionType end
struct Imputation <: ActionType end
struct ParamUpdate <: ActionType end
struct Readjust <: ActionType end

#TODO implement Jeffrey's priors
"""
    ImproperPrior

Flat prior
"""
struct ImproperPrior end
logpdf(::ImproperPrior, θ) = 0.0

valtype(::Val{T}) where T = T

"""
    idx(::Val{T}) where {T}

Return a tuple containing indices of parameters selected by Val{T}

# Examples
```julia-repl
julia> idx(Val((true, false, false, true, true)))
(1, 4, 5)
```
"""
function idx(::Val{T}) where T
    tuple((i for i in 1:length(T) if T[i])...)
end

"""
    moveToProperPlace(ϑ, θ, ::Val{T}) where {T}

Update parameter vector `θ` at indices specified by `Val{T}` with the
values collected in `ϑ`.

# Examples
```julia-repl
julia> moveToProperPlace([10,20,30], [1,2,3,4,5],
                         Val((true, false, true, false, true)))
5-element Array{Float64,1}:
 10.0
  2.0
 20.0
  4.0
 30.0
```
"""
function move_to_proper_place(ϑ, θ, ::Val{v}) where {v}
    θᵒ = copy(θ)
    idxNew = [i for i in 1:length(v) if v[i]]
    θᵒ[idxNew] = ϑ
    θᵒ
end

get_aux_flag(::Any) = nothing

# NOTE quick, but potentially dangerous way to generically handle the new
# cloning functionality, change later!
clone(A, B, C, D) = clone(A, B, C)


text: # Put convenience functions that are not strictly wrappers here


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

# simple array based priority queue, interface similar to
# DataStructures.jl prioirty queue, but simplified
# disabled checking for presence of keys when inserting

# Binary heap indexing
heapleft(i::Integer)   = i << 1        #2i
heapright(i::Integer)  = (i << 1) + 1 #2i + 1
heapparent(i::Integer) = i >> 1      #div(i, 2)

mutable struct SimplePQ{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Vector{Pair{K,V}}

    # Map elements to their index in xs
    index::Vector{Int}

    function SimplePQ{K,V}(itr) where {K,V}
        xs    = Vector{Pair{K,V}}(undef, length(itr))
        index = collect(1:length(itr))
        @inbounds for (i,p) in enumerate(itr)
            xs[i] = p
        end
        pq = new{K,V}(xs, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            heapify(pq, i)
        end

        pq
    end
end

SimplePQ(kv::Vector{T}) where {K,V,T <: Pair{K,V}} = SimplePQ{K,V}(kv)
length(pq::SimplePQ)  = length(pq.xs)
isempty(pq::SimplePQ) = isempty(pq.xs)
peek(pq::SimplePQ)    = (@inbounds return pq.xs[1])
@inline function swap(pq::SimplePQ{K,V}, i::Int, j::Int) where {K,V}
    @inbounds pq.xs[i], pq.xs[j] = pq.xs[j], pq.xs[i]
    @inbounds ki = pq.xs[j].first
    @inbounds kj = pq.xs[i].first
    @inbounds pq.index[ki], pq.index[kj] = pq.index[kj], pq.index[ki]
end

@inbounds function update_pq_at(pq::SimplePQ{K,V}, i::Int) where {K,V}
    pidx = heapparent(i)
    xs   = pq.xs
    if (pidx > 0) && xs[i].second < xs[pidx].second
        swap(pq, i, pidx)
        update_pq_at(pq, pidx)
    else
        ridx = heapright(i)
        lidx = heapleft(i)
        len  = length(pq)
        if ridx <= len
            if (lidx <= len) && (xs[lidx].second < xs[ridx].second)
                ridx = lidx
            end
        else
            if lidx > len
                return
            end
            ridx = lidx
        end

        if xs[ridx].second < xs[i].second
            swap(pq, i, ridx)
            update_pq_at(pq, ridx)
        end
    end
end

@inbounds function heapify(pq::SimplePQ{K,V}, i::Int) where {K,V}
    l   = heapleft(i)
    r   = heapright(i)
    len = length(pq)
    xs  = pq.xs

    smallest = ((l <= len) && (xs[l].second < xs[i].second)) ? l : i

    if (r <= len) && (xs[r].second < xs[smallest].second)
        smallest = r
    end

    if smallest != i
        swap(pq, i, smallest)
        heapify(pq, smallest)
    end
end

function getindex(pq::SimplePQ{K,V}, key) where {K,V}
    @inbounds return pq.xs[pq.index[key]].second
end

# Change the priority of an existing element
function setindex!(pq::SimplePQ{K, V}, value, key) where {K,V}
    @inbounds i        = pq.index[key]
    @inbounds pq.xs[i] = Pair{K,V}(key, value)
    update_pq_at(pq, i)
    value
end

# Unordered iteration through key value pairs in a SimplePQ
start(pq::SimplePQ)   = start(pq.index)
done(pq::SimplePQ, i) = done(pq.index, i)
function next(pq::SimplePQ{K,V}, i) where {K,V}
    (k, idx), i = next(pq.index, i)
    return (pq.xs[idx], i)
end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
@testset "Variable utils" begin
   # TODO fill
end


text:
export makeAxisArray, converttometer, imcenter

converttometer(x::NTuple{D,T}) where {D,T<:Quantity} = ustrip.(uconvert.(u"m",x))
converttometer(x::NTuple{D,T}) where {D,T<:AbstractFloat} = x

imcenter(img::AxisArray) = map(x->(0.5*(last(x)+first(x))),
                               ImageAxes.filter_space_axes(AxisArrays.axes(img), axisvalues(img)))
imcenter(img::ImageMeta) = imcenter(data(img))

function makeAxisArray(array::Array{T,5}, pixelspacing, offset, dt) where T
  N = size(array)
  im = AxisArray(array, Axis{:color}(1:N[1]),
		 Axis{:x}(range(offset[1],step=pixelspacing[1],length=N[2])),
		 Axis{:y}(range(offset[2],step=pixelspacing[2],length=N[3])),
		 Axis{:z}(range(offset[3],step=pixelspacing[3],length=N[4])),
		 Axis{:time}(range(0*unit(dt),step=dt,length=N[5])))
  return im
end

function makeAxisArray(array::Array{T,4}, pixelspacing, offset) where T
  N = size(array)
  im = AxisArray(array, Axis{:color}(1:N[1]),
		 Axis{:x}(range(offset[1],step=pixelspacing[1],length=N[2])),
		 Axis{:y}(range(offset[2],step=pixelspacing[2],length=N[3])),
		 Axis{:z}(range(offset[3],step=pixelspacing[3],length=N[4])))
  return im
end

function makeAxisArray(array::Array{T,3}, pixelspacing, offset) where T
  N = size(array)
  im = AxisArray(array,
		 Axis{:x}(range(offset[1],step=pixelspacing[1],length=N[1])),
		 Axis{:y}(range(offset[2],step=pixelspacing[2],length=N[2])),
		 Axis{:z}(range(offset[3],step=pixelspacing[3],length=N[3])))
  return im
end


text:
"""
	random_initial_guess(lx::Vector,ux::Vector,large_val::Float64=100.0)

Generate random initial guess
"""
function random_initial_guess(lx::Vector,ux::Vector,large_val::Float64=100.0)
	nx = length(lx)
	rand_arr = rand(nx)
	x0 = []
	for i = 1:nx
		if lx[i] == -Inf
			lx_i = -large_val
		else
			lx_i = lx[i]
		end
		if ux[i] == Inf
			ux_i = large_val
		else
			ux_i = ux[i]
		end
		push!(x0, lx_i + rand_arr[i]*(ux_i-lx_i))
	end
	return x0
end

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise