text:
"""
    number_of_states(model)

Return the total number of state variables corresponding to a model.
"""
number_of_states

number_of_states(submodel::Submodel) = val(submodel.nx)

number_of_states(coupling::Coupling) = val(coupling.nx)

number_of_states(model::CoupledModel) = number_of_states(model.coupling)

"""
    number_of_inputs(model)

Return the total number of inputs corresponding to a model.
"""
number_of_inputs

number_of_inputs(submodel::Submodel) = val(submodel.ny)

number_of_inputs(coupling::Coupling) = val(coupling.ny)

number_of_inputs(model::CoupledModel) = number_of_inputs(model.coupling)

"""
    number_of_parameters(model)

Return the total number of parameters corresponding to a model.
"""
number_of_parameters

number_of_parameters(submodel::Submodel) = val(submodel.np)

number_of_parameters(coupling::Coupling) = val(coupling.np)

number_of_parameters(model::CoupledModel) = number_of_parameters(model.coupling)

"""
    number_of_additional_parameters(model)

Return the total number of additional parameters corresponding to a coupled model.
"""
number_of_additional_parameters

number_of_additional_parameters(coupling::Coupling) = val(coupling.npc)

number_of_additional_parameters(model::CoupledModel) = number_of_additional_parameters(model.coupling)

"""
    rate_indices(model::CoupledModel)

Return the rate variable indices for each submodel in `model`
"""
rate_indices(model) = state_indices(model)

"""
    state_indices(model::CoupledModel)

Return the state variable indices for each submodel in `model`
"""
state_indices(model::CoupledModel) = state_indices(model.submodels)

function state_indices(submodels)
    nx = number_of_states.(submodels)
    ix2 = cumsum(nx)
    ix1 = ix2 .- nx .+ 1
    return UnitRange.(ix1, ix2)
end

"""
    input_indices(model::CoupledModel)

Return the input variable indices for each submodel in `model`
"""
input_indices(model::CoupledModel) = input_indices(model.submodels)

function input_indices(submodels)
    ny = number_of_inputs.(submodels)
    iy2 = cumsum(ny)
    iy1 = iy2 .- ny .+ 1
    return UnitRange.(iy1, iy2)
end

"""
    parameter_indices(model::CoupledModel)

Return the parameter indices for each submodel in `model`
"""
parameter_indices(model::CoupledModel) = parameter_indices(model.submodels)

function parameter_indices(submodels)
    np = number_of_parameters.(submodels)
    ip2 = cumsum(np)
    ip1 = ip2 .- np .+ 1
    return UnitRange.(ip1, ip2)
end

"""
    additional_parameter_indices(model::CoupledModel)

Return the indices corresponding to additional parameters introduced by the model coupling
"""
function additional_parameter_indices(model::CoupledModel)
    submodels = model.submodels
    ip1 = sum(number_of_parameters.(submodels)) + 1
    ip2 = number_of_parameters(model)
    return UnitRange(ip1, ip2)
end

"""
    set_rates(submodel::Submodel; kwargs...)

Set the rate vector values for `submodel` using the values in `kwargs`.
"""
set_rates(submodel; kwargs...) = set_states(submodel; kwargs...)

"""
    set_states(submodel::Submodel; kwargs...)

Set the state vector values for `submodel` using the values in `kwargs`.
"""
function set_states(submodel::Submodel; kwargs...)
    x = zeros(number_of_states(submodel))
    set_states!(x, submodel; kwargs...)
    return x
end

"""
    set_inputs(submodel::Submodel; kwargs...)

Set the input vector values for `submodel` using the values in `kwargs`.
"""
function set_inputs(submodel::Submodel; kwargs...)
    y = zeros(number_of_inputs(submodel))
    set_inputs!(y, submodel; kwargs...)
    return y
end

"""
    set_parameters(submodel::Submodel; kwargs...)

Set the parameter vector values for `submodel` using the values in `kwargs`.
"""
function set_parameters(submodel::Submodel; kwargs...)
    p = zeros(number_of_parameters(submodel))
    set_parameters!(p, submodel; kwargs...)
    return p
end

"""
    set_parameters(coupling::Coupling; kwargs...)

Set the additional parameter vector values for `coupling` using the values in `kwargs`.
"""
function set_parameters(coupling::Coupling; kwargs...)
    p = zeros(number_of_additional_parameters(coupling))
    set_parameters!(p, coupling; kwargs...)
    return p
end

"""
    set_additional_parameters(model::CoupledModel; kwargs...)

Set the additional parameter vector values for `model` using the values in `kwargs`.
"""
function set_additional_parameters(model::CoupledModel; kwargs...)
    return set_parameters(model.coupling; kwargs...)
end

"""
    set_rates!(dx, submodel::Submodel; kwargs...)

In-place version of [`set_rates`](@ref)
"""
set_rates!(dx, submodel::Submodel; kwargs...) = set_states!(dx, submodel; kwargs...)

"""
    set_states!(x, submodel::Submodel; kwargs...)

In-place version of [`set_states`](@ref)
"""
set_states!(x, submodel::Submodel; kwargs...) = submodel.setstate(x; kwargs...)

"""
    set_inputs!(y, submodel::Submodel; kwargs...)

In-place version of [`set_inputs`](@ref)
"""
set_inputs!(y, submodel::Submodel; kwargs...) = submodel.setinput(y; kwargs...)

"""
    set_parameters!(p, submodel::Submodel; kwargs...)

In-place version of [`set_parameters`](@ref)
"""
set_parameters!(p, submodel::Submodel; kwargs...) = submodel.setparam(p; kwargs...)

set_parameters!(p, coupling::Coupling; kwargs...) = coupling.setparam(p; kwargs...)

"""
    set_states!(x, model::CoupledModel, i; kwargs...)

Set the elements of the state vector `x` corresponding to the `i`th submodel of `model`
to the values in `kwargs`
"""
function set_states!(x, model::CoupledModel, i; kwargs...)

    return set_states!(view(x, state_indices(model)[i]), model.submodels[i]; kwargs...)
end

"""
    set_inputs!(y, model::CoupledModel, i; kwargs...)

Set the elements of the input vector `y` corresponding to the `i`th submodel of `model`
to the values in `kwargs`
"""
function set_inputs!(y, model::CoupledModel, i; kwargs...)

    return set_inputs!(view(y, input_indices(model)[i]), model.submodels[i]; kwargs...)
end

"""
    set_parameters!(p, model::CoupledModel, i; kwargs...)

Set the elements of the parameter vector `p` corresponding to the `i`th submodel of `model`
to the values in `kwargs`
"""
function set_parameters!(p, model::CoupledModel, i; kwargs...)

    return set_parameters!(view(p, parameter_indices(model)[i]), model.submodels[i]; kwargs...)
end

"""
    set_additional_parameters!(p, model::CoupledModel; kwargs...)

Set the elements of the parameter vector `p` corresponding to additional parameters
introduced by the model coupling to the values specified in `kwargs`
"""
function set_additional_parameters!(p, model::CoupledModel; kwargs...)

    return set_parameters!(view(p, additional_parameter_indices(model)), model.coupling; kwargs...)
end

"""
    separate_rates(model, dx)

Separate the rate vector entries in `dx`.
"""
separate_rates(model, dx) = separate_states(model, dx)

"""
    separate_states(model, x)

Separate the state vector entries in `x`.
"""
separate_states

separate_states(submodel::Submodel, x) = submodel.sepstate(x)

function separate_states(model::CoupledModel, x)

    xs = view.(Ref(x), state_indices(model))

    return separate_states.(model.submodels, xs)
end

"""
    separate_inputs(model, y)

Separate the input vector entries in `y`
"""
separate_inputs

separate_inputs(submodel::Submodel, y) = submodel.sepinput(y)

function separate_inputs(model::CoupledModel, y)

    ys = view.(Ref(y), input_indices(model))

    return separate_inputs.(model.submodels, ys)
end

"""
    separate_parameters(model, p)

Separate the parameter vector entries in `p`
"""
separate_parameters

separate_parameters(submodel::Submodel, p) = submodel.sepparam(p)

separate_parameters(coupling::Coupling, p) = coupling.sepparam(p)

function separate_parameters(model::CoupledModel, p)

    pmodels = view.(Ref(p), parameter_indices(model))
    pmodels_sep = separate_parameters.(model.submodels, pmodels)

    padd = view(p, additional_parameter_indices(model))
    padd_sep = separate_parameters(model.coupling, padd)

    return (pmodels_sep..., padd_sep)
end

"""
    get_residual(model, dx, x, y, p, t)

Calculate the residual for the specified model.
"""
function get_residual(model, dx, x, y, p, t)

    return _get_residual(inplaceness(model), model, dx, x, y, p, t)
end

# dispatch to an in-place function
function _get_residual(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    resid = similar(dx, TF)
    get_residual!(resid, model, dx, x, y, p, t)
    return resid
end

# calculate residual for a submodel
_get_residual(::OutOfPlace, submodel::Submodel, dx, x, y, p, t) = submodel.f(dx, x, y, p, t)

# calculate residual for a coupled model
function _get_residual(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    return vcat(get_residual.(model.submodels, dxs, xs, ys, ps, t)...)
end

"""
    get_residual!(resid, model, dx, x, y, p, t)

In-place version of [`get_residual`](@ref)
"""
function get_residual!(resid, model, dx, x, y, p, t)
    return _get_residual!(resid, inplaceness(model), model, dx, x, y, p, t)
end

# dispatch to an out-of-place function
function _get_residual!(resid, ::OutOfPlace, model, dx, x, y, p, t)
    return resid .= get_residual(model, dx, x, y, p, t)
end

# calculate residual for a submodel
function _get_residual!(resid, ::InPlace, submodel::Submodel, dx, x, y, p, t)
   return submodel.f(resid, dx, x, y, p, t)
end

# calculate residual for a coupled model
function _get_residual!(resid, ::InPlace, model::CoupledModel, dx, x, y, p, t)

    ix = state_indices(models)
    iy = input_indices(models)
    ip = parameter_indices(models)

    resids = view(Ref(resid), ix[i])
    dxs = view.(Ref(dx), ix[i])
    xs = view.(Ref(xs), ix[i])
    ys = view.(Ref(ys), iy[i])
    ps = view.(Ref(ps), ip[i])

    get_residual!.(resids, models, dxs, xs, ys, ps, t)

    return resid
end

"""
    get_rate_jacobian(model)
    get_rate_jacobian(model, p)
    get_rate_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the residual function for `model` with respect to the
state rates
"""
function get_rate_jacobian(model, args...)
    return _get_rate_jacobian(inplaceness(model), model, args...)
end

# dispatch to an in-place function
function _get_rate_jacobian(::InPlace, model)
    Nx = number_of_states(model)
    J = Matrix{Float64}(undef, Nx, Nx)
    return get_rate_jacobian!(J, model)
end

# dispatch to an in-place function
function _get_rate_jacobian(::InPlace, model, p)
    TF = promote_type(eltype(p))
    Nx = number_of_states(model)
    J = Matrix{TF}(undef, Nx, Nx)
    return get_rate_jacobian!(J, model, p)
end

# dispatch to an in-place function
function _get_rate_jacobian(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    Nx = number_of_states(model)
    J = Matrix{TF}(undef, Nx, Nx)
    return get_rate_jacobian!(J, model, dx, x, y, p, t)
end

# calculate rate jacobian for a submodel
function _get_rate_jacobian(::OutOfPlace, submodel::Submodel, args...)
    return jacobian(submodel.ratejac, submodel.nx, submodel.nx, args...)
end

# calculate rate jacobian for a coupled model
function _get_rate_jacobian(::OutOfPlace, model::CoupledModel)

    # partial of residual expression wrt state rates
    r_dx = get_rate_partial(model)

    # partial of residual expression wrt to inputs
    r_y = get_input_jacobian(model)

    # partial of input expression wrt state rates
    y_dx = get_coupling_rate_jacobian(model)

    # calculate jacobian (chain rule)
    return r_dx + r_y*y_dx
end

# calculate rate jacobian for a coupled model
function _get_rate_jacobian(::OutOfPlace, model::CoupledModel, p)

    # partial of residual expression wrt state rates
    r_dx = get_rate_partial(model, p)

    # partial of residual expression wrt to inputs
    r_y = get_input_jacobian(model, p)

    # partial of input expression wrt state rates
    y_dx = get_coupling_rate_jacobian(model, p)

    # calculate jacobian (chain rule)
    return r_dx + r_y*y_dx
end

# calculate rate jacobian for a coupled model
function _get_rate_jacobian(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    # partial of residual expression wrt state rates
    r_dx = get_rate_partial(model, dx, x, y, p, t)

    # partial of residual expression wrt to inputs
    r_y = get_input_jacobian(model, dx, x, y, p, t)

    # partial of input expression wrt state rates
    y_dx = get_coupling_rate_jacobian(model, dx, x, p, t)

    return r_dx + r_y*y_dx
end

function get_rate_partial(model::CoupledModel)

    # initialize diagonal blocks
    Jii = get_rate_jacobian.(model.submodels)

    # assemble into block diagonal matrix
    return block_diagonal(Jii)
end

function get_rate_partial(model::CoupledModel, p)

    # get parameters for each submodel
    ip = parameter_indices(model)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_rate_jacobian.(model.submodels, ps)

    # assemble into block diagonal matrix
    return block_diagonal(Jii)
end

function get_rate_partial(model::CoupledModel, dx, x, y, p, t)

    # get variable indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get views into variables for each submodel
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_rate_jacobian.(model.submodels, dxs, xs, ys, ps, t)

    # assemble into block diagonal matrix
    return block_diagonal(Jii)
end

"""
    get_rate_jacobian!(J, model)
    get_rate_jacobian!(J, model, p)
    get_rate_jacobian!(J, model, dx, x, y, p, t)

In-place version of [`get_rate_jacobian`](@ref).
"""
function get_rate_jacobian!(J, model, args...; kwargs...)
    return _get_rate_jacobian!(J, inplaceness(model), model, args...; kwargs...)
end

# dispatch to an out-of-place function
function _get_rate_jacobian!(J, ::OutOfPlace, model, args...; kwargs...)
    return J .= get_rate_jacobian(model, args...)
end

# calculate rate jacobian for a submodel
function _get_rate_jacobian!(J, ::InPlace, submodel::Submodel, args...)
    return jacobian!(J, submodel.ratejac, args...)
end

# calculate rate jacobian for a coupled model
function _get_rate_jacobian!(J, ::InPlace, model::CoupledModel;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dyddx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state and input indices
    ix = state_indices(model)
    iy = input_indices(model)

    # calculate input mass matrix
    get_coupling_rate_jacobian!(dyddx_cache, model)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        get_input_jacobian!(ri_yi, submodels[i])
        for j = 1:length(submodels)
            ri_dxj = view(J, ix[i], ix[j])
            yi_dxj = view(dyddx_cache, iy[i], ix[j])
            if i == j
                get_rate_jacobian!(ri_dxj, submodels[i])
                mul!(ri_dxj, ri_yi, yi_dxj, 1, 1)
            else
                mul!(ri_dxj, ri_yi, yi_dxj)
            end
        end
    end

    return J
end

# calculate rate jacobian for a coupled model
function _get_rate_jacobian!(J, ::InPlace, model::CoupledModel, p;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dyddx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # calculate input rate jacobian
    get_coupling_rate_jacobian!(dyddx_cache, model, p)

    # calculate rate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        pi = view(p, ip[i])
        get_input_jacobian!(ri_yi, submodels[i], pi)
        for j = 1:length(submodels)
            ri_dxj = view(J, ix[i], ix[j])
            yi_dxj = view(dyddx_cache, iy[i], ix[j])
            if i == j
                get_rate_jacobian!(ri_dxj, submodels[i], pi)
                mul!(ri_dxj, ri_yi, yi_dxj, 1, 1)
            else
                mul!(ri_dxj, ri_yi, yi_dxj)
            end
        end
    end

    return J
end

# calculate rate jacobian for a combination of models
function _get_rate_jacobian!(J, ::InPlace, model::CoupledModel,
    dx, x, y, p, t;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dyddx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # calculate input rate jacobian
    get_coupling_rate_jacobian!(dyddx_cache, model, dx, x, p, t)

    # calculate rate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        dxi = view(dx, ix[i])
        xi = view(x, ix[i])
        yi = view(y, iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, submodels[i], dxi, xi, yi, pi, t)
        for j = 1:length(submodels)
            ri_dxj = view(J, ix[i], ix[j])
            yi_dxj = view(dyddx_cache, iy[i], ix[j])
            if i == j
                get_rate_jacobian!(ri_dxj, submodels[i], dxi, xi, yi, pi, t)
                mul!(ri_dxj, ri_yi, yi_dxj, 1, 1)
            else
                mul!(ri_dxj, ri_yi, yi_dxj)
            end
        end
    end

    return J
end

"""
    get_state_jacobian(model)
    get_state_jacobian(model, p)
    get_state_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the residual function for `model` with respect to the
state variables
"""
function get_state_jacobian(model, args...)
    return _get_state_jacobian(inplaceness(model), model, args...)
end

# dispatch to an in-place function
function _get_state_jacobian(::InPlace, model)
    Nx = number_of_states(model)
    J = Matrix{Float64}(undef, Nx, Nx)
    return get_state_jacobian!(J, model)
end

# dispatch to an in-place function
function _get_state_jacobian(::InPlace, model, p)
    TF = eltype(p)
    Nx = number_of_states(model)
    J = Matrix{TF}(undef, Nx, Nx)
    return get_state_jacobian!(J, model, p)
end

# dispatch to an in-place function
function _get_state_jacobian(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    Nx = number_of_states(model)
    J = Matrix{TF}(undef, Nx, Nx)
    return get_state_jacobian!(J, model, dx, x, y, p, t)
end

# calculate state jacobian for a submodel
function _get_state_jacobian(::OutOfPlace, submodel::Submodel, args...)
    return jacobian(submodel.statejac, submodel.nx, submodel.nx, args...)
end

# calculate state jacobian for a coupled model
function _get_state_jacobian(::OutOfPlace, model::CoupledModel)

    # initialize mass matrix
    r_x = get_state_partial(model)

    # calculate input jacobian
    r_y = get_input_jacobian(model)

    # calculate input jacobian
    y_x = get_coupling_state_jacobian(model)

    return r_x + r_y*y_x
end

# calculate state jacobian for a coupled model
function _get_state_jacobian(::OutOfPlace, model::CoupledModel, p)

    # initialize mass matrix
    r_x = get_state_partial(model, p)

    # calculate input jacobian
    r_y = get_input_jacobian(model, p)

    # calculate input jacobian
    y_x = get_coupling_state_jacobian(model, p)

    return r_x + r_y*y_x
end

# calculate state jacobian for a coupled model
function _get_state_jacobian(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    # initialize mass matrix
    r_x = get_state_partial(model, dx, x, y, p, t)

    # calculate input jacobian
    r_y = get_input_jacobian(model, dx, x, y, p, t)

    # calculate input jacobian
    y_x = get_coupling_state_jacobian(model, dx, x, p, t)

    return r_x + r_y*y_x
end

function get_state_partial(model::CoupledModel)

    submodels = model.submodels

    # initialize diagonal blocks
    Jii = get_state_jacobian.(submodels)

    # assemble block diagonal jacobian
    return block_diagonal(Jii)
end

function get_state_partial(model::CoupledModel, p)

    submodels = model.submodels

    ip = parameter_indices(model)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_state_jacobian.(submodels, ps)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

function get_state_partial(model::CoupledModel, dx, x, y, p, t)

    submodels = model.submodels

    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_state_jacobian.(submodels, dxs, xs, ys, ps, t)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

"""
    get_state_jacobian!(J, model)
    get_state_jacobian!(J, model, p)
    get_state_jacobian!(J, model, dx, x, y, p, t)

In-place version of [`get_state_jacobian`](@ref)
"""
function get_state_jacobian!(J, model, args...; kwargs...)
    return _get_state_jacobian!(J, inplaceness(model), model, args...; kwargs...)
end

# dispatch to an out-of-place function
function _get_state_jacobian!(J, ::OutOfPlace, model, args...; kwargs...)
    return J .= get_state_jacobian(model, args...)
end

# calculate state jacobian for a submodel
function _get_state_jacobian!(J, ::InPlace, submodel::Submodel, args...)
    return jacobian!(J, submodel.statejac, args...)
end

# calculate state jacobian for a coupled model
function _get_state_jacobian!(J, ::InPlace, model::CoupledModel;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state and input indices
    ix = state_indices(model)
    iy = input_indices(model)

    # get input jacobian
    get_coupling_state_jacobian!(dydx_cache, model)

    # calculate jacobian
    for i = 1:length(model.submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        get_input_jacobian!(ri_yi, model.submodels[i])
        for j = 1:length(model.submodels)
            ri_xj = view(J, ix[i], ix[j])
            yi_xj = view(dydx_cache, iy[i], ix[j])
            if i == j
                get_state_jacobian!(ri_xj, model.submodels[i])
                mul!(ri_xj, ri_yi, yi_xj, 1, 1)
            else
                mul!(ri_xj, ri_yi, yi_xj)
            end
        end
    end

    return J
end

# calculate state jacobian for a coupled model
function _get_state_jacobian!(J, ::InPlace, model::CoupledModel, p;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input state jacobian
    get_coupling_state_jacobian!(dydx_cache, model, p)

    # calculate jacobian
    for i = 1:length(model.submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, model.submodels[i], pi)
        for j = 1:length(model.submodels)
            ri_xj = view(J, ix[i], ix[j])
            yi_xj = view(dydx_cache, iy[i], ix[j])
            if i == j
                get_state_jacobian!(ri_xj, model.submodels[i], pi)
                mul!(ri_xj, ri_yi, yi_xj, 1, 1)
            else
                mul!(ri_xj, ri_yi, yi_xj)
            end
        end
    end

    return J
end

# calculate state jacobian for a coupled model
function _get_state_jacobian!(J, ::InPlace, model::CoupledModel, dx, x, y, p, t;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydx_cache = similar(J, number_of_inputs(model), number_of_states(model))
    )

    # get state, input, and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input state jacobian
    get_coupling_state_jacobian!(dydx_cache, model, dx, x, p, t)

    # calculate jacobian
    for i = 1:length(model.submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        dxi = view(dx, ix[i])
        xi = view(x, ix[i])
        yi = view(y, iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, model.submodels[i], dxi, xi, yi, pi, t)
        for j = 1:length(model.submodels)
            ri_xj = view(J, ix[i], ix[j])
            yi_xj = view(dydx_cache, iy[i], ix[j])
            if i == j
                get_state_jacobian!(ri_xj, model.submodels[i], dxi, xi, yi, pi, t)
                mul!(ri_xj, ri_yi, yi_xj, 1, 1)
            else
                mul!(ri_xj, ri_yi, yi_xj)
            end
        end
    end

    return J
end

"""
    get_input_jacobian(model)
    get_input_jacobian(model, p)
    get_input_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the residual function for `model` with respect to the
input variables
"""
function get_input_jacobian(model, args...)
    return _get_input_jacobian(inplaceness(model), model, args...)
end

# dispatch to an in-place function
function _get_input_jacobian(::InPlace, model)
    Nx = number_of_states(model)
    Ny = number_of_inputs(model)
    J = Matrix{Float64}(undef, Nx, Ny)
    return get_input_jacobian!(J, model)
end

# dispatch to an in-place function
function _get_input_jacobian(::InPlace, model, p)
    TF = eltype(p)
    Nx = number_of_states(model)
    Ny = number_of_inputs(model)
    J = Matrix{TF}(undef, Nx, Ny)
    return get_input_jacobian!(J, model, p)
end

# dispatch to an in-place function
function _get_input_jacobian(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    Nx = number_of_states(model)
    Ny = number_of_inputs(model)
    J = Matrix{TF}(undef, Nx, Ny)
    return get_input_jacobian!(J, model, dx, x, y, p, t)
end

# calculate input jacobian for a submodel
function _get_input_jacobian(::OutOfPlace, submodel::Submodel, args...)
    return jacobian(submodel.inputjac, submodel.nx, submodel.ny, args...)
end

# calculate input jacobian for a coupled model
function _get_input_jacobian(::OutOfPlace, model::CoupledModel)

    submodels = model.submodels

    # initialize diagonal blocks
    Jii = get_input_jacobian.(submodels)

    # assemble block diagonal jacobian
    return block_diagonal(Jii)
end

# calculate input jacobian for a coupled model
function _get_input_jacobian(::OutOfPlace, model::CoupledModel, p)

    submodels = model.submodels

    ip = parameter_indices(model)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_input_jacobian.(submodels, ps)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

# calculate input jacobian for a coupled model
function _get_input_jacobian(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    submodels = model.submodels

    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_input_jacobian.(submodels, dxs, xs, ys, ps, t)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

"""
    get_input_jacobian!(J, model)
    get_input_jacobian!(J, model, p)
    get_input_jacobian!(J, model, dx, x, y, p, t)

In-place version of [`get_input_jacobian`](@ref)
"""
function get_input_jacobian!(J, model, args...; kwargs...)
    return _get_input_jacobian!(J, inplaceness(model), model, args...; kwargs...)
end

# dispatch to an out-of-place function
function _get_input_jacobian!(J, ::OutOfPlace, model, args...; kwargs...)
    return J .= get_input_jacobian(model, args...; kwargs...)
end

# calculate input jacobian for a submodel
function _get_input_jacobian!(J, ::InPlace, submodel::Submodel, args...)
    return jacobian!(J, submodel.inputjac, args...)
end

# calculate input jacobian for a coupled model
function _get_input_jacobian!(J, ::InPlace, model::CoupledModel)

    submodels = model.submodels

    # get state and input indices
    ix = state_indices(model)
    iy = input_indices(model)

    # calculate jacobian
    J .= 0
    for i = 1:length(submodels)
        Jii = view(J, ix[i], iy[j])
        get_state_jacobian!(Jii, submodels[i])
    end

    return J
end

# calculate input jacobian for a coupled model
function _get_input_jacobian!(J, ::Constant, ::InPlace, model::CoupledModel, p)

    submodels = model.submodels

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # calculate jacobian
    J .= 0
    for i = 1:length(submodels)
        Jii = view(J, ix[i], iy[j])
        pi = view(p, ip[i])
        get_input_jacobian!(Jii, submodels[i], pi)
    end

    return J
end

# calculate input jacobian for a coupled model
function _get_input_jacobian!(J, ::InPlace, model::CoupledModel, dx, x, y, p, t)

    submodels = model.submodels

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # calculate jacobian
    J .= 0
    for i = 1:length(submodels)
        Jii = view(J, ix[i], iy[j])
        dxi = view(dx, ix[i])
        xi = view(x, ix[i])
        yi = view(y, iy[i])
        pi = view(p, ip[i])
        get_input_jacobian!(Jii, submodels[i], dxi, xi, yi, pi, t)
    end

    return J
end

"""
    get_parameter_jacobian(model)
    get_parameter_jacobian(model, p)
    get_parameter_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the residual function for `model` with respect to the
parameters
"""
function get_parameter_jacobian(model, args...)
    return _get_parameter_jacobian(inplaceness(model), model, args...)
end

# dispatch to an in-place function
function _get_parameter_jacobian(::InPlace, model)
    Nx = number_of_states(model)
    Np = number_of_parameters(model)
    J = Matrix{Float64}(undef, Nx, Np)
    return get_parameter_jacobian!(J, model)
end

# dispatch to an in-place function
function _get_parameter_jacobian(::InPlace, model, p)
    TF = eltype(p)
    Nx = number_of_states(model)
    Np = number_of_parameters(model)
    J = Matrix{TF}(undef, Nx, Np)
    return get_parameter_jacobian!(J, model, p)
end

# dispatch to an in-place function
function _get_parameter_jacobian(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    Nx = number_of_states(model)
    Np = number_of_parameters(model)
    J = Matrix{TF}(undef, Nx, Np)
    return get_parameter_jacobian!(J, model, dx, x, y, p, t)
end

# calculate parameter jacobian for a submodel
function _get_parameter_jacobian(::OutOfPlace, submodel::Submodel, args...)
    return jacobian(submodel.paramjac, submodel.nx, submodel.np, args...)
end

# calculate parameter jacobian for a coupled model
function _get_parameter_jacobian(::OutOfPlace, model::CoupledModel)

    # initialize mass matrix
    r_p = get_parameter_partial(model)

    # calculate input jacobian
    r_y = get_input_jacobian(model)

    # calculate input jacobian
    y_p = get_coupling_parameter_jacobian(model)

    return r_p + r_y*y_p
end

# calculate parameter jacobian for a coupled model
function _get_parameter_jacobian(::OutOfPlace, model::CoupledModel, p)

    # initialize parameter jacobian
    r_p = get_parameter_partial(model, p)

    # calculate input jacobian
    r_y = get_input_jacobian(model, p)

    # calculate coupling function parameter jacobian
    y_p = get_coupling_parameter_jacobian(model, p)

    return r_p + r_y*y_p
end

# calculate parameter jacobian for a coupled model
function _get_parameter_jacobian(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    # initialize parameter jacobian
    r_p = get_parameter_partial(model, dx, x, y, p, t)

    # calculate input jacobian
    r_y = get_input_jacobian(model, dx, x, y, p, t)

    # calculate coupling function parameter jacobian
    y_p = get_coupling_parameter_jacobian(model, dx, x, p, t)

    return r_p + r_y*y_p
end

function get_parameter_partial(model::CoupledModel)

    submodels = model.submodels

    # initialize diagonal blocks
    Jii = get_parameter_jacobian.(submodels)

    # assemble block diagonal jacobian
    return block_diagonal(Jii)
end

function get_parameter_partial(model::CoupledModel, p)

    submodels = model.submodels

    ip = parameter_indices(model)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_parameter_jacobian.(submodels, ps)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

function get_parameter_partial(model::CoupledModel, dx, x, y, p, t)

    submodels = model.submodels

    ix = parameter_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    # initialize diagonal blocks
    Jii = get_parameter_jacobian.(submodels, dxs, xs, ys, ps, t)

    # return block diagonal matrix
    return block_diagonal(Jii)
end

"""
    get_parameter_jacobian!(J, model)
    get_parameter_jacobian!(J, model, p)
    get_parameter_jacobian!(J, model, dx, x, y, p, t)

In-place version of [`get_parameter_jacobian`](@ref)
"""
function get_parameter_jacobian!(J, model, args...; kwargs...)
    return _get_parameter_jacobian!(J, inplaceness(model), model, args...; kwargs...)
end

# dispatch to an out-of-place function
function _get_parameter_jacobian!(J, ::OutOfPlace, model, args...; kwargs...)
    return J .= get_parameter_jacobian(model, args...; kwargs...)
end

# calculate parameter jacobian for a submodel
function _get_parameter_jacobian!(J, ::InPlace, submodel::Submodel, args...)
    return jacobian!(J, submodel.paramjac, args...)
end

# calculate parameter jacobian for a combination of models
function _get_parameter_jacobian!(J, ::InPlace, model::CoupledModel;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydp_cache = similar(J, number_of_inputs(model), number_of_parameters(model)))

    # get state and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input parameter jacobian
    get_coupling_parameter_jacobian!(dydp_cache, model)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        get_input_jacobian!(ri_yi, submodels[i])
        for j = 1:length(submodels)
            ri_pj = view(J, ix[i], ip[j])
            yi_pj = view(dydp_cache, iy[i], ip[j])
            if i == j
                get_parameter_jacobian!(ri_pj, submodels[i])
                mul!(ri_pj, ri_yi, yi_pj, 1, 1)
            else
                mul!(ri_pj, ri_yi, yi_pj)
            end
        end
    end

    return J
end

# calculate parameter jacobian for a combination of models
function _get_parameter_jacobian!(J, ::Constant, ::InPlace,
    model::CoupledModel, dx, x, y, p, t;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydp_cache = similar(J, number_of_inputs(model), number_of_parameters(model))
    )

    # get state, input, and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input jacobian
    get_coupling_parameter_jacobian!(dydp_cache, model, p)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, submodels[i], pi)
        for j = 1:length(submodels)
            ri_pj = view(J, ix[i], ip[j])
            yi_pj = view(dydp_cache, iy[i], ip[j])
            if i == j
                get_parameter_jacobian!(ri_pj, submodels[i], pi)
                mul!(ri_pj, ri_yi, yi_pj, 1, 1)
            else
                mul!(ri_pj, ri_yi, yi_pj)
            end
        end
    end

    return J
end

# calculate parameter jacobian for a combination of model
function _get_parameter_jacobian!(J, ::Union{Linear,Nonlinear}, ::InPlace,
    model::CoupledModel, dx, x, y, p, t;
    drdy_cache = similar(J, number_of_states(model), number_of_inputs(model)),
    dydp_cache = similar(J, number_of_inputs(model), number_of_parameters(model))
    )

    # get state, input, and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input jacobian
    get_coupling_parameter_jacobian!(dydp_cache, model, dx, x, p, t)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        dxi = view(dx, ix[i])
        xi = view(x, ix[i])
        yi = view(y, iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, submodels[i], dxi, xi, yi, pi, t)
        for j = 1:length(submodels)
            ri_pj = view(J, ix[i], ip[j])
            yi_pj = view(dydp_cache, iy[i], ip[j])
            if i == j
                get_parameter_jacobian!(ri_pj, submodels[i], dxi, xi, yi, pi, t)
                mul!(ri_pj, ri_yi, yi_pj, 1, 1)
            else
                mul!(ri_pj, ri_yi, yi_pj)
            end
        end
    end

    return J
end

"""
    get_time_gradient(model)
    get_time_gradient(model, p)
    get_time_gradient(model, dx, x, y, p, t)

Calculate the derivative of the residual function for `model` with respect to time
"""
function get_time_gradient(model, args...)
    return _get_time_gradient(inplaceness(model), model, args...)
end

# dispatch to an in-place function
function _get_time_gradient(::InPlace, model)
    Nx = number_of_states(model)
    dT = Vector{Float64}(undef, Nx)
    return get_time_gradient!(dT, model)
end

# dispatch to an in-place function
function _get_time_gradient(::InPlace, model, p)
    TF = eltype(p)
    Nx = number_of_states(model)
    dT = Vector{TF}(undef, Nx)
    return get_time_gradient!(dT, model, p)
end

# dispatch to an in-place function
function _get_time_gradient(::InPlace, model, dx, x, y, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(y), eltype(p), typeof(t))
    Nx = number_of_states(model)
    dT = Vector{TF}(undef, Nx)
    return get_time_gradient!(dT, model, dx, x, y, p, t)
end

# calculate time gradient for a submodel
function _get_time_gradient(::OutOfPlace, submodel::Submodel, args...)
    return gradient(submodel.tgrad, submodel.nx, args...)
end

# calculate time derivative for a combination of models
function _get_time_gradient(::OutOfPlace, model::CoupledModel)

    # initialize mass matrix
    r_t = get_time_partial(model)

    # calculate input jacobian
    r_y = get_input_jacobian(model)

    # calculate input jacobian
    y_t = get_coupling_time_gradient(model)

    return r_t + r_y*y_t
end

# calculate parameter jacobian for a combination of models
function _get_time_gradient(::OutOfPlace, model::CoupledModel, p)

    # initialize parameter jacobian
    r_t = get_time_partial(model, p)

    # calculate input jacobian
    r_y = get_input_jacobian(model, p)

    # calculate coupling function parameter jacobian
    y_t = get_coupling_time_gradient(model, p)

    return r_t + r_y*y_t
end

# calculate parameter jacobian for a combination of models
function _get_time_gradient(::OutOfPlace, model::CoupledModel, dx, x, y, p, t)

    # initialize parameter jacobian
    r_t = get_time_partial(model, dx, x, y, p, t)

    # calculate input jacobian
    r_y = get_input_jacobian(model, dx, x, y, p, t)

    # calculate coupling function parameter jacobian
    y_t = get_coupling_time_gradient(model, dx, x, p, t)

    return r_t + r_y*y_t
end

function get_time_partial(model::CoupledModel)
    return vcat(get_time_gradient.(model.submodels)...)
end

function get_time_partial(model::CoupledModel, p)

    ip = parameter_indices(model)

    ps = view.(Ref(ps), ip[i])

    return vcat(get_time_gradient.(model.submodels, ps)...)
end

function get_time_partial(model::CoupledModel, dx, x, y, p, t)

    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    dxs = view.(Ref(dx), ix[i])
    xs = view.(Ref(xs), ix[i])
    ys = view.(Ref(ys), iy[i])
    ps = view.(Ref(ps), ip[i])

    return vcat(get_time_gradient.(model.submodels, dxs, xs, ys, ps, t)...)
end

"""
    get_time_gradient!(dT, model)
    get_time_gradient!(dT, model, p)
    get_time_gradient!(dT, model, dx, x, y, p, t)

In-place version of [`get_time_gradient`](@ref)
"""
function get_time_gradient!(dT, model, args...; kwargs...)
    return _get_time_gradient!(dT, inplaceness(model), model; kwargs...)
end

# dispatch to an out-of-place function
function _get_time_gradient!(dT, ::OutOfPlace, model, args...; kwargs...)
    return dT .= get_time_gradient(model, args...; kwargs...)
end

# calculate time gradient for a submodel
function _get_time_gradient!(dT, ::InPlace, submodel::Submodel, args...)
    return gradient!(dT, submodel.tgrad, args...)
end

# calculate parameter jacobian for a coupled model
function _get_time_gradient!(dT, ::InPlace, model::CoupledModel;
    drdy_cache = similar(dT, number_of_states(model), number_of_inputs(model)),
    dydt_cache = similar(dT, number_of_inputs(model))
    )

    # get state and input indices
    ix = state_indices(model)
    iy = input_indices(model)

    # get input jacobian
    get_coupling_time_gradient!(dydt_cache, model)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        get_input_jacobian!(ri_yi, submodels[i])
        for j = 1:length(submodels)
            ri_t = view(dT, ix[i])
            yi_t = view(dydt_cache, iy[i])
            if i == j
                get_parameter_jacobian!(Jij, submodels[i])
                mul!(ri_t, ri_yi, yi_t, 1, 1)
            else
                mul!(ri_t, ri_yi, yi_t)
            end
        end
    end

    return dT
end

# calculate parameter jacobian for a coupled model
function _get_time_gradient!(dT, ::InPlace, model::CoupledModel, p;
    drdy_cache = similar(dT, number_of_states(model), number_of_inputs(model)),
    dydt_cache = similar(dT, number_of_inputs(model))
    )

    # get state, input, and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input jacobian
    get_coupling_time_gradient!(dydt_cache, model, p)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, submodels[i], pi)
        for j = 1:length(submodels)
            ri_t = view(dT, ix[i], ix[j])
            yi_t = view(dydt_cache, iy[i], ix[j])
            if i == j
                get_parameter_jacobian!(ri_t, submodels[i], pi)
                mul!(ri_t, ri_yi, yi_t, 1, 1)
            else
                mul!(ri_t, ri_yi, yi_t)
            end
        end
    end

    return dT
end

# calculate parameter jacobian for a coupled model
function _get_time_gradient!(dT, ::InPlace, model::CoupledModel, dx, x, y, p, t;
    drdy_cache = similar(dT, number_of_states(model), number_of_inputs(model)),
    dydt_cache = similar(dT, number_of_inputs(model))
    )

    # get state, input, and parameter indices
    ix = state_indices(model)
    iy = input_indices(model)
    ip = parameter_indices(model)

    # get input jacobian
    get_coupling_time_gradient!(dydt_cache, model, dx, x, p, t)

    # calculate jacobian
    submodels = model.submodels
    for i = 1:length(submodels)
        ri_yi = view(drdy_cache, ix[i], iy[i])
        dxi = view(dx, ix[i])
        xi = view(x, ix[i])
        yi = view(y, iy[i])
        pi = view(p, ip[i])
        ri_yi = get_input_jacobian!(ri_yi, submodels[i], dxi, xi, yi, pi, t)
        for j = 1:length(submodels)
            ri_t = view(dT, ix[i], ix[j])
            yi_t = view(dydt_cache, iy[i], ix[j])
            if i == j
                get_parameter_jacobian!(ri_t, submodels[i], dxi, xi, yi, pi, t)
                mul!(ri_t, ri_yi, yi_t, 1, 1)
            else
                mul!(ri_t, ri_yi, yi_t)
            end
        end
    end

    return dT
end

"""
    get_coupling_inputs(model::CoupledModel, dx, x, p, t)

Calculate the coupling inputs for a coupled model.
"""
function get_coupling_inputs(model::CoupledModel, dx, x, p, t)
    get_coupling_inputs(model.coupling, dx, x, p, t)
end

function get_coupling_inputs(coupling::Coupling, dx, x, p, t)
    return _get_coupling_inputs(inplaceness(coupling), coupling, dx, x, p, t)
end

# dispatch to an in-place function
function _get_coupling_inputs(::InPlace, coupling, dx, x, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(p), typeof(t))
    Ny = number_of_inputs(coupling)
    y = zeros(TF, Ny)
    get_coupling_inputs!(y, coupling, dx, x, p, t)
    return y
end

# calculate out-of-place coupling inputs
function _get_coupling_inputs(::OutOfPlace, coupling, dx, x, p, t)
    return coupling.g(dx, x, p, t)
end

"""
    get_coupling_inputs!(y, models::CoupledModel, dx, x, p, t)

In-place version of [`get_coupling_inputs`](@ref)
"""
function get_coupling_inputs!(y, model::CoupledModel, dx, x, p, t)
    get_coupling_inputs!(y, model.coupling, dx, x, p, t)
end

function get_coupling_inputs!(y, coupling::Coupling, dx, x, p, t)
    return _get_coupling_inputs!(y, inplaceness(coupling), coupling, dx, x, p, t)
end

# dispatch to an out-of-place function
function _get_coupling_inputs!(y, ::OutOfPlace, coupling, dx, x, p, t)
    return y .= get_coupling_inputs(coupling, dx, x, p, t)
end

# calculate in-place coupling inputs
function _get_coupling_inputs!(y, ::InPlace, coupling, dx, x, p, t)
    return coupling.g(y, dx, x, p, t)
end

"""
    get_coupling_rate_jacobian(model)
    get_coupling_rate_jacobian(model, p)
    get_coupling_rate_jacobian(model, dx, x, p, t)

Calculate the jacobian of the coupling function with respect to the state rates
"""
function get_coupling_rate_jacobian(model::CoupledModel, args...)
    return get_coupling_rate_jacobian(model.coupling, args...)
end

function get_coupling_rate_jacobian(coupling::Coupling, args...)
    return _get_coupling_rate_jacobian(inplaceness(coupling), coupling, args...)
end

# dispatch to an in-place function
function _get_coupling_rate_jacobian(::InPlace, coupling)
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{Float64}(undef, Ny, Nx)
    return get_coupling_rate_jacobian!(J, coupling)
end

# dispatch to an in-place function
function _get_coupling_rate_jacobian(::InPlace, coupling, p)
    TF = eltype(p)
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{TF}(undef, Ny, Nx)
    return get_coupling_rate_jacobian!(J, coupling, p)
end

# dispatch to an in-place function
function _get_coupling_rate_jacobian(::InPlace, coupling, dx, x, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(p), typeof(t))
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{TF}(undef, Ny, Nx)
    return get_coupling_rate_jacobian!(J, coupling, dx, x, p, t)
end

# calculate out-of-place coupling function rate jacobian
function _get_coupling_rate_jacobian(::OutOfPlace, coupling, args...)

    return jacobian(coupling.ratejac, coupling.ny, coupling.nx, args...)
end

"""
    get_coupling_rate_jacobian!(J, model)
    get_coupling_rate_jacobian!(J, model, p)
    get_coupling_rate_jacobian!(J, model, dx, x, p, t)

In-place version of [`get_coupling_rate_jacobian`](@ref).
"""
function get_coupling_rate_jacobian!(J, model::CoupledModel, args...)
    return get_coupling_rate_jacobian!(J, model.coupling, args...)
end

function get_coupling_rate_jacobian!(J, coupling::Coupling, args...)
    return _get_coupling_rate_jacobian!(J, inplaceness(coupling), coupling, args...)
end

# dispatch to an out-of-place function
function _get_coupling_rate_jacobian!(J, ::OutOfPlace, coupling, args...)
    return J .= get_coupling_rate_jacobian(coupling, args...)
end

# calculate in-place coupling function rate jacobian
function _get_coupling_rate_jacobian!(J, ::InPlace, coupling, args...)

    return jacobian!(J, coupling.ratejac, args...)
end

"""
    get_coupling_state_jacobian(model)
    get_coupling_state_jacobian(model, p)
    get_coupling_state_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the coupling function for `model` with respect to the
state variables
"""
function get_coupling_state_jacobian(model::CoupledModel, args...)
    return get_coupling_state_jacobian(model.coupling, args...)
end

function get_coupling_state_jacobian(coupling::Coupling, args...)
    return _get_coupling_state_jacobian(inplaceness(coupling), coupling, args...)
end

# dispatch to an in-place function
function _get_coupling_state_jacobian(::InPlace, coupling)
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{Float64}(undef, Ny, Nx)
    return get_coupling_state_jacobian!(J, coupling)
end

# dispatch to an in-place function
function _get_coupling_state_jacobian(::InPlace, coupling, p)
    TF = eltype(p)
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{TF}(undef, Ny, Nx)
    return get_coupling_state_jacobian!(J, coupling, p)
end

# dispatch to an in-place function
function _get_coupling_state_jacobian(::InPlace, coupling, dx, x, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(p), typeof(t))
    Ny = number_of_inputs(coupling)
    Nx = number_of_states(coupling)
    J = Matrix{TF}(undef, Ny, Nx)
    return get_coupling_state_jacobian!(J, coupling, dx, x, p, t)
end

# calculate out-of-place coupling function state jacobian
function _get_coupling_state_jacobian(::OutOfPlace, coupling, args...)

    return jacobian(coupling.statejac, coupling.ny, coupling.nx, args...)
end

"""
    get_coupling_state_jacobian!(J, model)
    get_coupling_state_jacobian!(J, model, p)
    get_coupling_state_jacobian!(J, model, dx, x, p, t)

In-place version of [`get_coupling_state_jacobian`](@ref)
"""
function get_coupling_state_jacobian!(J, model::CoupledModel, args...)
    return get_coupling_state_jacobian!(J, model.coupling, args...)
end

function get_coupling_state_jacobian!(J, coupling, args...)
    return _get_coupling_state_jacobian!(J, inplaceness(coupling), coupling, args...)
end

# dispatch to an out-of-place function
function _get_coupling_state_jacobian!(J, ::OutOfPlace, coupling, args...)
    return J .= get_coupling_state_jacobian(models, args...)
end

# calculate in-place coupling function state jacobian
function _get_coupling_state_jacobian!(J, ::InPlace, coupling, args...)

    return jacobian!(J, coupling.statejac, args...)
end

"""
    get_coupling_parameter_jacobian(model)
    get_coupling_parameter_jacobian(model, p)
    get_coupling_parameter_jacobian(model, dx, x, y, p, t)

Calculate the jacobian of the coupling function for `model` with respect to the
parameters
"""
function get_coupling_parameter_jacobian(model::CoupledModel, args...)
    return get_coupling_parameter_jacobian(model.coupling, args...)
end

function get_coupling_parameter_jacobian(coupling::Coupling, args...)
    return _get_coupling_parameter_jacobian(inplaceness(coupling), coupling, args...)
end

# dispatch to an in-place function
function _get_coupling_parameter_jacobian(::InPlace, coupling)
    Ny = number_of_inputs(coupling)
    Np = number_of_parameters(coupling)
    J = Matrix{Float64}(undef, Ny, Np)
    return get_coupling_parameter_jacobian!(J, coupling)
end

# dispatch to an in-place function
function _get_coupling_parameter_jacobian(::InPlace, coupling, p)
    TF = eltype(p)
    Ny = number_of_inputs(coupling)
    Np = number_of_parameters(coupling)
    J = Matrix{TF}(undef, Ny, Np)
    return get_coupling_parameter_jacobian!(J, coupling)
end

# dispatch to an in-place function
function _get_coupling_parameter_jacobian(::InPlace, coupling, dx, x, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(p), typeof(t))
    Ny = number_of_inputs(coupling)
    Np = number_of_parameters(coupling)
    J = Matrix{TF}(undef, Ny, Np)
    return get_coupling_parameter_jacobian!(J, coupling, dx, x, p, t)
end

# calculate out-of-place coupling function parameter jacobian
function _get_coupling_parameter_jacobian(::OutOfPlace, coupling, args...)

    return jacobian(coupling.paramjac, coupling.ny, coupling.np, args...)
end

"""
    get_coupling_parameter_jacobian!(J, model)
    get_coupling_parameter_jacobian!(J, model, p)
    get_coupling_parameter_jacobian!(J, model, dx, x, p, t)

In-place version of [`get_coupling_parameter_jacobian`](@ref)
"""
function get_coupling_parameter_jacobian!(J, model::CoupledModel, args...)
    return get_coupling_parameter_jacobian!(J, model.coupling, args...)
end

function get_coupling_parameter_jacobian!(J, coupling::Coupling, args...)
    return _get_coupling_parameter_jacobian!(J, inplaceness(coupling), coupling, args...)
end

# dispatch to an out-of-place function
function _get_coupling_parameter_jacobian!(J, ::OutOfPlace, coupling, args...)
    return J .= get_coupling_parameter_jacobian(coupling, args...)
end

# calculate in-place coupling function parameter jacobian
function _get_coupling_parameter_jacobian!(J, ::InPlace, coupling, args...)

    return jacobian!(J, coupling.paramjac, args...)
end

"""
    get_coupling_time_gradient(model)
    get_coupling_time_gradient(model, p)
    get_coupling_time_gradient(model, dx, x, y, p, t)

Calculate the derivative of the coupling function for `model` with respect to time
"""
function get_coupling_time_gradient(model::CoupledModel, args...)
    return get_coupling_time_gradient(model.coupling, args...)
end

function get_coupling_time_gradient(coupling::Coupling, args...)
    return _get_coupling_time_gradient(inplaceness(coupling), coupling, args...)
end

# dispatch to an in-place function
function _get_coupling_time_gradient(::InPlace, coupling)
    Ny = number_of_inputs(coupling)
    dT = Vector{Float64}(undef, Ny)
    return get_coupling_time_gradient!(dT, coupling)
end

# dispatch to an in-place function
function _get_coupling_time_gradient(::InPlace, coupling, p)
    TF = eltype(p)
    Ny = number_of_inputs(coupling)
    dT = Vector{TF}(undef, Ny)
    return get_coupling_time_gradient!(dT, coupling, p)
end

# dispatch to an in-place function
function _get_coupling_time_gradient(::InPlace, coupling, dx, x, p, t)
    TF = promote_type(eltype(dx), eltype(x), eltype(p), typeof(t))
    Ny = number_of_inputs(coupling)
    dT = Vector{TF}(undef, Ny)
    return get_coupling_time_gradient!(dT, coupling, dx, x, p, t)
end

# calculate out-of-place coupling function time gradient
function _get_coupling_time_gradient(::OutOfPlace, coupling, args...)

    return gradient(coupling.tgrad, coupling.ny, args...)
end

"""
    get_coupling_time_gradient!(dT, model)
    get_coupling_time_gradient!(dT, model, p)
    get_coupling_time_gradient!(dT, model, dx, x, p, t)

In-place version of [`get_coupling_time_gradient`](@ref)
"""
function get_coupling_time_gradient!(dT, model::CoupledModel, args...)
    return get_coupling_time_gradient!(dT, model.coupling, args...)
end

function get_coupling_time_gradient!(dT, coupling::Coupling, args...)
    return _get_coupling_time_gradient!(dT, inplaceness(coupling), coupling, args...)
end

# dispatch to an out-of-place function
function _get_coupling_time_gradient!(dT, ::OutOfPlace, coupling, args...)
    return dT .= get_coupling_time_gradient(coupling, args...)
end

# calculate in-place coupling function time gradient
function _get_coupling_time_gradient!(dT, ::InPlace, coupling, args...)

    return gradient!(dT, coupling.tgrad, args...)
end


text:
# setup path -
_PATH_TO_PROJECT = pwd()
_PATH_TO_SRC = joinpath(_PATH_TO_PROJECT, "src")
_PATH_TO_CONFIG = joinpath(_PATH_TO_PROJECT, "config")
_PATH_TO_OUTPUT = joinpath(_PATH_TO_PROJECT, "output")

# import PKG, activate the project and instaintiate -
import Pkg
Pkg.activate(_PATH_TO_PROJECT)
Pkg.instantiate()

# load packages -
using DifferentialEquations
using Plots
using LaTeXStrings
using TOML
using Logging

# my codes -

# Problem.jl
path_to_problem = joinpath(_PATH_TO_SRC,"Problem.jl")
include("$(path_to_problem)")

# Balances.jl
path_to_balances = joinpath(_PATH_TO_SRC,"Balances.jl")
include("$(path_to_balances)")

# Control.jl
path_to_control = joinpath(_PATH_TO_SRC,"Control.jl")
include("$(path_to_control)")

# Kinetics.jl
path_to_kinetics = joinpath(_PATH_TO_SRC,"Kinetics.jl")
include("$(path_to_kinetics)")

# Solver.jl
path_to_solver = joinpath(_PATH_TO_SRC,"Solver.jl")
include("$(path_to_solver)")

ext: jl


text:
module OpenTelemetryInstrumentationDistributed

using OpenTelemetryAPI
using Distributed

using TOML
const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
const INSTRUMENTATION_INFO = InstrumentationInfo(
    name = string(@__MODULE__),
    version = PKG_VERSION,
    schema_url = "https://oolong.dev/OpenTelemetry.jl/dev/OpenTelemetryInstrumentationDistributed/",
)

const DISTRIBUTED_TRACER = Ref{Tracer}()
const REMOTE_CALLMSG_COUNT = Ref{Counter{UInt}}()

function Distributed.local_remotecall_thunk(f, args::Tuple, kwargs::Base.Pairs)
    () -> begin
        f_name = applicable(nameof, f) ? string(nameof(f)) : "UNKNOWN"
        with_span(
            "Distributed local call [$f_name]",
            DISTRIBUTED_TRACER[];
            kind = SPAN_KIND_INTERNAL,
            attributes = Dict{String,TAttrVal}(
                "Distributed.from_worker_id" => myid(),
                "Distributed.target_worker_id" => myid(),
                "Distributed.remote_call_name" => f_name,
            ),
        ) do
            invokelatest(f, args...; kwargs...)
            REMOTE_CALLMSG_COUNT[](;
                from_worker_id = myid(),
                target_worker_id = myid(),
                remote_call_name = f_name,
                mode = "",  # no mode info in local call
            )
        end
    end
end

struct CallMsgWrapper
    f::Any
    mode::Symbol
    from::Int
    span_ctx::Union{SpanContext,Nothing}
end

function (w::CallMsgWrapper)(args...; kw...)
    name = applicable(nameof, w.f) ? string(nameof(w.f)) : "UNKNOWN"
    with_context(;
        OpenTelemetryAPI.SPAN_KEY_IN_CONTEXT =>
            OpenTelemetryAPI.NonRecordingSpan("", w.span_ctx, nothing),
    ) do
        with_span(
            "Distributed remote call [$name]",
            DISTRIBUTED_TRACER[];
            kind = SPAN_KIND_INTERNAL,
            attributes = Dict{String,TAttrVal}(
                "Distributed.from_worker_id" => w.from,
                "Distributed.target_worker_id" => myid(),
                "Distributed.remote_call_name" => name,
            ),
        ) do
            invokelatest(w.f, args...; kw...)
            REMOTE_CALLMSG_COUNT[](;
                from_worker_id = w.from,
                target_worker_id = myid(),
                remote_call_name = name,
                mode = string(w.mode),
            )
        end
    end
end

struct DummyWrapper
    kwargs::Any
end

Base.iterate(w::DummyWrapper, args...) = iterate(w.kwargs, args...)

function Distributed.CallMsg{M}(f, args::Tuple, kwargs::Base.Pairs) where {M}
    Distributed.CallMsg{M}(
        CallMsgWrapper(f, M, myid(), span_context()),
        args,
        DummyWrapper(kwargs),
    )
end

"""
## Metrics

| Meter Name    | Instrument Name | Instrument Type | Unit | Dimensions                                               | Description             |
|:------------- |:--------------- |:--------------- |:---- |:-------------------------------------------------------- |:----------------------- |
| `Distributed` | `n_remote_call` | `Counter{UInt}` |      | `from_worker_id`, `target_worker_id`, `remote_call_name` | Number of remote calls. |

## Spans

The following attributes are added on span creation:

  - `Distributed.from_worker_id`
  - `Distributed.target_worker_id`, the value is the same with `Distributed.from_worker_id` when the function is executed on local processor.
  - `Distributed.remote_call_name`, which is inferred by `Base.nameof`, fallback to `UNKNOWN` when not implemented.
"""
function init(;
    tracer_provider = global_tracer_provider(),
    meter_provider = global_meter_provider(),
)
    REMOTE_CALLMSG_COUNT[] = Counter{UInt}(
        "n_remote_call",
        Meter(
            "Distributed";
            provider = meter_provider,
            instrumentation_info = INSTRUMENTATION_INFO,
        );
        unit = "",
        description = "Number of remote calls executed.",
    )
    DISTRIBUTED_TRACER[] =
        Tracer(; provider = tracer_provider, instrumentation_info = INSTRUMENTATION_INFO)
end

function __init__()
    init()
end

end # module


text:
"""
    GBDTs

Grammar-Based Decision Tree (GBDT) is a machine learning model that can be used for the interpretable classification and categorization of heterogeneous multivariate time series data.
"""
module GBDTs

export
        GBDTNode,
        GBDT,
        induce_tree,
        partition,
        members_by_bool,
        classify,
        node_members,
        children_id,
        id,
        label,
        gbes_result,
        isleaf,
        children,
        gini,
        gini_loss

using Discretizers
using Reexport
using StatsBase
using TikzGraphs, LightGraphs
@reexport using AbstractTrees
@reexport using ExprRules
@reexport using ExprOptimization

"""
    GBDTNode

Node object of a GBDT.
"""
struct GBDTNode
    id::Int
    label::Int
    gbes_result::Union{Nothing,ExprOptResult}
    children::Vector{GBDTNode}
end
function GBDTNode(id::Int, label::Int)
    GBDTNode(id, label, nothing, GBDTNode[])
end

"""
    GBDT

GBDT model produced by induce_tree.
"""
struct GBDT
    tree::GBDTNode
    catdisc::Union{Nothing,CategoricalDiscretizer}
end

"""
    Counter

Mutable int counter
"""
mutable struct Counter
    i::Int
end

"""
    id(node::GBDTNode)

Returns the node id.
"""
id(node::GBDTNode) = node.id
"""
    label(node::GBDTNode)

Returns the node label.
"""
label(node::GBDTNode) = node.label
"""
    gbes_result(node::GBDTNode)

Returns the result from GBES.
"""
gbes_result(node::GBDTNode) = node.gbes_result
"""
    isleaf(node::GBDTNode)

Returns true if node is a leaf.
"""
isleaf(node::GBDTNode) = isempty(node.children)

#AbstractTrees interface
AbstractTrees.children(node::GBDTNode) = node.children
function AbstractTrees.printnode(io::IO, node::GBDTNode)
    print(io, "$(node.id): label=$(node.label)")
    r = node.gbes_result
    if r != nothing
        print(io, ", loss=$(round(r.loss;digits=2)), $(r.expr)")
    end
end

"""
    ishomogeneous(v::AbstractVector{T}) where T

Returns true if all elements in v are the same.
"""
ishomogeneous(v::AbstractVector{T}) where T = length(unique(v)) == 1

"""
    gini_loss(node::RuleNode, grammar::Grammar, X::AbstractVector{T}, y_truth::AbstractVector{Int},
                     members::AbstractVector{Int}, eval_module::Module;
                     w1::Float64=100.0,
                     w2::Float64=0.1) where T

Default loss function based on gini impurity and number of nodes in the derivation tree.
See Lee et al. "Interpretable categorization of heterogeneous time series data"
"""
function gini_loss(node::RuleNode, grammar::Grammar, X::AbstractVector{T}, y_truth::AbstractVector{Int},
                     members::AbstractVector{Int}, eval_module::Module;
                     w1::Float64=100.0,
                     w2::Float64=0.1) where T
    ex = get_executable(node, grammar)
    y_bool = partition(X, members, ex, eval_module)
    members_true, members_false = members_by_bool(members, y_bool)
    return w1*gini(y_truth[members_true], y_truth[members_false]) + w2*length(node)
end

"""
    induce_tree(grammar::Grammar, typ::Symbol, p::ExprOptAlgorithm, X::AbstractVector{XT},
                        y::AbstractVector{YT}, max_depth::Int, loss::Function=gini_loss,
                        eval_module::Module=Main) where {XT,YT}

Learn a GBDT from labeled data.  Categorical labels are converted to integers.
# Arguments:
- `grammar::Grammar`: grammar
- `typ::Symbol`: start symbol
- `p::ExprOptAlgorithm`: Parameters for ExprOptimization algorithm
- `X::AbstractVector{XT}`: Input data features, e.g., a MultivariateTimeSeries
- `y::AbstractVector{YT}`: Input (class) labels.
- `max_depth::Int`: Maximum depth of GBDT.
- `loss::Function`: Loss function.  See gini_loss() for function signature.
- `eval_module::Module`: Module in which expressions are evaluated.
"""
function induce_tree(grammar::Grammar, typ::Symbol, p::ExprOptAlgorithm, X::AbstractVector{XT},
                        y::AbstractVector{YT}, max_depth::Int, loss::Function=gini_loss,
                        eval_module::Module=Main; kwargs...) where {XT,YT}
    catdisc = CategoricalDiscretizer(y)
    y_truth = encode(catdisc, y)
    induce_tree(grammar, typ, p, X, y_truth, max_depth, loss, eval_module;
                catdisc=(catdisc), kwargs...)
end
"""
    induce_tree(grammar::Grammar, typ::Symbol, p::ExprOptAlgorithm, X::AbstractVector{T},
                        y_truth::AbstractVector{Int}, max_depth::Int, loss::Function=gini_loss,
                        eval_module::Module=Main;
                        catdisc::Union{Nothing,CategoricalDiscretizer}=nothing),
                        verbose::Bool=false) where T

Learn a GBDT from labeled data.
# Arguments:
- `grammar::Grammar`: grammar
- `typ::Symbol`: start symbol
- `p::ExprOptAlgorithm`: Parameters for ExprOptimization algorithm
- `X::AbstractVector{XT}`: Input data features, e.g., a MultivariateTimeSeries
- `y_truth::AbstractVector{Int}`: Input (class) labels.
- `max_depth::Int`: Maximum depth of GBDT.
- `loss::Function`: Loss function.  See gini_loss() for function signature.
- `eval_module::Module`: Module in which expressions are evaluated.
- `catdisc::Union{Nothing,CategoricalDiscretizer}`: Discretizer used for converting the labels.
- `min_members_per_branch::Int`: Minimum number of members for a valid branch.
- `prevent_same_label::Bool`: Prevent split if both branches have the same dominant label
- `verbose::Bool`: Verbose outputs
"""
function induce_tree(grammar::Grammar, typ::Symbol, p::ExprOptAlgorithm, X::AbstractVector{T},
                        y_truth::AbstractVector{Int}, max_depth::Int, loss::Function=gini_loss,
                        eval_module::Module=Main;
                        catdisc::Union{Nothing,CategoricalDiscretizer}=nothing,
                        min_members_per_branch::Int=0,
                        prevent_same_label::Bool=true,
                        verbose::Bool=false) where T
    verbose && println("Starting...")
    @assert length(X) == length(y_truth)
    members = collect(1:length(y_truth))
    node_count = Counter(0)
    node = _split(node_count, grammar, typ, p, X, y_truth, members, max_depth, loss, eval_module,
                 min_members_per_branch=min_members_per_branch,
                 prevent_same_label=prevent_same_label, verbose=verbose)
    return GBDT(node, catdisc)
end
function _split(node_count::Counter, grammar::Grammar, typ::Symbol, p::ExprOptAlgorithm,
                       X::AbstractVector{T}, y_truth::AbstractVector{Int},
                       members::AbstractVector{Int},
                       d::Int, loss::Function, eval_module::Module;
                       min_members_per_branch::Int=0,
                       prevent_same_label::Bool=true,
                       verbose::Bool=false) where T
    id = node_count.i += 1  #assign ids in preorder
    if d == 0 || ishomogeneous(y_truth[members])
        return GBDTNode(id, mode(y_truth[members]))
    end

    #gbes
    gbes_result = optimize(p, grammar, typ, (node,grammar)->loss(node, grammar, X, y_truth,
        members, eval_module); verbose=verbose)
    y_bool = partition(X, members, gbes_result.expr, eval_module)
    members_true, members_false = members_by_bool(members, y_bool)

    #don't create split if split doesn't result in two valid groups
    if length(members_true) <= min_members_per_branch || length(members_false) <= min_members_per_branch
        return GBDTNode(id, mode(y_truth[members]))
    end

    #don't create split if both sides of the split have the same dominant label
    if prevent_same_label && (mode(y_truth[members_true]) == mode(y_truth[members_false]))
        return GBDTNode(id, mode(y_truth[members]))
    end

    child_true = _split(node_count, grammar, typ, p, X, y_truth, members_true, d-1,
        loss, eval_module;
        min_members_per_branch=min_members_per_branch,
        prevent_same_label=prevent_same_label,
        verbose=verbose)
    child_false = _split(node_count, grammar, typ, p, X, y_truth, members_false, d-1,
        loss, eval_module;
        min_members_per_branch=min_members_per_branch,
        prevent_same_label=prevent_same_label,
        verbose=verbose)

    return GBDTNode(id, mode(y_truth[members]), gbes_result, [child_true, child_false])
end

"""
    partition(X::AbstractVector{T}, members::AbstractVector{Int}, expr, eval_module::Module) where T

Returns a Boolean vector of length members containing the results of evaluating expr on each member.  Expressions are evaluated in eval_module.
"""
function partition(X::AbstractVector{T}, members::AbstractVector{Int}, expr, eval_module::Module) where T
    y_bool = Vector{Bool}(undef, length(members))
    for i in eachindex(members)
        @eval eval_module x = $(X[members[i]])
        y_bool[i] = Core.eval(eval_module, expr) #use x in expression
    end
    y_bool
end

"""
    members_by_bool(members::AbstractVector{Int}, y_bool::AbstractVector{Bool})

Returns a tuple containing the results of splitting members by the Boolean values in y_bool.
"""
function members_by_bool(members::AbstractVector{Int}, y_bool::AbstractVector{Bool})
    @assert length(y_bool) == length(members)
    return members[findall(y_bool)], members[findall(!,y_bool)]
end

"""
    gini(v1::AbstractVector{T}, v2::AbstractVector{T}) where T

Returns the gini impurity of v1 and v2 weighted by number of elements.
"""
function gini(v1::AbstractVector{T}, v2::AbstractVector{T}) where T
    N1, N2 = length(v1), length(v2)
    return (N1*gini(v1) + N2*gini(v2)) / (N1+N2)
end
"""
    gini(v::AbstractVector{T}) where T

Returns the Gini impurity of v.  Returns 0.0 if empty.
"""
function gini(v::AbstractVector{T}) where T
    isempty(v) && return 0.0
    return 1.0 - sum(abs2, proportions(v))
end

"""
    Base.length(model::GBDT)

Returns the number of vertices in the GBDT.
"""
Base.length(model::GBDT) = length(model.tree)
"""
    Base.length(root::GBDTNode)

Returns the number of vertices in the tree rooted at root.
"""
function Base.length(root::GBDTNode)
    retval = 1
    for c in root.children
        retval += length(c)
    end
    return retval
end

Base.show(io::IO, model::GBDT) = Base.show(io::IO, model.tree)
Base.show(io::IO, tree::GBDTNode) = print_tree(io, tree)

"""
    Base.display(model::GBDT; edgelabels::Bool=false)

Returns a TikzGraphs plot of the tree.  Turn off edgelabels for cleaner plot.  Left branch is true, right branch is false.
"""
function Base.display(model::GBDT; kwargs...)
    display(model.tree, model.catdisc; kwargs...)
end
"""
    Base.display(root::GBDTNode, catdisc::Uniont{Nothing,CategoricalDiscretizer}=nothing;
                     edgelabels::Bool=false)

Returns a TikzGraphs plot of the tree.  Turn off edgelabels for cleaner plot.  Left branch is true, right branch is false.
If catdisc is supplied, use it to decode the labels.
"""
function Base.display(root::GBDTNode, catdisc::Union{Nothing,CategoricalDiscretizer}=nothing;
                     edgelabels::Bool=false)
    n_nodes = length(root)
    g = DiGraph(n_nodes)
    text_labels, edge_labels = Vector{String}(undef,n_nodes), Dict{Tuple{Int,Int},String}()
    for node in PreOrderDFS(root)
        if node.gbes_result != nothing
            r = node.gbes_result
            text_labels[node.id] = string("$(node.id): $(verbatim(string(r.expr)))")
        else
            label = catdisc != nothing ?  decode(catdisc, node.label) : node.label
            text_labels[node.id] = string("$(node.id): label=$(label)")
        end
        for (i, ch) in enumerate(node.children)
            add_edge!(g, node.id, ch.id)
            edge_labels[(node.id, ch.id)] = i==1 ? "True" : "False"
        end
    end
    if edgelabels
        return TikzGraphs.plot(g, text_labels; edge_labels=edge_labels)
    else
        return TikzGraphs.plot(g, text_labels)
    end
end
#Stay in text mode, escape some latex characters
function verbatim(s::String)
    s = replace(s, "_"=>"\\_")
end

"""
    classify(model::GBDT, X::AbstractVector{T}, members::AbstractVector{Int},
                     eval_module::Module=Main) where T

Predict classification label of each member using GBDT model.  Evaluate expressions in eval_module.
"""
function classify(model::GBDT, X::AbstractVector{T},
    members::AbstractVector{Int}=collect(1:length(X)),
    eval_module::Module=Main) where T

    classify(model.tree, X, members, eval_module; catdisc=model.catdisc)
end
"""
    classify(node::GBDTNode, X::AbstractVector{T}, members::AbstractVector{Int}, eval_module::Module=Main;
                     catdisc::Union{Nothing,CategoricalDiscretizer}=nothing) where T

Predict classification label of each member using GBDT tree.  Evaluate expressions in eval_module.  If catdisc is available, use discretizer to decode labels.
"""
function classify(node::GBDTNode, X::AbstractVector{T},
    members::AbstractVector{Int}=collect(1:length(X)),
    eval_module::Module=Main;
    catdisc::Union{Nothing,CategoricalDiscretizer}=nothing) where T

    y_pred = Vector{Int}(undef,length(members))
    for i in eachindex(members)
        @eval eval_module x=$(X[i])
        y_pred[i] = _classify(node, eval_module)
    end
    if catdisc == nothing
        return y_pred
    else
        return decode(catdisc, y_pred)
    end
end
function _classify(node::GBDTNode, eval_module::Module)
    isleaf(node) && return node.label

    ex = get_expr(node.gbes_result)
    ch =  Core.eval(eval_module, ex) ? node.children[1] : node.children[2]
    return _classify(ch, eval_module)
end

"""
    node_members(model::GBDT, X::AbstractVector{T}, members::AbstractVector{Int},
                      eval_module::Module=Main) where T

Returns the members of each node in the tree.
"""
function node_members(model::GBDT, X::AbstractVector{T}, members::AbstractVector{Int},
                      eval_module::Module=Main) where T
    node_members(model.tree, X, members, eval_module)
end
"""
    node_members(node::gbdtnode, x::abstractvector{t}, members::abstractvector{int},
                      eval_module::module=main) where T

returns the members of each node in the tree.
"""
function node_members(node::GBDTNode, X::AbstractVector{T}, members::AbstractVector{Int},
                      eval_module::Module=Main) where T
    mvec = Vector{Vector{Int}}(undef,length(node))
    _node_members!(mvec, node, X, members, eval_module)
    mvec
end
function _node_members!(mvec::Vector{Vector{Int}}, node::GBDTNode, X::AbstractVector{T},
                        members::AbstractVector{Int}, eval_module::Module) where T
    mvec[node.id] = deepcopy(members)
    isleaf(node) && return

    ex = get_expr(node.gbes_result)
    y_bool = partition(X, members, ex, eval_module)
    members_true, members_false = members_by_bool(members, y_bool)
    _node_members!(mvec, node.children[1], X, members_true, eval_module)
    _node_members!(mvec, node.children[2], X, members_false, eval_module)
end

"""
    Base.getindex(model::GBDT, id::Int)

returns node with id
"""
function Base.getindex(model::GBDT, id::Int)
    for node in PreOrderDFS(model.tree)
        if node.id == id
            return node
        end
    end
    error("node id not found")
end

"""
    children_id(node::GBDTNode)

returns a vector that contains the node ids of the children of node
"""
children_id(node::GBDTNode) = Int[c.id for c in children(node)]

"""
   get_expr(node::GBDTNode)

returns the expression of the node
"""
ExprOptimization.get_expr(node::GBDTNode) = get_expr(node.gbes_result)

end # module


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# rwd2tests_2.jl

# 2nd block
# Test the new functions of ReadWriteDlm2

@testset "2_      dfregex" begin
    @test true  == occursin(dfregex("HH:MM"), "22:00")
    @test true  == occursin(dfregex("H:M"), "1:30")
    @test true  == occursin(dfregex("H:M"), "10:30")
    @test false == occursin(dfregex("H:M"), "10:62")
    @test true  == occursin(dfregex("HHhMM"), "22h00")
    @test false == occursin(dfregex("HHhMM"), "22h")
    @test true  == occursin(dfregex("Hh"), "22h")
    @test true  == occursin(dfregex("HhM\\min"), "22h15min")
    @test true  == occursin(dfregex("yyyy-mm-dd"), "2018-01-23")
    @test true  == occursin(dfregex("yyyy.mm.dd"), "2018.01.23")
    @test false == occursin(dfregex("yyyy-mm-dd"), "2018.01.23")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM:SS"), "2018-01-23T23:52:00")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MM:SS"), "2018-01-23T23:72:00")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMz"), "2018-01-23T23:52+01:00")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMz"), "2018-01-23T23:52+0100")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMz"), "2018-01-23T23:52-10:00")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MM z"), "2018-01-23T23:52+01:00")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MM z"), "2018-01-23T23:52 ")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM z"), "2018-01-23T23:52 10:00")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMZ"), "2018-01-23T23:52UTC")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMZ"), "2018-01-23T23:52CET")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMZ"), "2018-01-23T23:52CEST")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MM Z"), "2018-01-23T23:52CET")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM Z"), "2018-01-23T23:52 CEST")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MM Z"), "2018-01-23T23:52CEST")
    @test false == occursin(dfregex("yyyy-mm-ddTHH:MMZ"), "2018-01-23T23:52CE")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MMZZZZZZ"), "2018-01-23T23:52CEST")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM ZZZZ"), "2018-01-23T23:52 CET")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM zzzzzz"), "2018-01-23T23:52 +11:00")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM zzzz"), "2018-01-23T23:52 -09:00")
    @test true  == occursin(dfregex("yyyy-mm-ddTHH:MM\\Z"), "2018-01-23T23:52Z")
end

@testset "2_new functions" begin
    let data = "2015-01-01;5,1;Text1\n10;19e6;4\n"
        @test isequaldlm(readdlm2(IOBuffer(data)), [Date(2015) 5.1 "Text1";10 190.0e5 4.0], Any)
        @test isequaldlm(readdlm2(IOBuffer(data)), [Date(2015) 5.1 "Text1";10 190.0e5 4.0], Any)
    end

    a = [Date(2015) 5.1 "Text1";10 190.0e5 4.0]
    writedlm2("test.csv", a, decimal='€')
    @test read("test.csv", String) == "2015-01-01;5€1;Text1\n10;1€9e7;4€0\n"
    b = readdlm2("test.csv", rs=(r"(\d)€(\d)", s"\1.\2"))
    rm("test.csv")
    @test b[1] == Date(2015)

    a = DateTime(2017)
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "2017-01-01T00:00:00.0\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test b[1] == DateTime(2017)

    a = Any[Date(2017,1,14) DateTime(2017,2,15,23,40,59)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "2017-01-14;2017-02-15T23:40:59.0\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test isequaldlm(a, b, Any)

    a = Any[Date(2017, 1, 1) DateTime(2017, 2, 15, 23, 0, 0)]
    writedlm2("test.csv", a, dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    @test read("test.csv", String) == "01/2017;15.02.2017/23.h\n"
    writedlm2("test.csv", a, dfs="", dtfs="")
    @test read("test.csv", String) == "2017-01-01;2017-02-15T23:00:00\n"
    writedlm2("test.csv", a, decimal='.')
    @test read("test.csv", String) == "2017-01-01;2017-02-15T23:00:00.0\n"
    writedlm2("test.csv", a, decimal='.', dfs="\\Y\\ear: yyyy", dtfs="\\Y\\ear: yyyy")
    @test read("test.csv", String) == "Year: 2017;Year: 2017\n"
    writedlm2("test.csv", a, dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    @test read("test.csv", String) == "01/2017;15.02.2017/23.h\n"
    b = readdlm2("test.csv", dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    @test isequaldlm(a, b, Any)
    b = readdlm2("test.csv", rs=(), dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    rm("test.csv")
    @test isequaldlm(a, b, Any)

    a = Date(2017,5,1)
    writedlm2("test.csv", a, dfs="dd.mm.yyyy")
    @test read("test.csv", String) == "01.05.2017\n"
    b = readdlm2("test.csv", dfs="dd.mm.yyyy")
    rm("test.csv")
    @test b[1] == a

    a = "ABC"
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "A\nB\nC\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test b[1] == "A" && b[3] == "C"

    D = rand(1:9, 10, 5)
    writedlm2("test.csv", D)
    @test length(read("test.csv", String)) == 100
    b = readdlm2("test.csv")
    rm("test.csv")
    @test D == b

    a = 10.9
    writedlm2("test.csv", a, decimal='€')
    @test read("test.csv", String) == "10€9\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test b[1] == "10€9"

    a = 10.9
    writedlm2("test.csv", a, decimal='.')
    @test read("test.csv", String) == "10.9\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test b[1] == 10.9

    a = [10.9 12.5; Date(2017) DateTime(2017)]
    writedlm2("test.csv", a, ',', decimal='€')
    b = readdlm2("test.csv", ',', decimal='€')
    rm("test.csv")
    @test a == b

    a = [10.9 12.5; Date(2017) DateTime(2017)]
    writedlm2("test.csv", a, ',', decimal='.')
    b = readdlm2("test.csv", ',', decimal='.')
    rm("test.csv")
    @test a == b

    a = [10.9 12.5; Date(2017) DateTime(2017)]
    writedlm2("test.csv", a)
    b = readdlm2("test.csv")
    rm("test.csv")
    @test a == b
end

@testset "2_  date format" begin
    # test date format strings with variable length
    a = DateTime(2017,5,1,5,59,1)
    writedlm2("test.csv", a, dtfs="E, dd.mm.yyyy H:M:S")
    @test read("test.csv", String) == "Monday, 01.05.2017 5:59:1\n"
    b = readdlm2("test.csv", dtfs="E, dd.mm.yyyy H:M:S")
    rm("test.csv")
    @test a == b[1]

    a = DateTime(2017,5,1,5,59,1,898)
    writedlm2("test.csv", a, dtfs="E, d.u yyyy H:M:S,s")
    @test read("test.csv", String) == "Monday, 1.May 2017 5:59:1,898\n"
    b = readdlm2("test.csv", dtfs="E, d.u yyyy H:M:S.s")
    @test a == b[1]
    b = readdlm2("test.csv", dtfs="E, d.u yyyy H:M:S,s")
    @test a == b[1]
    rm("test.csv")

    # test date format strings with fix length
    a = [DateTime(2017,5,1,5,59,1,898) 1.0 1.1 1.222e7 1 true]
    writedlm2("test.csv", a, dtfs="yyyyymmmdddTHHHMMMSSSsss")
    @test read("test.csv", String) == "02017005001T005059001898;1,0;1,1;1,222e7;1;true\n"
    b = readdlm2("test.csv", dtfs="yyyyymmmdddTHHHMMMSSSsss")
    rm("test.csv")
    @test b == a
end

@testset "2_ Time parsing" begin
    a = [Dates.Time(23,55,56,123,456,789) Dates.Time(23,55,56,123,456) Dates.Time(23,55,56,123) Dates.Time(12,45) Dates.Time(11,23,11)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "23:55:56,123456789;23:55:56,123456;23:55:56,123;12:45:00;11:23:11\n"
    b = readdlm2("test.csv")
    @test b == a
    write("test.csv", "23:55:56.123456789;23:55:56,123456;23:55:56,123;12:45;11:23:11\n")
    b = readdlm2("test.csv")
    rm("test.csv")
    @test b == a

    a = [Dates.Time(23,55,56,123,456) Dates.Time(12,45) Dates.Time(11,23,11)]
    writedlm2("test.csv", a, decimal='.')
    @test read("test.csv", String) == "23:55:56.123456;12:45:00;11:23:11\n"
    @test readdlm2("test.csv", dtfs="", dfs="") == ["23:55:56.123456" "12:45:00" "11:23:11"]
    @test readdlm2("test.csv") == a
    rm("test.csv")
end

@testset "2_   DateLocale" begin # Test locale for french and german
    Dates.LOCALES["french"] = Dates.DateLocale(
        ["janvier", "février", "mars", "avril", "mai", "juin",
            "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
        ["janv", "févr", "mars", "avril", "mai", "juin",
            "juil", "août", "sept", "oct", "nov", "déc"],
        ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
        ["lu", "ma", "me", "je", "ve", "sa", "di"],
        )

    Dates.LOCALES["german"] = Dates.DateLocale(
        ["Januar", "Februar", "März", "April", "Mai", "Juni",
            "Juli", "August", "September", "Oktober", "November", "Dezember"],
        ["Jan", "Feb", "Mar", "Apr", "Mai", "Jun",
            "Jul", "Aug", "Sep", "Okt", "Nov", "Dez"],
        ["Montag", "Dienstag", "Mittwoch", "Donnerstag", "Freitag", "Samstag", "Sonntag"],
        ["Mo", "Di", "Mi", "Do", "Fr", "Sa", "So"],
        )

    a = DateTime(2017,5,1,5,59,1)
    writedlm2("test.csv", a, dtfs="E, dd.mm.yyyy H:M:S", locale="french")
    @test read("test.csv", String) == "lundi, 01.05.2017 5:59:1\n"
    b = readdlm2("test.csv", dtfs="E, dd.mm.yyyy H:M:S", locale="french")
    rm("test.csv")
    @test b[1] == a

    a = DateTime(2017,1,1,5,59,1,898)
    writedlm2("test.csv", a, dtfs="E, d.u yyyy H:M:S,s", locale="french")
    @test read("test.csv", String) == "dimanche, 1.janv 2017 5:59:1,898\n"
    b = readdlm2("test.csv", dtfs="E, d.u yyyy H:M:S.s", locale="french")
    @test b[1] == a
    b = readdlm2("test.csv", dtfs="E, d.u yyyy H:M:S,s", locale="french")
    @test b[1] == a
    rm("test.csv")

    a = DateTime(2017,8,1,5,59,1)
    writedlm2("test.csv", a, dtfs="E, dd.mm.yyyy H:M:S", locale="german")
    @test read("test.csv", String) == "Dienstag, 01.08.2017 5:59:1\n"
    b = readdlm2("test.csv", dtfs="E, dd.mm.yyyy H:M:S", locale="german")
    rm("test.csv")
    @test b[1] == a

    a = DateTime(2017,11,1,5,59,1,898)
    writedlm2("test.csv", a, dtfs="E, d. U yyyy H:M:S,s", locale="german")
    @test read("test.csv", String) == "Mittwoch, 1. November 2017 5:59:1,898\n"
    b = readdlm2("test.csv", dtfs="E, d. U yyyy H:M:S.s", locale="german")
    @test b[1] == a
    b = readdlm2("test.csv", dtfs="E, d. U yyyy H:M:S,s", locale="german")
    @test b[1] == a
    rm("test.csv")
end


text:
import Base: sizeof, print, show, convert, read, write, bswap, ==

const Fchar = Int8

"""
    FString{L}

Datatype for reading and writing character strings from `FortranFile`s.
The type parameter `L` signifies the length of the string.
This is the equivalent to the Fortran datatype `CHARACTER(len=L)`.
"""
struct FString{L}
   data :: Array{Fchar,1}
end

sizeof(::Type{FString{L}}) where {L} = L
sizeof(::FString{L}) where {L} = L
sizeof(a::Array{FString{L}}) where {L} = L*length(a)

# print(io::IO, s::FString{L}) where {L} = print(io, trimstring(s))
show(io::IO, s::FString{L}) where {L} = begin print(io, "FString($L,"); show(io, trimstring(s)); print(io, ")") end

==(a::FString, b::FString) = trimstring(a)==trimstring(b)

bswap(s::FString{L}) where {L} = s # no conversion needed for byte-based strings

function convert(::Type{FString{L}}, s::String) where {L}
   spc  = Fchar(' ')
   data = fill(spc, L)
   for (i,c) in enumerate(s)
      if i>L; break; end
      data[i] = Fchar(c)
   end
   FString{L}(data)
end

"""
    FString(L, s::String)

Convert the Julia `String` `s` to an `FString{L}`.
`s` must contain only ASCII characters.
As in Fortran, the string will be padded with spaces or truncated in order to reach the desired length.
"""
FString(L, s::String) = convert( FString{L}, s )

String(s::FString{L}) where {L} = String(map(Char,s.data))
convert(::Type{String}, s::FString{L}) where {L} = String(s)

function read( io::IO, t::Type{FString{L}} ) where {L}
   s = read!(io, Array{Fchar}(undef, L))
   FString{L}(s)
end

function write( io::IO, s::FString{L} ) where {L}
   write(io, s.data)
end

"""
    trimlen(s::FString)

Returns the length of the `FString` `s` with trailing spaces ignored.
"""
function trimlen( s::FString{L} ) where {L}
   l = L
   while l>0
      if s.data[l] != Fchar(' '); break; end
      l -= 1
   end
   return l
end

"""
    trim(s::FString)

Returns a truncated copy of the `FString` `s` where all trailing
spaces are removed.
"""
function trim( s::FString{L} ) where {L}
   l = trimlen(s)
   FString{l}(s.data[1:l])
end

"""
    trimstring(s::FString)

Convert the `FString` `s` into a Julia `String`, where
trailing spaces are removed. Use `String(s)` to keep the spaces.
"""
trimstring( s::FString{L} ) where {L} =
   String( map(Char, s.data[1:trimlen(s)]) )


text:
using Dates
# using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

@testset "RegistryCI.jl" begin
    @testset "RegistryCI.jl unit tests" begin
        @info("Running the RegistryCI.jl unit tests")
        include("registryci_registry_testing.jl")
    end

    @testset "TagBot.jl unit tests" begin
        if v"1.0" <= VERSION < VersionNumber(1, 5, typemax(UInt32))
            @info("Running the TagBot.jl unit tests", VERSION)
            include("tagbot-unit.jl")
        else
            @warn("Skipping the TagBot.jl unit tests", VERSION)
        end
    end

    @testset "AutoMerge.jl unit tests" begin
        @info("Running the AutoMerge.jl unit tests")
        include("automerge-unit.jl")
    end

    AUTOMERGE_RUN_INTEGRATION_TESTS = get(ENV, "AUTOMERGE_RUN_INTEGRATION_TESTS", "")::String
    if AUTOMERGE_RUN_INTEGRATION_TESTS == "true"
        @testset "AutoMerge.jl integration tests" begin
            @info("Running the AutoMerge.jl integration tests")
            include("automerge-integration.jl")
        end
    end
end


text:
# # Cells in Literate Julia
#
# The `{:cell}` syntax for [executable cells](#) is not limited to just
# markdown files. It will work the same across other source types as well.
# This page illustrates this using literate Julia.
#
# To use cells in a Julia file just add `{:cell}` at the end of a comment block
# before the code that you would like to be evaluated. Like so:
#
# {:cell}

f(x) = x + 1
f(1)

# The cells behave the same across all source types and so you will need to
# provide named cells to be able to share values between cells.
#
# {cell=example}

f(x::Integer) = x - 1
length(methods(f))

# As we can see above the function `f` only has 1 method defined since it does
# not know about the previous cell. Below we will define another cell but use
# the same name as before so that we can access it's definitions.
#
# {cell=example}

f(1)

# When the resulting value is all that you care about when writing a cell use
# can just discard the output stream and the cell itself using a little CSS
# by adding a `style="display: none"` to the cell. For example below we use
# the attributes `{cell=example output=false style="display: none"}` to hide
# all but the result. The cell content below is `methods(f)`:
#
# {cell=example output=false style="display: none"}

methods(f)

# Finally, we can also create types that have a `show` method defined for
# `text/markdown`. This allows us to embed arbitrary generated markdown
# within a parent document.
#
# {cell=embedded}

struct Boldify
    x
end

Base.show(io::IO, ::MIME"text/markdown", b::Boldify) = print(io, "**", b.x, "**")

Boldify("This text will appear in bold.")

# {cell=embedded style="display: none"}

using Markdown

# Inline cells are also possible. The source that produces them will not be
# displayed though, only the resulting value. `stdout` and `stderr` will not
# show either. If the resulting value can be displayed as `text/markdown` it's
# AST will be embedded, like above with the block cells, otherwise it will
# appear as a ```md"`text/plain`"```{cell=embedded} wrapped in inline code. The
# previous `text/plain` was embedded in the resulting AST using an inline cell
# that evaluted the expression ```Markdown.md"`text/plain`"```. Values, such as
# `1+2`{:cell} or `pi`{:cell}, will just display as inline code.


text:
using Random
using DABenchmarks
using LaTeXStrings

Random.seed!(0)

using Distributed
machine = [("MACHINENAME", 2)]
addprocs(machine, dir = "/home/USER", exename="/path/to/julia", exeflags="--project")
# addprocs(2, exeflags="--project") # For local workers

@everywhere using StochasticPrograms
@everywhere using ProgressiveHedgingSolvers
@everywhere using HydroModels
@everywhere using Gurobi

function prepare_dbenchmark(nscenarios::Integer, nsamples::Integer; timeout::Integer = 1000)
    ph = ProgressiveHedgingSolver(()->GurobiSolver(OutputFlag=0), execution = :synchronous, penalty = :adaptive, ζ = 0.01, α = 0.75, log = false, τ = 5e-6)
    solvers = [ph]
    solvernames = ["Synchronous progressive-hedging"]
    # Create Day-ahead benchmark
    return Scaling(solvernames, solvers, nscenarios, nsamples; timeout = timeout)
end
@info "Running progressive-hedging benchmarks on 2 workers..."
dab = prepare_dbenchmark(1000, 5); warmup_benchmarks!(dab); run_benchmarks!(dab); save_results!(dab, "ph_2.json")
exit()

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
function chebyshev_polynomial(order::S,x::T) where {T<:Number,S<:Integer}

  polynomial    = Array{T}(undef,1,order+1)
  polynomial[1] = one(T)

  for i = 2:order+1
    if i == 2
      polynomial[i] = x
    else
      polynomial[i] = 2*x*polynomial[i-1]-polynomial[i-2]
    end
  end

  return polynomial

end

function chebyshev_polynomial(order::S,x::AbstractArray{T,1}) where {T<:Number,S<:Integer}

  polynomial      = Array{T}(undef,length(x),order+1)
  polynomial[:,1] = ones(T,length(x))

  for i = 2:order+1
    for j = 1:length(x)
      if i == 2
        polynomial[j,i] = x[j]
      else
        polynomial[j,i] = 2*x[j]*polynomial[j,i-1]-polynomial[j,i-2]
      end
    end
  end

  return polynomial

end

function derivative_of_chebyshev_polynomial(order::S,x::T) where {T<:Number,S<:Integer}

  polynomial    = Array{T}(undef,1,order+1)
  poly_deriv    = Array{T}(undef,1,order+1)
  polynomial[1] = one(T)
  poly_deriv[1] = zero(T)

  for i = 2:order+1
    if i == 2
      polynomial[i] = x
      poly_deriv[i] = one(T)
    else
      polynomial[i] = 2*x*polynomial[i-1]-polynomial[i-2]
      poly_deriv[i] = 2*polynomial[i-1]+2*x*poly_deriv[i-1]-poly_deriv[i-2]
    end
  end

  return poly_deriv

end


text:
# linear regression use case

using LinearAlgebra

include("Observations.jl")

# We have to implement
# (observations: Observations, position: Array{Float64}, term : Int64) -> Float64

"""
# function make\\_linear\\_reg\\_term\\_value(observations:: Observations)

return a closure :

    function linear_reg_term_value(position:: Array{Float64}, term :: Int64)

```math
term value = 1/2 * (y -  observation[term] \\cdot position)^2
```
Our position data is encoded in an array of dimension (1 + length of an observation)

"""
function make_linear_reg_term_value(observations:: Observations)
    function linear_reg_term_value(position:: Array{Float64,1}, term :: Int64)
        x, y = observations.datas[term], observations.value_at_data[term]
        @assert length(position) == length(x) "inequal vector length"
        0.5 * (y - dot(position, x))^2
    end
end

"""
# function make\\_linear\\_reg\\_term\\_gradient(observations::Observations)

returns a closure:

    function linear_reg_term_gradient(position::Array{Float64,1}, term::Int64, gradient::Array{Float64,1})

just take the gradient of the expression for linear\\_reg_term_value
"""
function make_linear_reg_term_gradient(observations::Observations)
    function linear_reg_term_gradient(position::Array{Float64,1}, term::Int64, gradient::Array{Float64,1})
        x, y = observations.datas[term], observations.value_at_data[term]
        e = y - dot(x, position)
    #        @debug "term_gradient e , x term"  e x term
        # it is the .= that forces affectation! Just = does not write into gradient. copy! is clear
        copy!(gradient, -e * x)
    end
end


text:
## This is a working test file for the julia packages installed.

t = Template(;
    user = "keithjlee",
    license = "MIT",
    authors = ["Keith J. Lee"],
    julia_version = v"1.5"
)

generate("TopOptTools", t)


text:
__precompile__()
module TestMetadata
using OpenEphysLoader
using Main.TestUtilities

using Test, Dates, Random

test_dir = joinpath(dirname(@__FILE__), "data")

@testset "Metadata" begin
    info_ans = (
        v"0.3.5",
        v"0",
        DateTime(2015, 7, 21, 14, 50, 12),
        "Windows 7",
        "RIGDESK"
    )
    plugin_info_ans = (
        v"0.4.1",
        v"3",
        DateTime(2016, 10, 21, 14, 19, 52),
        "Windows 8",
        "GLWORKSTATION"
    )

    rhythm_ans = (
        100, # id
        1, # lowcut
        7500, # highcut
        false, # adcs_on
        false, # noiseslicer
        true, # ttl_fastsettle
        false, # dac_ttl
        true, # dac_hpf
        true, # dsp_offset
        0 # dsp_cutoff
    )
    BITVOLTS = 0.19499999284744263
    channel_1_ans = (
        "CH1",
        0,
        BITVOLTS,
        [1024],
        "100_CH1.continuous"
    )
    channel_2_ans = (
        "CH2",
        1,
        BITVOLTS,
        [1024],
        "100_CH2.continuous"
    )

    expermeta_ans = (
        v"0.4.0",
        1,
        false
    )

    recording_ans = (
        0,
        30000
    )

    # pre plugin settings
    expermeta = metadata(test_dir)
    settings = expermeta.settings
    test_fields(settings, (info_ans), recording_chain = false)
    @test ! isempty(settings.recording_chain.nodes)
    rhythmnode = settings.recording_chain.nodes[1].content
    test_fields(rhythmnode, rhythm_ans..., channels = false)
    test_fields(rhythmnode.channels[1], channel_1_ans...)
    test_fields(rhythmnode.channels[2], channel_2_ans...)
    @test expermeta.recordings[1].recording_processors[1] == rhythmnode
    test_fields(expermeta, expermeta_ans..., recordings = false, settings = false)
    @test ! isempty(expermeta.recordings)
    test_fields(expermeta.recordings[1], recording_ans..., recording_processors = false)
    @test (show(devnull, settings.recording_chain); true) # Test that this does not error
    @test (show(devnull, settings); true)
    @test (show(devnull, expermeta.recordings[1]); true)
    @test (show(devnull, expermeta); true)

    #plug in settings
    @test_throws CorruptedException metadata(
        test_dir,
        settingsfile = "plugin_settings.xml",
    )
    @test_throws ErrorException metadata(
        test_dir,
        settingsfile = randstring()
    )
    @test_throws ErrorException metadata(
        test_dir,
        continuousmeta = randstring()
    )

    plugexpermeta = metadata(
        test_dir,
        settingsfile = "plugin_settings.xml",
        continuousmeta = "plugin_Continuous_Data.openephys"
    )
    @test plugexpermeta.settings.info.plugin_api_version == v"3"
end

end # module TestMetaData