text:
"""
    Transducers.ReduceSplitBy(f, rf, init)

Split chunks by elements evaluated to `true` by `f` and reduce each chunk by
reducing function `rf`.

!!! note
    This is an internal implementation detail of [`SplitBy`](@ref) for now.

The reducing function `rf` receives either a `Bulk(x)` if `!f(x)` or a
`End(x)` if `f(x)` returns `true`. Just just after `rf` is called with
`End(x)`, its accumulator is finalized by [`complete`](@ref) and then passed
to the downstream transducer/reducing function.

# Examples
```jldoctest
julia> using Transducers
       using Transducers: ReduceSplitBy, Bulk, End

julia> 1:5 |> ReduceSplitBy(isodd, Map(getindex)'(string), "") |> collect
3-element Vector{String}:
 "1"
 "23"
 "45"

julia> function just_bulk(x)
           if x isa Bulk
               Some(x[])
           else
               nothing
           end
       end;

julia> 1:5 |> ReduceSplitBy(isodd, KeepSomething(just_bulk)'(string), "") |> collect
3-element Vector{String}:
 ""
 "2"
 "4"
```
"""
struct ReduceSplitBy{F,RF,Init} <: Transducer
    f::F
    rf::RF
    init::Init
end

abstract type _Element end
Base.getindex(x::_Element) = x.value

struct Bulk{T} <: _Element
    value::T
end

struct End{T} <: _Element
    value::T
end

struct Chunk{Right}
    right::Right
end
is_prelude(::Chunk) = true

struct Vacant{Left,Right}
    left::Left
    right::Right
end

start(rf::R_{ReduceSplitBy}, acc) = wrap(rf, Chunk(nothing), start(inner(rf), acc))

next(rf::R_{ReduceSplitBy}, acc, input) =
    wrapping(rf, acc) do state, iacc
        if state.right === nothing
            acc0 = start(xform(rf).rf, xform(rf).init)
        else
            acc0 = something(state.right)
        end
        if xform(rf).f(input)
            acc1 = next(xform(rf).rf, acc0, End(input))
            # @show acc1
            if state isa Vacant
                y = complete(xform(rf).rf, acc1)
                iacc′ = next(inner(rf), iacc, y)
                left = state.left
            else
                iacc′ = iacc
                left = acc1
            end
            (Vacant(left, nothing), iacc′)
        else
            right = Some(next(xform(rf).rf, acc0, Bulk(input)))
            ((@set state.right = right), iacc)
        end
    end

function complete(rf::R_{ReduceSplitBy}, acc)
    state, iacc1 = unwrap(rf, acc)
    if state.right === nothing
        iacc2 = iacc1
    else
        yr = complete(xform(rf).rf, something(state.right))
        iacc2 = next(inner(rf), iacc1, yr)
    end
    if state isa Vacant
        # @show state.left
        yl = complete(inner(rf), state.left)
        # TODO: Is it OK to always use DefaultInit? Is there a better way?
        iacc0 = next(inner(rf), start(inner(rf), DefaultInit), yl)
        iacc3 = combine(inner(rf), iacc0, iacc2)
    else
        iacc3 = iacc2
    end
    return complete(inner(rf), iacc3)
end

@inline maybe_combine(_, ::Nothing, ::Nothing) = nothing
@inline maybe_combine(_, x::Some, ::Nothing) = x
@inline maybe_combine(_, ::Nothing, x::Some) = x
@inline maybe_combine(rf::RF, a::Some, b::Some) where {RF} =
    Some(combine(rf, something(a), something(b)))

function combine(rf::R_{ReduceSplitBy}, a, b)
    a1, a2 = unwrap(rf, a)
    b1, b2 = unwrap(rf, b)
    if a1 isa Chunk
        if b1 isa Chunk
            c1 = Chunk(maybe_combine(xform(rf).rf, a1.right, b1.right))
        else
            left = something(maybe_combine(xform(rf).rf, a1.right, Some(b1.left)))
            c1 = Vacant(left, b1.right)
        end
    else
        if b1 isa Chunk
            c1 = Vacant(a1.left, maybe_combine(xform(rf).rf, a1.right, b1.right))
        else
            acc = something(maybe_combine(xform(rf).rf, a1.right, Some(b1.left)))
            y = complete(xform(rf).rf, acc)
            a2 = next(inner(rf), a2, y)
            c1 = Vacant(a1.left, b1.right)
        end
    end
    c2 = combine(inner(rf), a2, b2)
    return wrap(rf, c1, c2)
end

"""
    SplitBy(f; [keepend = false,] [keepempty = false,])

Split input collection into chunks delimited by the elements on which `f`
returns `true`. This can be used to implement parallel and lazy versions of
functions like `eachline` and `split`.

If `keepend` is `true` (or `Val(true)`), include the "delimiter"/end element
at the end of each chunk. If `keepempty` is `true` (or `Val(true)`), include
empty chunks. When `keepend` is `true`, the value of `keepempty` is
irrelevant since the chunks cannot be empty (i.e., it at least contains the
end).

The input collection (`xs` in `SplitBy(...)(xs)`) has to support `eachindex`
and `view` or `SubString`.

# Extended Help

## Examples

For demonstration, consider the following input stream and `SplitBy(iszero;
...)` used with the following options:

```
input     keepend=false       keepend=false        keepend=true
          keepempty=false     keepempty=true

1           `.                  `.                  `.
2            | y1                | y1                | y1
3           ,'                  ,'                   |
0                                                   ,'
1           `.                  `.                  `.
2            | y2                | y2                | y2
3            |                   |                   |
4           ,'                  ,'                   |
0                               __ y3               ,'
0                                                    ] y3
1           `.                  `.                  `.
2            | y3                | y4                | y4
```

In the above diagram, `yi` (`i = 1, 2, 3, 4`) are the output of `SplitBy`. This
can be checked in the REPL easily as follows. (Note: we are using
`Map(collect)` for cleaner printing; it's not required unless the arrays is
mutated in the downstream.)

```jldoctest
julia> using Transducers

julia> xs = [1, 2, 3, 0, 1, 2, 3, 4, 0, 0, 1, 2];  # input

julia> xs |> SplitBy(iszero) |> Map(collect) |> collect
3-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [1, 2, 3, 4]
 [1, 2]

julia> xs |> SplitBy(iszero; keepempty = true) |> Map(collect) |> collect
4-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [1, 2, 3, 4]
 []
 [1, 2]

julia> xs |> SplitBy(iszero; keepend = true) |> Map(collect) |> collect
4-element Vector{Vector{Int64}}:
 [1, 2, 3, 0]
 [1, 2, 3, 4, 0]
 [0]
 [1, 2]
```
"""
struct SplitBy{F,KeepEnd<:ValBool,KeepEmpty<:ValBool} # <: Transducer  # TODO: should I?
    f::F
    keepend::KeepEnd
    keepempty::KeepEmpty
end

SplitBy(
    f;
    keepend::Union{Bool,ValBool} = Val(false),
    keepempty::Union{Bool,ValBool} = Val(false),
) = SplitBy(f, asval(keepend), asval(keepempty))

struct SplitByReducingFunction{KeepEnd<:ValBool,KeepEmpty<:ValBool}
    keepend::KeepEnd
    keepempty::KeepEmpty
end

function next(rf::SplitByReducingFunction, a, b)
    # @show a b
    j = b[]
    if a === nothing
        i = j
        k = j - 1
    else
        i, k, = a::Tuple{Any,Any}
    end
    if b isa Bulk
        return (i, j)
    else
        b::End
        if rf.keepend === Val(true)
            return (i, j)
        elseif rf.keepempty === Val(true)
            return (i, k)
        elseif a === nothing
            return nothing
        else
            return (i, k)
        end
    end
end

function combine(::SplitByReducingFunction, a, b)
    a === nothing && return b
    b === nothing && return a
    return (a[1], b[2])
end

complete(::SplitByReducingFunction, a) = a

(f::SplitBy)(xs) =
    xs |>
    eachindex |>
    ReduceSplitBy(
        i -> f.f(@inbounds xs[i]),
        SplitByReducingFunction(f.keepend, f.keepempty),
        nothing,
    ) |>
    NotA(Nothing) |>
    Map(sb_viewer(xs))
# TODO: Check if it is better to use `pairs`; note that it'd require
# `pairs(::AbstractString)` to return some structured iterator.

sb_viewer(xs) = view_any((i, j),) = view(xs, i:j)
sb_viewer(xs::AbstractString) = view_str((i, j),) = SubString(xs, i, j)


text:
# Knitro model attributes

##################################################
# objective
##################################################
# set objective sense
"Set objective goal."
function KN_set_obj_goal(m::Model, objgoal::Cint)
    ret = @kn_ccall(set_obj_goal, Cint, (Ptr{Cvoid}, Cint),
                    m.env, objgoal)
    _checkraise(ret)
end

"""
Add linear structure to the objective function.
Each component i of arrays indexVars and coefs adds a linear term
   coefs[i]*x[indexVars[i]]
to the objective.

"""
function KN_add_obj_linear_struct(m::Model,
                                  objIndices::Vector{Cint},
                                  objCoefs::Vector{Cdouble})
    nnz = length(objIndices)
    @assert nnz == length(objCoefs)
    ret = @kn_ccall(add_obj_linear_struct, Cint,
                    (Ptr{Cvoid}, KNLONG, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    objIndices,
                    objCoefs)
    _checkraise(ret)
end
KN_add_obj_linear_struct(m::Model, objindex::Int, objCoefs::Cdouble) =
    KN_add_obj_linear_struct(m, Cint[objindex], [objCoefs])

# quadratic part of objective
"""
Add quadratic structure to the objective function.
Each component i of arrays indexVars1, indexVars2 and coefs adds a quadratic
term
   coefs[i]*x[indexVars1[i]]*x[indexVars2[i]]
to the objective.

Note: if indexVars2[i] is < 0 then it adds a linear term
      coefs[i]*x[indexVars1[i]] instead.

"""
function KN_add_obj_quadratic_struct(m::Model,
                                     indexVars1::Vector{Cint},
                                     indexVars2::Vector{Cint},
                                     coefs::Vector{Cdouble})
    nnz = length(indexVars1)
    @assert nnz == length(indexVars2) == length(coefs)
    ret = @kn_ccall(add_obj_quadratic_struct, Cint,
                    (Ptr{Cvoid}, KNLONG, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    indexVars1,
                    indexVars2,
                    coefs)
    _checkraise(ret)
end

function KN_add_obj_constant(m::Model, constant::Cdouble)
    ret = @kn_ccall(add_obj_constant, Cint, (Ptr{Cvoid}, Cdouble),
                    m.env, constant)
    _checkraise(ret)
end

function KN_set_obj_scaling(m::Model, objScaleFactor::Cdouble)
    ret = @kn_ccall(set_obj_scaling, Cint, (Ptr{Cvoid}, Cdouble),
                    m.env, objScaleFactor)
    _checkraise(ret)
end

"""
Specify some properties of the objective and constraint functions.
Note: use bit-wise specification of the features:
bit value   meaning
  0     1   KN_OBJ_CONVEX/KN_CON_CONVEX
  1     2   KN_OBJ_CONCAVE/KN_CON_CONCAVE
  2     4   KN_OBJ_CONTINUOUS/KN_CON_CONTINUOUS
  3     8   KN_OBJ_DIFFERENTIABLE/KN_CON_DIFFERENTIABLE
  4    16   KN_OBJ_TWICE_DIFFERENTIABLE/KN_CON_TWICE_DIFFERENTIABLE
  5    32   KN_OBJ_NOISY/KN_CON_NOISY
  6    64   KN_OBJ_NONDETERMINISTIC/KN_CON_NONDETERMINISTIC

"""
function KN_set_obj_property(m::Model, objProperty::Cint)
    ret = @kn_ccall(set_obj_property, Cint, (Ptr{Cvoid}, Cint),
                    m.env, objProperty)
    _checkraise(ret)
end

function KN_set_obj_name(m::Model, name::AbstractString)
    ret = @kn_ccall(set_obj_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}),
                    m.env, name)
    _checkraise(ret)
end

##################################################
# Modifiers
##################################################
if KNITRO_VERSION >= v"12.4"
    function KN_del_obj_constant(kc::Model)
        ret = @kn_ccall(del_obj_constant, Cint, (Ptr{Cvoid},), kc.env)
        _checkraise(ret)
    end

    function KN_chg_obj_constant(kc::Model, constant::Cdouble)
        ret = @kn_ccall(chg_obj_constant, Cint, (Ptr{Cvoid}, Cdouble), kc.env, constant)
        _checkraise(ret)
    end

    function KN_add_obj_linear_term(kc::Model, indexVar::Cint, coef::Cdouble)
        ret = @kn_ccall(add_obj_linear_term, Cint, (Ptr{Cvoid}, Cint, Cdouble), kc.env, indexVar, coef)
        _checkraise(ret)
    end

    function KN_chg_obj_linear_term(kc::Model, indexVar::Cint, coef::Cdouble)
        ret = @kn_ccall(chg_obj_linear_term, Cint, (Ptr{Cvoid}, Cint, Cdouble), kc.env, indexVar, coef)
        _checkraise(ret)
    end

    function KN_del_obj_linear_term(kc::Model, indexVar::Cint)
        ret = @kn_ccall(del_obj_linear_term, Cint, (Ptr{Cvoid}, Cint), kc.env, indexVar)
        _checkraise(ret)
    end

    function KN_del_obj_linear_struct(kc::Model, indexVars::Vector{Cint})
        ret = @kn_ccall(del_obj_linear_struct, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}),
                        kc.env, length(indexVars), indexVars)
        _checkraise(ret)
    end

    function KN_chg_obj_linear_struct(kc::Model, indexVars::Vector{Cint}, coefs::Vector{Cdouble})
        nnz = length(indexVars)
        @assert nnz == length(coefs)
        ret = @kn_ccall(chg_obj_linear_struct, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                        kc.env, nnz, indexVars, coefs)
        _checkraise(ret)
    end

    function KN_add_obj_quadratic_term(kc::Model, indexVar1::Cint, indexVar2::Cint, coef::Cdouble)
        ret = @kn_ccall(add_obj_quadratic_term, Cint, (Ptr{Cvoid}, Cint, Cint, Cdouble), kc.env, indexVar1, indexVar2, coef)
        _checkraise(ret)
    end
end

##################################################
# Generic getters
##################################################
function KN_get_number_vars(kc::Ptr{Cvoid})
    num_vars = Cint[0]
    ret = @kn_ccall(get_number_vars, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    kc, num_vars)
    _checkraise(ret)
    return num_vars[1]
end
KN_get_number_vars(model::Model) = KN_get_number_vars(model.env.ptr_env)

function KN_get_number_cons(kc::Ptr{Cvoid})
    num_cons = Cint[0]
    ret = @kn_ccall(get_number_cons, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    kc, num_cons)
    _checkraise(ret)
    return num_cons[1]
end
KN_get_number_cons(model::Model) = KN_get_number_cons(model.env.ptr_env)

function KN_get_obj_value(m::Model)
    obj = Cdouble[0]
    ret = @kn_ccall(get_obj_value, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, obj)
    _checkraise(ret)
    return obj[1]
end

function KN_get_obj_type(m::Model)
    obj_type = Cint[0]
    ret = @kn_ccall(get_obj_type, Cint, (Ptr{Cvoid}, Ptr{Cint}), m.env, obj_type)
    _checkraise(ret)
    return obj_type[1]
end

##################################################
# Constraints getters
##################################################
function KN_get_con_values(m::Model)
    nc = KN_get_number_cons(m)
    consvals = zeros(Cdouble, nc)
    ret = @kn_ccall(get_con_values_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, consvals)
    _checkraise(ret)
    return consvals
end

function KN_get_con_values(m::Model, cIndex::Integer)
    nc = 1
    consvals = zeros(Cdouble, nc)
    ret = @kn_ccall(get_con_value, Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}),
                    m.env, cIndex, consvals)
    _checkraise(ret)
    return consvals[1]
end

function KN_get_con_types(m::Model)
    nc = KN_get_number_cons(m)
    constypes = zeros(Cint, nc)
    ret = @kn_ccall(get_con_types_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, constypes)
    _checkraise(ret)
    return constypes
end

##################################################
# Continuous optimization results
##################################################
"""
Return the number of iterations made by KN_solve in "numIters".
"""
function KN_get_number_iters(m::Model)
    num_iters = Cint[0]
    ret = @kn_ccall(get_number_iters, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, num_iters)
    _checkraise(ret)
    return num_iters[1]
end

"""
Return the number of conjugate gradient (CG) iterations made by
KN_solve in "numCGiters".

"""
function KN_get_number_cg_iters(m::Model)
    num_iters = Cint[0]
    ret = @kn_ccall(get_number_cg_iters, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, num_iters)
    _checkraise(ret)
    return num_iters[1]
end

"""
Return the absolute feasibility error at the solution in "absFeasError".
Refer to the Knitro manual section on Termination Tests for a
detailed definition of this quantity.

"""
function KN_get_abs_feas_error(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_abs_feas_error, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the relative feasibility error at the solution in "relFeasError".
Refer to the Knitro manual section on Termination Tests for a
detailed definition of this quantity.

"""
function KN_get_rel_feas_error(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_rel_feas_error, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the absolute optimality error at the solution in "absOptError".
Refer to the Knitro manual section on Termination Tests for a
detailed definition of this quantity.

"""
function KN_get_abs_opt_error(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_abs_opt_error, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the relative optimality error at the solution in "relOptError".
Refer to the Knitro manual section on Termination Tests for a
detailed definition of this quantity.

"""
function KN_get_rel_opt_error(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_rel_opt_error, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

##################################################
# Fetch solution utils
##################################################
#--------------------
# Objective gradient
#--------------------
function KN_get_objgrad_nnz(m::Model)
    res = Cint[0]
    ret = @kn_ccall(get_objgrad_nnz, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_objgrad_values(m::Model)
    nnz = KN_get_objgrad_nnz(m)
    indexVars = zeros(Cint, nnz)
    objGrad = zeros(Cdouble, nnz)
    ret = @kn_ccall(get_objgrad_values, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, indexVars, objGrad)
    _checkraise(ret)
    return indexVars, objGrad
end

#--------------------
# Jacobian
#--------------------
function KN_get_jacobian_nnz(m::Model)
    res = KNLONG[0]
    ret = @kn_ccall(get_jacobian_nnz, Cint, (Ptr{Cvoid}, Ptr{KNLONG}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_jacobian_values(m::Model)
    nnz = KN_get_jacobian_nnz(m)
    jacvars = zeros(Cint, nnz)
    jaccons = zeros(Cint, nnz)
    jaccoef = zeros(Cdouble, nnz)
    ret = @kn_ccall(get_jacobian_values, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, jacvars, jaccons, jaccoef)
    _checkraise(ret)
    return jacvars, jaccons, jaccoef
end

#--------------------
# Rsd Jacobian
#--------------------
function KN_get_rsd_jacobian_nnz(m::Model)
    res = KNLONG[0]
    ret = @kn_ccall(get_rsd_jacobian_nnz, Cint, (Ptr{Cvoid}, Ptr{KNLONG}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the values of the residual Jacobian in "indexRsds", "indexVars",
and "rsdJac".  The Jacobian values returned correspond to the non-zero
sparse Jacobian indices provided by the user.

"""
function KN_get_rsd_jacobian_values(m::Model)
    nnz = KN_get_rsd_jacobian_nnz(m)
    jacvars = zeros(Cint, nnz)
    jaccons = zeros(Cint, nnz)
    jaccoef = zeros(Cdouble, nnz)
    ret = @kn_ccall(get_rsd_jacobian_values, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, jacvars, jaccons, jaccoef)
    _checkraise(ret)
    return jacvars, jaccons, jaccoef
end

#--------------------
# Hessian
#--------------------
function KN_get_hessian_nnz(m::Model)
    res = KNLONG[0]
    ret = @kn_ccall(get_hessian_nnz, Cint, (Ptr{Cvoid}, Ptr{KNLONG}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_hessian_values(m::Model)
    nnz = KN_get_hessian_nnz(m)
    indexVars1 = zeros(Cint, nnz)
    indexVars2 = zeros(Cint, nnz)
    hess = zeros(Cdouble, nnz)
    ret = @kn_ccall(get_hessian_values, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, indexVars1, indexVars2, hess)
    _checkraise(ret)
    return indexVars1, indexVars2, hess
end

# Getters for CPU time are implemented only for Knitro version >= 12.0.
if KNITRO_VERSION >= v"12.0"
    function KN_get_solve_time_cpu(m::Model)
        tcpu = zeros(Cdouble, 1)
        ret = @kn_ccall(get_solve_time_cpu, Cint,
                        (Ptr{Cvoid}, Ptr{Cdouble}), m.env, tcpu)
        _checkraise(ret)
        return tcpu[1]
    end

    function KN_get_solve_time_real(m::Model)
        treal = zeros(Cdouble, 1)
        ret = @kn_ccall(get_solve_time_real, Cint,
                        (Ptr{Cvoid}, Ptr{Cdouble}), m.env, treal)
        _checkraise(ret)
        return treal[1]
    end
end
##################################################
# MIP utils
##################################################
#--------------------
# Getters
#--------------------
"""
Return the number of nodes processed in the MIP solve
in "numNodes".

"""
function KN_get_mip_number_nodes(m::Model)
    res = Cint[0]
    ret = @kn_ccall(get_mip_number_nodes, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the number of continuous subproblems processed in the
MIP solve in "numSolves".

"""
function KN_get_mip_number_solves(m::Model)
    res = Cint[0]
    ret = @kn_ccall(get_mip_number_solves, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the final absolute integrality gap in the MIP solve
in "absGap". Refer to the Knitro manual section on Termination
Tests for a detailed definition of this quantity. Set to
KN_INFINITY if no incumbent (i.e., integer feasible) point found.

"""
function KN_get_mip_abs_gap(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_abs_gap, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the final absolute integrality gap in the MIP solve
int "relGap". Refer to the Knitro manual section on Termination
Tests for a detailed definition of this quantity.  Set to
KN_INFINITY if no incumbent (i.e., integer feasible) point found.

"""
function KN_get_mip_rel_gap(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_rel_gap, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the objective value of the MIP incumbent solution in
"incumbentObj". Set to KN_INFINITY if no incumbent (i.e., integer
feasible) point found.

"""
function KN_get_mip_incumbent_obj(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_incumbent_obj, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the value of the current MIP relaxation bound in "relaxBound".

"""
function KN_get_mip_relaxation_bnd(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_relaxation_bnd, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"""
Return the objective value of the most recently solved MIP
node subproblem in "lastNodeObj".

"""
function KN_get_mip_lastnode_obj(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_lastnode_obj, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

"Return the MIP incumbent solution in 'x' if one exists."
function KN_get_mip_incumbent_x(m::Model)
    res = Cdouble[0]
    ret = @kn_ccall(get_mip_incumbent_x, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, res)
    _checkraise(ret)
    return res[1]
end

#--------------------
# Branching priorities
#--------------------
"""
Set the branching priorities for integer variables. Must first
set the types of variables (e.g. by calling KN_set_var_types) before
calling this function. Priorities must be positive numbers
(variables with non-positive values are ignored).  Variables with
higher priority values will be considered for branching before
variables with lower priority values.  When priorities for a subset
of variables are equal, the branching rule is applied as a tiebreaker.
Values for continuous variables are ignored.  Knitro makes a local
copy of all inputs, so the application may free memory after the call.

"""
function KN_set_mip_branching_priorities(m::Model, nindex::Integer, xPriorities::Cint)
    ret = @kn_ccall(set_mip_branching_priority, Cint,
                    (Ptr{Cvoid}, Cint, Cint),
                    m.env, nindex, xPriorities)
    _checkraise(ret)
end

function KN_set_mip_branching_priorities(m::Model, xIndex::Vector{Cint}, xPriorities::Vector{Cint})
    nvar = length(xIndex)
    @assert nvar == length(xPriorities)
    ret = @kn_ccall(set_mip_branching_priorities, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}),
                    m.env, nvar, xIndex, xPriorities)
    _checkraise(ret)
end

function KN_set_mip_branching_priorities(m::Model, xPriorities::Vector{Cint})
    ret = @kn_ccall(set_mip_branching_priorities_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, xPriorities)
    _checkraise(ret)
end

#--------------------
# Intvar strategies
#--------------------
"""
Set strategies for dealing with individual integer variables. Possible
strategy values include:
  KN_MIP_INTVAR_STRATEGY_NONE    0 (default)
  KN_MIP_INTVAR_STRATEGY_RELAX   1
  KN_MIP_INTVAR_STRATEGY_MPEC    2 (binary variables only)
indexVars should be an index value corresponding to an integer variable
(nothing is done if the index value corresponds to a continuous variable),
and xStrategies should correspond to one of the strategy values listed above.

"""
function KN_set_mip_intvar_strategies(m::Model, nindex::Integer, xStrategies::Cint)
    ret = @kn_ccall(set_mip_intvar_strategy, Cint,
                    (Ptr{Cvoid}, Cint, Cint),
                    m.env, nindex, xStrategies)
    _checkraise(ret)
end

function KN_set_mip_intvar_strategies(m::Model, xIndex::Vector{Cint}, xStrategies::Vector{Cint})
    nvar = length(xIndex)
    @assert nvar == length(xStrategies)
    ret = @kn_ccall(set_mip_intvar_strategies, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}),
                    m.env, nvar, xIndex, xStrategies)
    _checkraise(ret)
end

function KN_set_mip_intvar_strategies(m::Model, xStrategies::Vector{Cint})
    ret = @kn_ccall(set_mip_intvar_strategies_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, xStrategies)
    _checkraise(ret)
end

##################################################
# Parameters
##################################################
#------------------------------
# Setters
#------------------------------
# Int params
function KN_set_param(m::Model, id::Integer, value::Integer)
    ret = @kn_ccall(set_int_param, Cint, (Ptr{Cvoid}, Cint, Cint),
                    m.env, id, value)
    _checkraise(ret)
end

function KN_set_param(m::Model, param::AbstractString, value::Integer)
    ret = @kn_ccall(set_int_param_by_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Cint),
                    m.env, param, value)
    _checkraise(ret)
end

# Double params
function KN_set_param(m::Model, id::Integer, value::Cdouble)
    ret = @kn_ccall(set_double_param, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, id, value)
    _checkraise(ret)
end

function KN_set_param(m::Model, param::AbstractString, value::Cdouble)
    ret = @kn_ccall(set_double_param_by_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Cdouble),
                    m.env, param, value)
    _checkraise(ret)
end

# Char params
function KN_set_param(m::Model, id::Integer, value::AbstractString)
    ret = @kn_ccall(set_char_param, Cint, (Ptr{Cvoid}, Cint, Ptr{Cchar}),
                    m.env, id, value)
    _checkraise(ret)
end

function KN_set_param(m::Model, param::AbstractString, value::AbstractString)
    ret = @kn_ccall(set_char_param_by_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cchar}),
                    m.env, param, value)
    _checkraise(ret)
end

#------------------------------
# Getters
#------------------------------

# Int params
function KN_get_int_param(m::Model, id::Integer)
    res = Cint[0]
    ret = @kn_ccall(get_int_param, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}),
                    m.env, id, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_int_param(m::Model, param::AbstractString)
    res = Cint[0]
    ret = @kn_ccall(get_int_param_by_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cint}),
                    m.env, param, res)
    _checkraise(ret)
    return res[1]
end

# Double params
function KN_get_double_param(m::Model, id::Integer)
    res = Cdouble[0.]
    ret = @kn_ccall(get_double_param, Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}),
                    m.env, id, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_double_param(m::Model, param::AbstractString)
    res = Cdouble[0.]
    ret = @kn_ccall(get_double_param_by_name, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cdouble}),
                    m.env, param, res)
    _checkraise(ret)
    return res[1]
end

#------------------------------
# Params information
#------------------------------
function KN_get_param_name(m::Model, id::Integer)
    output_size = 128
    res = " "^output_size
    ret = @kn_ccall(get_param_name, Cint, (Ptr{Cvoid}, Cint, Ptr{Cchar}, Csize_t),
                    m.env, id, res, output_size)
    _checkraise(ret)
    return format_output(res)
end

function KN_get_param_doc(m::Model, id::Integer)
    output_size = 128
    res = " "^output_size
    ret = @kn_ccall(get_param_doc, Cint, (Ptr{Cvoid}, Cint, Ptr{Cchar}, Csize_t),
                    m.env, id, res, output_size)
    _checkraise(ret)
    return format_output(res)
end

function KN_get_param_type(m::Model, id::Integer)
    res = Cint[0]
    ret = @kn_ccall(get_param_type, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}),
                    m.env, id, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_num_param_values(m::Model, id::Integer)
    res = Cint[0]
    ret = @kn_ccall(get_num_param_values, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}),
                    m.env, id, res)
    _checkraise(ret)
    return res[1]
end

function KN_get_param_value_doc(m::Model, id::Integer, value_id::Integer)
    output_size = 128
    res = " "^output_size
    ret = @kn_ccall(get_param_value_doc, Cint,
                    (Ptr{Cvoid}, Cint, Cint, Ptr{Cchar}, Csize_t),
                    m.env, id, value_id, res, output_size)
    _checkraise(ret)
    return format_output(res)
end

function KN_get_param_id(m::Model, name::AbstractString)
    res = Cint[0]
    ret = @kn_ccall(get_param_id, Cint, (Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cint}),
                    m.env, name, res)
    _checkraise(ret)
    return res[1]
end

function Base.show(io::IO, m::Model)
    if is_valid(m)
        println(io, "$(get_release())")
        println(io, "-----------------------")
        println(io, "Problem Characteristics")
        println(io, "-----------------------")
        println(io, "Objective goal:  Minimize")
        println(io, "Objective type:  $(KN_get_obj_type(m))")
        println(io, "Number of variables:                             $(KN_get_number_vars(m))")
        println(io, "Number of constraints:                           $(KN_get_number_cons(m))")
        println(io, "Number of nonzeros in Jacobian:                  $(KN_get_jacobian_nnz(m))")
        println(io, "Number of nonzeros in Hessian:                   $(KN_get_hessian_nnz(m))")

    else
        println(io, "KNITRO Problem: NULL")
    end
end

ext: jl


text:
@testset "592.fraction-addition-and-subtraction.jl" begin
    @test fraction_addition("-1/2+1/2") == "0/1"
    @test fraction_addition("-1/2+1/2+1/3") == "1/3"
    @test fraction_addition("1/3-1/2") == "-1/6"
    @test fraction_addition("5/3+1/3") == "2/1"
end


text:
export TrivialGate, IdentityGate, igate

abstract type TrivialGate{D} <: PrimitiveBlock{D} end

mat(::Type{T}, d::TrivialGate{D}) where {T,D} = IMatrix{D^nqudits(d),T}()
Base.adjoint(g::TrivialGate) = g
occupied_locs(g::TrivialGate) = ()

"""
    IdentityGate{D} <: TrivialGate{D}

The identity gate.
"""
struct IdentityGate{D} <: TrivialGate{D}
    n::Int
end
nqudits(ig::IdentityGate) = ig.n

"""
    igate(n::Int; nlevel=2)

The constructor for identity gate.
"""
igate(n::Int; nlevel=2) = IdentityGate{nlevel}(n)


text:
### A Nbdev script file ###
### Autogenerated file. Don't modify. ###

module Make
#export
using OrderedCollections

#export
include("../src/Documenter.jl")

#export
import YAML

#export
"""
> build(nbsdir)--> extracts markdown items and builds documentation from the notebooks in the given directory. Also, creates creates a readme file with content in the project root.
"""
function build(nbsdir)
    Documenter.export2md(nbsdir)
	Documenter.export2readme()
end

#export
mkdocs = OrderedDict("site-name" => "your site name", "nav" => [OrderedDict("Home" => "index.md")],
"theme" => OrderedDict("name" => "material", "palette" => OrderedDict("primary" => "white"), "font" => OrderedDict("text" => "Roboto"), "logo" => "assets/asset image filename", "favicon" => "asset/ asset name", "features" => "navigation.tabs", "markdown_extensions" => ["meta", "toc"], "plugins" => ["search", "tags"]), "extra" => OrderedDict("social" => [OrderedDict("icon" => "fontawsome/brans/twitter", "link" => "twitter link"), OrderedDict("icon" => "fotawsome/brands/linkedin", "link" => "linkedin link")]), "repo_url" => "url to your git repo", "repo_name" => "name of your git repo")

#export
begin
"""
> new(projname)-> Scaffolds a new project directory with the necessary sub-directories and configuration files.
"""
#TODO: creation of assets, iamges, index.jl is not getting
#created when this function is called from terminal
function new(projname = "newproj")
	subdirs = ["$projname/docs", "nbs", "$projname/src", "$projname/docs/images", "$projname/docs/assets"]
	files = ["$projname/Manifest.toml", "$projname/Project.toml", "$projname/mkdocs.yml", "$projname/nbs/index.jl"]
	map(mkpath, subdirs)
	map(touch, files)
	YAML.write_file(files[3], mkdocs)
end
end

end


text:
using Test
using NeXLMatrixCorrection

include("xpp.jl")
include("citzaf.jl")
include("riveros.jl")
include("iterate.jl")
include("xppu.jl")
include("aspure.jl")


text:
__precompile__(true)
module Alazar

# Type aliases go here
export U32, U16, U8, S32, S16
export dsp_module_handle
export AlazarBits, Alazar8Bit, Alazar12Bit, Alazar16Bit
export AlazarFFTBits, U16Log, U16Amp2, U8Log, U8Amp2
export S32Real, S32Imag, FloatLog, FloatAmp2

export alazaropen

const U32 = Culong
const U16 = Cushort
const U8  = Cuchar
const S32 = Clong
const S16 = Cshort
const dsp_module_handle = Ptr{Void}

abstract type AlazarBits end
abstract type AlazarFFTBits <: AlazarBits end

struct Alazar8Bit <: AlazarBits
    b::U8
end
struct Alazar12Bit <: AlazarBits
    b::U16
end
struct Alazar16Bit <: AlazarBits
    b::U16
end
struct U16Log <: AlazarFFTBits
    b::U16
end
struct U16Amp2 <: AlazarFFTBits
    b::U16
end
struct U8Log <: AlazarFFTBits
    b::U8
end
struct U8Amp2 <: AlazarFFTBits
    b::U8
end
struct S32Real <: AlazarFFTBits
    b::S32
end
struct S32Imag <: AlazarFFTBits
    b::S32
end
struct FloatLog <: AlazarFFTBits
    b::Cfloat
end
struct FloatAmp2 <: AlazarFFTBits
    b::Cfloat
end

Base.convert(::Type{UInt8}, x::Alazar8Bit) = ltoh(convert(UInt8, x.b))
Base.convert(::Type{T}, x::Alazar8Bit) where {T <: Integer} = convert(T, convert(UInt8, x))
Base.convert(::Type{T}, x::Alazar8Bit) where {T <: AbstractFloat} = T(convert(UInt8,x)/0xFF*2-1)

Base.convert(::Type{UInt16}, x::Alazar12Bit) = ltoh(convert(UInt16, x.b)) >> 4
Base.convert(::Type{T}, x::Alazar12Bit) where {T <: Integer} = convert(T, convert(UInt16, x))
Base.convert(::Type{T}, x::Alazar12Bit) where {T <: AbstractFloat} = T(convert(UInt16,x)/0xFFF*2-1)

Base.convert(::Type{UInt16}, x::Alazar16Bit) = ltoh(convert(UInt16, x.b))
Base.convert(::Type{T}, x::Alazar16Bit) where {T <: Integer} = convert(T, convert(UInt16, x))
Base.convert(::Type{T}, x::Alazar16Bit) where {T <: AbstractFloat} = T(convert(UInt16,x)/0xFFFF*2-1)

Base.convert(::Type{S}, x::T) where {S <: Real,T <: AlazarBits} = convert(S, x.b)
Base.convert(::Type{S}, x::T) where {S <: Complex,T <: S32Real} = S(x.b,0)
Base.convert(::Type{S}, x::T) where {S <: Complex,T <: S32Imag} = S(0,x.b)
Base.convert(::Type{T}, x::T) where {T <: AlazarBits} = x
Base.convert(::Type{S}, x::T) where {S <: AlazarBits,T <: AlazarBits} = S(x.b)

# Constants and exceptions go here
include("constants.jl")

# Load libraries
const ats = is_windows() ? "ATSApi.dll" : "libATSApi.so"
const linux_libc = "libc.so.6"
# The library should only be loaded once. We don't load it automatically because
# if another worker process loads this module, there could be problems.
#
# Instead the library should be given the chance to load (if it hasn't already)
# whenever Julia objects representing instruments are created.
alazaropen() = begin
    @static is_windows() ? begin
        atsHandle = Libdl.dlopen(ats)
        atexit(()->Libdl.dlclose(atsHandle))
    end : (@static is_linux() ? begin
        atsHandle = Libdl.dlopen(ats)
        libcHandle = Libdl.dlopen(linux_libc)
        atexit(()->begin
            Libdl.dlclose(atsHandle)
            Libdl.dlclose(libcHandle)
        end)
    end : error("OS not supported"))
end

include("buffer.jl")
include("errors.jl")
include("api.jl")
include("dsp.jl")

end


text:
function areequal(theory::Vector{<:Rule}, exprs...;
    mod=@__MODULE__, params=SaturationParams())
    g = EGraph(exprs[1])
    areequal(g, theory, exprs...; params=params)
end

function areequal(g::EGraph, t::Vector{<:Rule}, exprs...;
    mod=@__MODULE__, params=SaturationParams())
    @log "Checking equality for " exprs
    if length(exprs) == 1; return true end
    # rebuild!(G)

    @log "starting saturation"

    goal = EqualityGoal(g, collect(exprs))

    # alleq = () -> (all(x -> in_same_set(G.uf, ids[1], x), ids[2:end]))

    params.goal = goal
    # params.stopwhen = alleq

    report = saturate!(g, t, params; mod=mod)

    # display(g.classes); println()
    if !(report.reason isa Saturated) && !reached(g, goal)
        return missing # failed to prove
    end
    return reached(g, goal)
end

import ..gettheory

macro areequal(theory, exprs...)
    t = gettheory(theory, __module__; compile=false)
    areequal(t, exprs...; mod=__module__)
end

macro areequalg(G, theory, exprs...)
    t = gettheory(theory, __module__; compile=false)
    areequal(getfield(__module__, G), t, exprs...; mod=__module__)
end


text:
module PBDS
# __precompile__(false)

using ComputedFieldTypes
using LinearAlgebra, ForwardDiff
using StaticArrays, StaticNumbers, Tullio
using RigidBodyDynamics, MeshCatMechanisms
using GeometryBasics, Rotations
using Base64

import RigidBodyDynamics: CustomCollections, CacheElement, isdirty, setdirty!
import GeometryBasics: min_euclidean, min_euclideansq, min_dist_dim

# types.jl
export TangentBundle, T, ProductManifold, PM, ProductTaskMap, PTM, Chart, ChartList, TaskList,
       TaskGDSList, EmbRep, ChartRep
# manifolds.jl
export dim, embdim, basedim, default_chart, choose_chart_emb, choose_chart_chart, base,
       emb_to_chart, chart_to_emb
# manifolds/
export ℝ, 𝕊, SterProj, SterProjSouth, SterProjNorth, Angleπ, Angle2π
# tasks.jl
export task_map_emb, task_map_emb_chart, task_map_chart_emb, task_map_chart,
       base_task_map_emb, base_task_map_emb_chart, base_task_map_chart_emb, base_task_map_chart,
       task_jacobian_emb, task_jacobian_emb_chart, task_jacobian_chart_emb, task_jacobian_chart,
       base_task_jacobian_emb, base_task_jacobian_emb_chart,
       base_task_jacobian_chart_emb, base_task_jacobian_chart,
       task_jacobian_emb_dot, task_jacobian_chart_dot,
       base_task_jacobian_emb_dot, base_task_jacobian_chart_dot
# tasks/
export AngularDistance, AngularDistanceT,
       AngularPositionAroundPoint, AngularPositionAroundPointT,
       DistanceFromPoint, DistanceFromPointT, DistanceFromPointDynamic,
       DistanceFromSphereSurface, DistanceFromSphereSurfaceT,
       DistanceSphereToSphere, DistanceSphereToSphereT,
       DistanceSphereToCylinder, DistanceSphereToCylinderT,
       DistanceSphereToCup, DistanceSphereToCupT, DistanceSphereToCupDynamic,
       DistanceSphereToBox, DistanceSphereToBoxT, DistanceSphereToBoxDynamic,
       DistanceSphereToLine, DistanceSphereToLineT,
       Identity, IdentityT,
       PositionAroundPoint, PositionAroundPointT,
       Coordinate, CoordinateT,
       XCoordinateDistance, XCoordinateDistanceT,
       YCoordinateDistance, YCoordinateDistanceT,
       LinkFrameTransform, LinkFrameTransformT,
       LinkSpherePosition, LinkSpherePositionT
# tasks/arm/
export @ext_standard_task_type, LinkSpherePairDistance, LinkSpherePairDistanceT,
       JointToLinkPosition, JointToLinkPositionT
# tasks/task_types.jl
export BlankTask, Damping, Attractor, CollisionAvoidance, VelocityLimiter,
       BlankTaskGDS, DampingGDS, AttractorGDS, CollisionAvoidanceGDS, VelocityLimiterGDS
# caches.jl
export ControllerCache
# caches/
export FramePointCache, LinkFrameCache
# trees.jl
export TreeNode, add_child!, check_tree, set_charts_emb!, set_charts_chart!
# forces.jl
export default_potential, default_dissipative_forces,
       potential_from_home_chart, dissipative_forces_from_home_chart
# metrics.jl
export default_metric, metric_from_home_chart, default_weight_metric, weight_metric_from_home_chart
# zero_order_tasks.jl
export single_task_acceleration, multiple_task_acceleration, task_acceleration
# horizontal_bundle_tasks.jl
export dissipative_term_emb, dissipative_term_chart
# trajectory.jl
export propagate_task, propagate_tasks, propagate_task_emb, propagate_tasks_emb
# utils.jl
export smooth_normalization, smooth_norm, exp_barrier_value, exp_barrier_slope, softplus_slope,
       smooth_max, smooth_min, smooth_clamp, smooth_abs, logistic_val, wrapTo2Pi, wrapToPi,
       line_segment_distance, aligned_box_distance, html_video

include("types.jl")
include("manifolds.jl")
include("metrics.jl")
include("connections.jl")
include("caches.jl")
include("tasks.jl")
include("forces.jl")
include("trees.jl")
include("pbds_tasks.jl")
include("pbds_tb_tasks.jl")
include("rmpflow_tasks.jl")
include("trajectory.jl")
include("utils.jl")

end #module


text:
@testset "485.max-consecutive-ones.jl" begin
    @test find_max_consecutive_ones([1, 1, 0, 1, 1, 1]) == 3
end


text:
# # Grid operations

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
There are a variety of (purely) grid-based operators that are useful for carrying
out calculations in immersed layer problems. We will demonstrate a few of them
here.
We will start by generating the cache, just as we did in [Immersed layer caches](@ref)
=#

using ImmersedLayers
using CartesianGrids
using Plots

#=
### Set up a grid and cache
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)

#=
We still generate a cache for these operations, but
now, we only supply the grid. There are no immersed surfaces
for this demonstration.
=#
cache = SurfaceScalarCache(g,scaling=GridScaling)

#=
To demonstrate, let's generate a Gaussian
=#
p = zeros_grid(cache)
xg, yg = x_grid(cache), y_grid(cache)
p .= exp.(-(xg∘xg)-(yg∘yg))

#=
Now, let's generate the gradient of these data
=#
v = zeros_gridgrad(cache)
grad!(v,p,cache)
plot(v,cache)

#=
And finally, let's compute the convective derivative,

$$\mathbf{v}\cdot\nabla\mathbf{v}$$

For this, we create a separate cache, using [`ConvectiveDerivativeCache`](@ref), which
can be constructed from the existing `cache`. This extra cache holds additional
memory for making the calculation of the convective derivative faster. We will
=#
cdcache = ConvectiveDerivativeCache(cache)
vdv = zeros_gridgrad(cache)
convective_derivative!(vdv,v,cache,cdcache) #hide
@time convective_derivative!(vdv,v,cache,cdcache)
nothing #hide

# Plot it
plot(vdv,cache)

#md # ## Surface-grid operator functions
#md # ```@docs
#md # divergence!
#md # grad!
#md # curl!
#md # convective_derivative!
#md # convective_derivative
#md # ConvectiveDerivativeCache
#md # laplacian!
#md # inverse_laplacian!
#md # ```


text:
module host

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Pkg
Pkg.add("Coverage")

using Coverage
Codecov.submit_local(process_folder())


text:
function MixtureMultivariateRealDistribution(arg0::List)
    return MixtureMultivariateRealDistribution((List,), arg0)
end

function MixtureMultivariateRealDistribution(arg0::RandomGenerator, arg1::List)
    return MixtureMultivariateRealDistribution((RandomGenerator, List), arg0, arg1)
end


text:
A = matrixdepot("rosser", Int, 8, 2, 1)

# B is the test matrix used in the paper
#. Rosser, Lanczos, Hestenes and Karush, J. Res. Natl. Bur. Stand. Vol. 47 (1951), pp291-297.
B = [611 196 -192 407 -8 -52 -49 29;
     196 899 113 -192 -71 -43 -8 -44;
     -192 113 899 196 61 49 8 52;
     407 -192 196 611 8 44 59 -23;
     -8 -71 61 8 411 -599 208 208;
     -52 -43 49 44 -599 411 208 208;
     -49 -8 8 59 208 208 99 -911;
     29 -44 52 -23 208 208 -911 99];

@test A == B

# test eigenvalues
e1 = eigvals(matrixdepot("rosser", 2, 2, 1))
e2 = [500, 510]
@test e1 ≈ e2

e1 = eigvals(matrixdepot("rosser", 4, 2, 1))
e2 = [0.1, 1000, 1019.9, 1020]
@test e1 ≈ e2 atol=1e-2

e1 = eigvals(matrixdepot("rosser", 16, 2, 1))
e2 = [-1020, -1020, 0, 0, 0.098, 0.098, 1000, 1000, 1000, 1000, 1000,
      1020, 1020, 1020, 1020, 1020]
@test e1 ≈ e2 atol=11

@test_throws ArgumentError matrixdepot("rosser", 0)
@test matrixdepot("rosser", 1) isa Matrix

println("'rosser' passed test...")


text:
using Base.Iterators: product

# --- Nodes and States ---

"""
    Node = Int16

Primitive type for node index. Alias for `Int16`.
"""
const Node = Int16

"""
    Name = String

Primitive type for node names. Alias for `String`.
"""
const Name = String

"""
    abstract type AbstractNode end

Node type for directed, acyclic graph.
"""
abstract type AbstractNode end

"""
    struct ChanceNode <: AbstractNode

A struct for chance nodes, includes the name, information set and states of the node
"""
struct ChanceNode <: AbstractNode
    name::Name
    I_j::Vector{Name}
    states::Vector{Name}
    function ChanceNode(name, I_j, states)
        return new(name, I_j, states)
    end

end

"""
    struct DecisionNode <: AbstractNode

A struct for decision nodes, includes the name, information set and states of the node
"""
struct DecisionNode <: AbstractNode
    name::Name
    I_j::Vector{Name}
    states::Vector{Name}
    function DecisionNode(name, I_j, states)
        return new(name, I_j, states)
    end
end

"""
    struct ValueNode <: AbstractNode

A struct for value nodes, includes the name and information set of the node
"""
struct ValueNode <: AbstractNode
    name::Name
    I_j::Vector{Name}
    function ValueNode(name, I_j)
        return new(name, I_j)
    end
end

"""
    const State = Int

Primitive type for the number of states. Alias for `Int16`.
"""
const State = Int16

"""
    struct States <: AbstractArray{State, 1}

States type. Works like `Vector{State}`.

# Examples
```julia-repl
julia> S = States(State.([2, 3, 2, 4]))
4-element States:
 2
 3
 2
 4
```
"""
struct States <: AbstractArray{State, 1}
    vals::Vector{State}
    function States(vals::Vector{State})
        if !all(vals .≥ 1)
            throw(DomainError("All states must be ≥ 1."))
        end
        new(vals)
    end
end

Base.size(S::States) = size(S.vals)
Base.IndexStyle(::Type{<:States}) = IndexLinear()
Base.getindex(S::States, i::Int) = getindex(S.vals, i)
Base.length(S::States) = length(S.vals)
Base.eltype(S::States) = eltype(S.vals)

# --- Paths ---

"""
    const Path{N} = NTuple{N, State} where N

Path type. Alias for `NTuple{N, State} where N`.
"""
const Path{N} = NTuple{N, State} where N

"""
    const ForbiddenPath = Tuple{Vector{Node}, Set{Path}}

ForbiddenPath type.

# Examples
```julia-repl
julia> ForbiddenPath(([1, 2], Set([(1, 2)])))
(Int16[1, 2], Set(Tuple{Vararg{Int16,N}} where N[(1, 2)])

julia> ForbiddenPath[
    ([1, 2], Set([(1, 2)])),
    ([3, 4, 5], Set([(1, 2, 3), (3, 4, 5)]))
]
2-element Array{Tuple{Array{Int16,1},Set{Tuple{Vararg{Int16,N}} where N}},1}:
 ([1, 2], Set([(1, 2)]))
 ([3, 4, 5], Set([(1, 2, 3), (3, 4, 5)]))
```
"""
const ForbiddenPath = Tuple{Vector{Node}, Set{Path}}

"""
    const FixedPath = Dict{Node, State}

FixedPath type.

# Examples
```julia-repl
julia> FixedPath(Dict(1=>1, 2=>3))
Dict{Int16,Int16} with 2 entries:
  2 => 3
  1 => 1
```
"""
const FixedPath = Dict{Node, State}

"""
    function paths(states::AbstractVector{State})

Iterate over paths in lexicographical order.

# Examples
```julia-repl
julia> states = States(State.([2, 3]))
2-element States:
 2
 3

julia> vec(collect(paths(states)))
6-element Array{Tuple{Int16,Int16},1}:
 (1, 1)
 (2, 1)
 (1, 2)
 (2, 2)
 (1, 3)
 (2, 3)
```
"""
function paths(states::AbstractVector{State})
    product(UnitRange.(one(eltype(states)), states)...)
end

"""
    function paths(states::AbstractVector{State}, fixed::FixedPath)

Iterate over paths with fixed states in lexicographical order.

# Examples
```julia-repl
julia> states = States(State.([2, 3]))
2-element States:
 2
 3

julia> vec(collect(paths(states, Dict(Node(1) => State(2)))))
3-element Array{Tuple{Int16,Int16},1}:
 (2, 1)
 (2, 2)
 (2, 3)
```
"""
function paths(states::AbstractVector{State}, fixed::FixedPath)
    iters = collect(UnitRange.(one(eltype(states)), states))
    for (i, v) in fixed
        iters[i] = UnitRange(v, v)
    end
    product(iters...)
end

# --- Probabilities ---

"""
    struct Probabilities{N} <: AbstractArray{Float64, N}

Construct and validate stage probabilities (probabilities for a single node).

# Examples
```julia-repl
julia> data = [0.5 0.5 ; 0.2 0.8]
2×2 Array{Float64,2}:
 0.5  0.5
 0.2  0.8

julia> X = Probabilities(Node(2), data)
2×2 Probabilities{2}:
 0.5  0.5
 0.2  0.8

julia> s = (1, 2)
(1, 2)

julia> X(s)
0.5
```
"""
struct Probabilities{N} <: AbstractArray{Float64, N}
    c::Node
    data::Array{Float64, N}
    function Probabilities(c::Node, data::Array{Float64, N}) where N
        for i in CartesianIndices(size(data)[1:end-1])
            if !(sum(data[i, :]) ≈ 1)
                throw(DomainError("Probabilities should sum to one."))
            end
        end
        new{N}(c, data)
    end
end

Base.size(P::Probabilities) = size(P.data)
Base.IndexStyle(::Type{<:Probabilities}) = IndexLinear()
Base.getindex(P::Probabilities, i::Int) = getindex(P.data, i)
Base.getindex(P::Probabilities, I::Vararg{Int,N}) where N = getindex(P.data, I...)

(X::Probabilities)(s::Path) = X[s...]

# --- Path Probability ---

"""
    abstract type AbstractPathProbability end

Abstract path probability type.
"""
abstract type AbstractPathProbability end

"""
    struct DefaultPathProbability <: AbstractPathProbability

Path probability obtained as a product of the probability values corresponding to path s in each chance node.

# Examples
```julia-repl
julia> C = [2]
1-element Array{Int64,1}:
 2

julia> I_j = [[1]]
1-element Array{Array{Int64,1},1}:
 [1]

julia> X = [Probabilities(Node(2), [0.5 0.5; 0.2 0.8])]
1-element Array{Probabilities{2},1}:
 [0.5 0.5; 0.2 0.8]

julia> P = DefaultPathProbability(C, I_j, X)
DefaultPathProbability(Int16[2], Array{Int16,1}[[1]], Probabilities[[0.5 0.5; 0.2 0.8]])

julia> s = Path((1, 2))
(1, 2)

julia> P(s)
0.5
```
"""
struct DefaultPathProbability <: AbstractPathProbability
    C::Vector{Node}
    I_c::Vector{Vector{Node}}
    X::Vector{Probabilities}
    function DefaultPathProbability(C, I_c, X)
        if length(C) == length(I_c)
            new(C, I_c, X)
        else
            throw(DomainError("The number of chance nodes and information sets given to DefaultPathProbability should be equal."))
        end
    end

end

function (P::DefaultPathProbability)(s::Path)
    prod(X(s[[I_c; c]]) for (c, I_c, X) in zip(P.C, P.I_c, P.X))
end

# --- Utilities ---

"""
    const Utility = Float32

Primitive type for utility. Alias for `Float32`.
"""
const Utility = Float32

"""
    struct Utilities{N} <: AbstractArray{Utility, N}

State utilities.

# Examples
```julia-repl
julia> vals = Utility.([1.0 -2.0; 3.0 4.0])
2×2 Array{Float32,2}:
 1.0  -2.0
 3.0   4.0

julia> Y = Utilities(Node(3), vals)
2×2 Utilities{2}:
 1.0  -2.0
 3.0   4.0

julia> s = Path((1, 2))
 (1, 2)

julia> Y(s)
-2.0f0
```
"""
struct Utilities{N} <: AbstractArray{Utility, N}
    v::Node
    data::Array{Utility, N}
    function Utilities(v::Node, data::Array{Utility, N}) where N
        if any(isinf(u) for u in data)
            throw(DomainError("A value should be defined for each element of a utility matrix."))
        end
        new{N}(v, data)
    end
end

Base.size(Y::Utilities) = size(Y.data)
Base.IndexStyle(::Type{<:Utilities}) = IndexLinear()
Base.getindex(Y::Utilities, i::Int) = getindex(Y.data, i)
Base.getindex(Y::Utilities, I::Vararg{Int,N}) where N = getindex(Y.data, I...)

(Y::Utilities)(s::Path) = Y[s...]

# --- Path Utility ---

"""
    abstract type AbstractPathUtility end

Abstract path utility type.
"""
abstract type AbstractPathUtility end

"""
    struct DefaultPathUtility <: AbstractPathUtility

Default path utility obtained as a sum of the utility values corresponding to path s in each value node.

# Examples
```julia-repl
julia> vals = Utility.([1.0 -2.0; 3.0 4.0])
2×2 Array{Float32,2}:
 1.0  -2.0
 3.0   4.0

julia> Y = [Utilities(Node(3), vals)]
1-element Array{Utilities{2},1}:
 [1.0 -2.0; 3.0 4.0]

julia> I_3 = [[1,2]]
1-element Array{Array{Int64,1},1}:
 [1, 2]

julia> U = DefaultPathUtility(I_3, Y)
DefaultPathUtility(Array{Int16,1}[[1, 2]], Utilities[[1.0 -2.0; 3.0 4.0]])

julia> s = Path((1, 2))
(1, 2)

julia> U(s)
-2.0f0

julia> t = Utility(-100.0)

julia> U(s, t)
-102.0f0
```
"""
struct DefaultPathUtility <: AbstractPathUtility
    I_v::Vector{Vector{Node}}
    Y::Vector{Utilities}
end

function (U::DefaultPathUtility)(s::Path)
    sum(Y(s[I_v]) for (I_v, Y) in zip(U.I_v, U.Y))
end

function (U::DefaultPathUtility)(s::Path, t::Utility)
    U(s) + t
end

# --- Influence diagram ---
"""
    mutable struct InfluenceDiagram
        Nodes::Vector{AbstractNode}
        Names::Vector{Name}
        I_j::Vector{Vector{Node}}
        States::Vector{Vector{Name}}
        S::States
        C::Vector{Node}
        D::Vector{Node}
        V::Vector{Node}
        X::Vector{Probabilities}
        Y::Vector{Utilities}
        P::AbstractPathProbability
        U::AbstractPathUtility
        translation::Utility
        function InfluenceDiagram()
            new(Vector{AbstractNode}())
        end
    end

Hold all information related to the influence diagram.

# Fields
- `Nodes::Vector{AbstractNode}`: Vector of added abstract nodes.
- `Names::Vector{Name}`: Names of nodes in order of their indices.
- `I_j::Vector{Vector{Node}}`: Information sets of nodes in order of their indices.
    Nodes of information sets identified by their indices.
- `States::Vector{Vector{Name}}`: States of each node in order of their indices.
- `S::States`: Vector showing the number of states each node has.
- `C::Vector{Node}`: Indices of chance nodes in ascending order.
- `D::Vector{Node}`: Indices of decision nodes in ascending order.
- `V::Vector{Node}`: Indices of value nodes in ascending order.
- `X::Vector{Probabilities}`: Probability matrices of chance nodes in order of chance
    nodes in C.
- `Y::Vector{Utilities}`: Utility matrices of value nodes in order of value nodes in V.
- `P::AbstractPathProbability`: Path probabilities.
- `U::AbstractPathUtility`: Path utilities.
- `translation::Utility`: Utility translation for storing the positive or negative
    utility translation.

# Examples
```julia
diagram = InfluenceDiagram()
```
"""
mutable struct InfluenceDiagram
    Nodes::Vector{AbstractNode}
    Names::Vector{Name}
    I_j::Vector{Vector{Node}}
    States::Vector{Vector{Name}}
    S::States
    C::Vector{Node}
    D::Vector{Node}
    V::Vector{Node}
    X::Vector{Probabilities}
    Y::Vector{Utilities}
    P::AbstractPathProbability
    U::AbstractPathUtility
    translation::Utility
    function InfluenceDiagram()
        new(Vector{AbstractNode}())
    end
end

# --- Adding nodes ---

function validate_node(diagram::InfluenceDiagram,
    name::Name,
    I_j::Vector{Name};
    value_node::Bool=false,
    states::Vector{Name}=Vector{Name}())

    if !allunique([map(x -> x.name, diagram.Nodes)..., name])
        throw(DomainError("All node names should be unique."))
    end

    if !allunique(I_j)
        throw(DomainError("All nodes in an information set should be unique."))
    end

    if !allunique([name, I_j...])
        throw(DomainError("Node should not be included in its own information set."))
    end

    if !value_node
        if length(states) < 2
            throw(DomainError("Each chance and decision node should have more than one state."))
        end
    end

    if value_node
        if isempty(I_j)
            @warn("Value node $name is redundant.")
        end
    end
end

"""
    function add_node!(diagram::InfluenceDiagram, node::AbstractNode)

Add node to influence diagram structure.

# Examples
```julia-repl
julia> add_node!(diagram, ChanceNode("O", [], ["lemon", "peach"]))
1-element Array{AbstractNode,1}:
 ChanceNode("O", String[], ["lemon", "peach"])
```
"""
function add_node!(diagram::InfluenceDiagram, node::AbstractNode)
    if !isa(node, ValueNode)
        validate_node(diagram, node.name, node.I_j, states = node.states)
    else
        validate_node(diagram, node.name, node.I_j, value_node = true)
    end
    push!(diagram.Nodes, node)
end

# --- Adding Probabilities ---
"""
    struct ProbabilityMatrix{N} <: AbstractArray{Float64, N}
        nodes::Vector{Name}
        indices::Vector{Dict{Name, Int}}
        matrix::Array{Float64, N}
    end

Construct probability matrix.
"""
struct ProbabilityMatrix{N} <: AbstractArray{Float64, N}
    nodes::Vector{Name}
    indices::Vector{Dict{Name, Int}}
    matrix::Array{Float64, N}
end

Base.size(PM::ProbabilityMatrix) = size(PM.matrix)
Base.getindex(PM::ProbabilityMatrix, I::Vararg{Int,N}) where N = getindex(PM.matrix, I...)
function Base.setindex!(PM::ProbabilityMatrix, p::T, I::Vararg{Union{String, Int},N}) where {N, T<:Real}
    I2 = []
    for i in 1:N
        if isa(I[i], String)
            if get(PM.indices[i], I[i], 0) == 0
                throw(DomainError("Node $(probability_matrix.nodes[i]) does not have state $(I[i])."))
            end
            push!(I2, PM.indices[i][I[i]])
        else
            push!(I2, I[i])
        end
    end
    PM.matrix[I2...] = p
end
function Base.setindex!(PM::ProbabilityMatrix{N}, P::Array{T}, I::Vararg{Union{String, Int, Colon}, N}) where {N, T<:Real}
    I2 = []
    for i in 1:N
        if isa(I[i], Colon)
            push!(I2, :)
        elseif isa(I[i], String)
            if get(PM.indices[i], I[i], 0) == 0
                throw(DomainError("Node $(probability_matrix.nodes[i]) does not have state $(I[i])."))
            end
            push!(I2, PM.indices[i][I[i]])
        else
            push!(I2, I[i])
        end
    end
    PM.matrix[I2...] = P
end

"""
    function ProbabilityMatrix(diagram::InfluenceDiagram, node::Name)

Initialise a probability matrix for a given chance node. The matrix is initialised with zeros.

# Examples
```julia-repl
julia> X_O = ProbabilityMatrix(diagram, "O")
2-element ProbabilityMatrix{1}:
 0.0
 0.0
```
"""
function ProbabilityMatrix(diagram::InfluenceDiagram, node::Name)
    if node ∉ diagram.Names
        throw(DomainError("Node $node should be added as a node to the influence diagram."))
    end
    if node ∉ diagram.Names[diagram.C]
        throw(DomainError("Only chance nodes can have probability matrices."))
    end

    # Find the node's indices and it's I_c nodes
    c = findfirst(x -> x==node, diagram.Names)
    nodes = [diagram.I_j[c]..., c]
    names = diagram.Names[nodes]

    indices = Vector{Dict{Name, Int}}()
    for j in nodes
        states = Dict{Name, Int}(state => i
            for (i, state) in enumerate(diagram.States[j])
        )
        push!(indices, states)
    end
    matrix = fill(0.0, diagram.S[nodes]...)

    return ProbabilityMatrix(names, indices, matrix)
end

"""
    function add_probabilities!(diagram::InfluenceDiagram, node::Name, probabilities::AbstractArray{Float64, N}) where N

Add probability matrix to influence diagram, specifically to its `X` vector.

# Examples
```julia
julia> X_O = ProbabilityMatrix(diagram, "O")
2-element ProbabilityMatrix{1}:
 0.0
 0.0

julia> X_O["lemon"] = 0.2
0.2

julia> add_probabilities!(diagram, "O", X_O)
ERROR: DomainError with Probabilities should sum to one.:

julia> X_O["peach"] = 0.8
0.2

julia> add_probabilities!(diagram, "O", X_O)
1-element Array{Probabilities,1}:
 [0.2, 0.8]
```
!!! note
    The function `generate_arcs!` must be called before probabilities or utilities can be added to the influence diagram.
"""
function add_probabilities!(diagram::InfluenceDiagram, node::Name, probabilities::AbstractArray{Float64, N}) where N
    c = findfirst(x -> x==node, diagram.Names)

    if c ∈ [j.c for j in diagram.X]
        throw(DomainError("Probabilities should be added only once for each node."))
    end

    if size(probabilities) == Tuple((diagram.S[j] for j in (diagram.I_j[c]..., c)))
        if isa(probabilities, ProbabilityMatrix)
            # Check that probabilities sum to one happesn in Probabilities
            push!(diagram.X, Probabilities(Node(c), probabilities.matrix))
        else
            push!(diagram.X, Probabilities(Node(c), probabilities))
        end
    else
        throw(DomainError("The dimensions of a probability matrix should match the node's states' and information states' cardinality. Expected $(Tuple((diagram.S[n] for n in (diagram.I_j[c]..., c)))) for node $name, got $(size(probabilities))."))
    end
end

# --- Adding Utilities ---

"""
    struct UtilityMatrix{N} <: AbstractArray{Utility, N}
        I_v::Vector{Name}
        indices::Vector{Dict{Name, Int}}
        matrix::Array{Utility, N}
    end

Construct utility matrix.
"""
struct UtilityMatrix{N} <: AbstractArray{Utility, N}
    I_v::Vector{Name}
    indices::Vector{Dict{Name, Int}}
    matrix::Array{Utility, N}
end

Base.size(UM::UtilityMatrix) = size(UM.matrix)
Base.getindex(UM::UtilityMatrix, I::Vararg{Int,N}) where N = getindex(UM.matrix, I...)
function Base.setindex!(UM::UtilityMatrix{N}, y::T, I::Vararg{Union{String, Int},N}) where {N, T<:Real}
    I2 = []
    for i in 1:N
        if isa(I[i], String)
            if get(UM.indices[i], I[i], 0) == 0
                throw(DomainError("Node $(probability_matrix.nodes[i]) does not have state $(I[i])."))
            end
            push!(I2, UM.indices[i][I[i]])
        else
            push!(I2, I[i])
        end
    end
    UM.matrix[I2...] = y
end
function Base.setindex!(UM::UtilityMatrix{N}, Y::Array{T}, I::Vararg{Union{String, Int, Colon}, N}) where {N, T<:Real}
    I2 = []
    for i in 1:N
        if isa(I[i], Colon)
            push!(I2, :)
        elseif isa(I[i], String)
            if get(UM.indices[i], I[i], 0) == 0
                throw(DomainError("Node $(probability_matrix.nodes[i]) does not have state $(I[i])."))
            end
            push!(I2, UM.indices[i][I[i]])
        else
            push!(I2, I[i])
        end
    end
    UM.matrix[I2...] = Y
end

"""
    function UtilityMatrix(diagram::InfluenceDiagram, node::Name)

Initialise a utility matrix for a value node. The matrix is initialised with `Inf` values.

# Examples
```julia-repl
julia> Y_V3 = UtilityMatrix(diagram, "V3")
2×3 UtilityMatrix{2}:
 Inf  Inf  Inf
 Inf  Inf  Inf
```
"""
function UtilityMatrix(diagram::InfluenceDiagram, node::Name)
    if node ∉ diagram.Names
        throw(DomainError("Node $node should be added as a node to the influence diagram."))
    end
    if node ∉ diagram.Names[diagram.V]
        throw(DomainError("Only value nodes can have consequence matrices."))
    end

    # Find the node's indexand it's I_v nodes
    v = findfirst(x -> x==node, diagram.Names)
    I_v = diagram.I_j[v]
    names = diagram.Names[I_v]

    indices = Vector{Dict{Name, Int}}()
    for j in I_v
        states = Dict{Name, Int}(state => i
            for (i, state) in enumerate(diagram.States[j])
        )
        push!(indices, states)
    end
    matrix = Array{Utility}(fill(Inf, diagram.S[I_v]...))

    return UtilityMatrix(names, indices, matrix)
end

"""
    function add_utilities!(diagram::InfluenceDiagram, node::Name, utilities::AbstractArray{T, N}) where {N,T<:Real}

Add utility matrix to influence diagram, specifically to its `Y` vector.

# Examples
```julia-repl
julia> Y_V3 = UtilityMatrix(diagram, "V3")
2×3 UtilityMatrix{2}:
 Inf  Inf  Inf
 Inf  Inf  Inf

julia> Y_V3["peach", :] = [-40, -20, 0]
3-element Array{Int64,1}:
 -40
 -20
   0

julia> Y_V3["lemon", :] = [-200, 0, 0]
3-element Array{Int64,1}:
 -200
    0
    0

julia> add_utilities!(diagram, "V3", Y_V3)
1-element Array{Utilities,1}:
 [-200.0 0.0 0.0; -40.0 -20.0 0.0]

julia> add_utilities!(diagram, "V1", [0, -25])
2-element Array{Utilities,1}:
 [-200.0 0.0 0.0; -40.0 -20.0 0.0]
 [0.0, -25.0]
```
!!! note
    The function `generate_arcs!` must be called before probabilities or utilities can be added to the influence diagram.
"""
function add_utilities!(diagram::InfluenceDiagram, node::Name, utilities::AbstractArray{T, N}) where {N,T<:Real}
    v = findfirst(x -> x==node, diagram.Names)

    if v ∈ [j.v for j in diagram.Y]
        throw(DomainError("Utilities should be added only once for each node."))
    end
    if any(u ==Inf for u in utilities)
        throw(DomainError("Utility values should be less than infinity."))
    end

    if size(utilities) == Tuple((diagram.S[j] for j in diagram.I_j[v]))
        if isa(utilities, UtilityMatrix)
            push!(diagram.Y, Utilities(Node(v), utilities.matrix))
        else
            # Conversion to Float32 using Utility(), since machine default is Float64
            push!(diagram.Y, Utilities(Node(v), [Utility(u) for u in utilities]))
        end
    else
        throw(DomainError("The dimensions of the utilities matrix should match the node's information states' cardinality. Expected $(Tuple((diagram.S[n] for n in diagram.I_j[v]))) for node $name, got $(size(utilities))."))
    end
end

# --- Generating Arcs ---

function validate_structure(Nodes::Vector{AbstractNode}, C_and_D::Vector{AbstractNode}, n_CD::Int, V::Vector{AbstractNode}, n_V::Int)
    # Validating node structure
    if n_CD == 0
        throw(DomainError("The influence diagram must have chance or decision nodes."))
    end
    if !(union((n.I_j for n in Nodes)...) ⊆ Set(n.name for n in Nodes))
        throw(DomainError("Each node that is part of an information set should be added as a node."))
    end
    # Checking the information sets of C and D nodes
    if !isempty(union((j.I_j for j in C_and_D)...) ∩ Set(v.name for v in V))
        throw(DomainError("Information sets should not include any value nodes."))
    end
    # Checking the information sets of V nodes
    if !isempty(V) && !isempty(union((v.I_j for v in V)...) ∩ Set(v.name for v in V))
        throw(DomainError("Information sets should not include any value nodes."))
    end
    # Check for redundant chance or decision nodes.
    last_CD_nodes = setdiff((j.name for j in C_and_D), (j.I_j for j in C_and_D)...)
    for i in last_CD_nodes
        if !isempty(V) && i ∉ union((v.I_j for v in V)...)
            @warn("Node $i is redundant.")
        end
    end
end

"""
    function generate_arcs!(diagram::InfluenceDiagram)

Generate arc structures using nodes added to influence diagram, by ordering nodes,
giving them indices and generating correct values for the vectors Names, I_j, states,
S, C, D, V in the influence digram. Abstraction is created and the names of the nodes
and states are only used in the user interface from here on.

# Examples
```julia
generate_arcs!(diagram)
```
"""
function generate_arcs!(diagram::InfluenceDiagram)

    # Chance and decision nodes
    C_and_D = filter(x -> !isa(x, ValueNode), diagram.Nodes)
    n_CD = length(C_and_D)
    # Value nodes
    V_nodes = filter(x -> isa(x, ValueNode), diagram.Nodes)
    n_V = length(V_nodes)

    validate_structure(diagram.Nodes, C_and_D, n_CD, V_nodes, n_V)

    # Declare vectors for results (final resting place InfluenceDiagram.Names and InfluenceDiagram.I_j)
    Names = Vector{Name}(undef, n_CD+n_V)
    I_j = Vector{Vector{Node}}(undef, n_CD+n_V)
    states = Vector{Vector{Name}}()
    S = Vector{State}(undef, n_CD)
    C = Vector{Node}()
    D = Vector{Node}()
    V = Vector{Node}()

    # Declare helper collections
    indices = Dict{Name, Node}()
    indexed_nodes = Set{Name}()
    # Declare index
    index = 1

    while true
        # Index nodes C and D that don't yet have indices but whose I_j have indices
        new_nodes = filter(j -> (j.name ∉ indexed_nodes && Set(j.I_j) ⊆ indexed_nodes), C_and_D)
        for j in new_nodes
            # Update helper collections
            push!(indices, j.name => index)
            push!(indexed_nodes, j.name)
            # Update results
            Names[index] = Name(j.name)    #TODO datatype conversion happens here, should we use push! ?
            I_j[index] = map(x -> Node(indices[x]), j.I_j)
            push!(states, j.states)
            S[index] = State(length(j.states))
            if isa(j, ChanceNode)
                push!(C, Node(index))
            else
                push!(D, Node(index))
            end
            # Increase index
            index += 1
        end

        # If no new nodes were indexed this iteration, terminate while loop
        if isempty(new_nodes)
            if index < n_CD
                throw(DomainError("The influence diagram should be acyclic."))
            else
                break
            end
        end
    end

    # Index value nodes
    for v in V_nodes
        # Update results
        Names[index] = Name(v.name)
        I_j[index] = map(x -> Node(indices[x]), v.I_j)
        push!(V, Node(index))
        # Increase index
        index += 1
    end

    diagram.Names = Names
    diagram.I_j = I_j
    diagram.States = states
    diagram.S = States(S)
    diagram.C = C
    diagram.D = D
    diagram.V = V
    # Declaring X and Y
    diagram.X = Vector{Probabilities}()
    diagram.Y = Vector{Utilities}()
end

# --- Generating Diagram ---
"""
    function generate_diagram!(diagram::InfluenceDiagram;
    default_probability::Bool=true,
    default_utility::Bool=true,
    positive_path_utility::Bool=false,
    negative_path_utility::Bool=false)

Generate complete influence diagram with probabilities and utilities as well.

# Arguments
- `default_probability::Bool=true`: Choice to use default path probabilities.
- `default_utility::Bool=true`: Choice to use default path utilities.
- `positive_path_utility::Bool=false`: Choice to use a positive path utility translation.
- `negative_path_utility::Bool=false`: Choice to use a negative path utility translation.

# Examples
```julia
generate_diagram!(diagram)
```

!!! note
    The influence diagram must be generated after probabilities and utilities are added
    but before creating the decision model.

!!! note
    If the default probabilities and utilities are not used, define `AbstractPathProbability`
    and `AbstractPathUtility` structures and define P(s), U(s) and U(s, t) functions
    for them. Add the `AbstractPathProbability` and `AbstractPathUtility` structures
    to the influence diagram fields P and U.
"""
function generate_diagram!(diagram::InfluenceDiagram;
    default_probability::Bool=true,
    default_utility::Bool=true,
    positive_path_utility::Bool=false,
    negative_path_utility::Bool=false)

    # Sort probabilities and consequences
    sort!(diagram.X, by = x -> x.c)
    sort!(diagram.Y, by = x -> x.v)

    # Declare P and U if defaults are used
    if default_probability
        diagram.P = DefaultPathProbability(diagram.C, diagram.I_j[diagram.C], diagram.X)
    end
    if default_utility
        diagram.U = DefaultPathUtility(diagram.I_j[diagram.V], diagram.Y)
        if positive_path_utility
            # Conversion to Float32 using Utility(), since machine default is Float64
            diagram.translation = 1 -  minimum(diagram.U(s) for s in paths(diagram.S))
        elseif negative_path_utility
            diagram.translation = -1 - maximum(diagram.U(s) for s in paths(diagram.S))
        else
            diagram.translation = 0
        end
    end

end

"""
    function index_of(diagram::InfluenceDiagram, node::Name)

Get the index of a given node.

# Example
```julia-repl
julia> idx_O = index_of(diagram, "O")
1
```
"""
function index_of(diagram::InfluenceDiagram, node::Name)
    idx = findfirst(isequal(node), diagram.Names)
    if isnothing(idx)
        throw(DomainError("Name $node not found in the diagram."))
    end
    return idx
end

"""
    function num_states(diagram::InfluenceDiagram, node::Name)

Get the number of states in a given node.

# Example
```julia-repl
julia> NS_O = num_states(diagram, "O")
2
```
"""
function num_states(diagram::InfluenceDiagram, node::Name)
    return diagram.S[index_of(diagram, node)]
end

# --- ForbiddenPath and FixedPath outer construction functions ---
"""
    function ForbiddenPath(diagram::InfluenceDiagram, nodes::Vector{Name}, paths::Vector{NTuple{N, Name}}) where N

ForbiddenPath outer construction function. Create ForbiddenPath variable.

# Arguments
- `diagram::InfluenceDiagram`: Influence diagram structure
- `nodes::Vector{Name}`: Vector of nodes involved in forbidden paths. Identified by their names.
- `paths`::Vector{NTuple{N, Name}}`: Vector of tuples defining the forbidden combinations of states. States identified by their names.

# Example
```julia
ForbiddenPath(diagram, ["R1", "R2"], [("high", "low"), ("low", "high")])
```
"""
function ForbiddenPath(diagram::InfluenceDiagram, nodes::Vector{Name}, paths::Vector{NTuple{N, Name}}) where N
    node_indices = Vector{Node}()
    for node in nodes
        j = findfirst(i -> i == node, diagram.Names)
        if isnothing(j)
            throw(DomainError("Node $node does not exist."))
        end
        push!(node_indices, j)
    end

    path_set = Set{Path}()
    for s in paths
        s_states = Vector{State}()
        for (i, s_i) in enumerate(s)
            s_i_index = findfirst(x -> x == s_i, diagram.States[node_indices[i]])
            if isnothing(s_i_index)
                throw(DomainError("Node $(nodes[i]) does not have a state called $s_i."))
            end

            push!(s_states, s_i_index)
        end
        push!(path_set, Path(s_states))
    end

    return ForbiddenPath((node_indices, path_set))
end

"""
    function FixedPath(diagram::InfluenceDiagram, fixed::Dict{Name, Name})

FixedPath outer construction function. Create FixedPath variable.

# Arguments
- `diagram::InfluenceDiagram`: Influence diagram structure
- `fixed::Dict{Name, Name}`: Dictionary of nodes and their fixed states. Order is node=>state, and both are idefied with their names.

# Example
```julia-repl
julia> fixedpath = FixedPath(diagram, Dict("O" => "lemon"))
Dict{Int16,Int16} with 1 entry:
  1 => 1

julia> vec(collect(paths(states, fixedpath)))
3-element Array{Tuple{Int16,Int16},1}:
 (1, 1)
 (1, 2)
 (1, 3)

```
"""
function FixedPath(diagram::InfluenceDiagram, fixed::Dict{Name, Name})
    fixed_paths = Dict{Node, State}()

    for (j, s_j) in fixed
        j_index = findfirst(i -> i == j, diagram.Names)
        if isnothing(j_index)
            throw(DomainError("Node $j does not exist."))
        end

        s_j_index = findfirst(s -> s == s_j, diagram.States[j_index])
        if isnothing(s_j_index)
            throw(DomainError("Node $j does not have a state called $s_j."))
        end
        push!(fixed_paths, Node(j_index) => State(s_j_index))
    end

    return FixedPath(fixed_paths)
end

# --- Local Decision Strategy ---

"""
    LocalDecisionStrategy{N} <: AbstractArray{Int, N}

Local decision strategy type.
"""
struct LocalDecisionStrategy{N} <: AbstractArray{Int, N}
    d::Node
    data::Array{Int, N}
    function LocalDecisionStrategy(d::Node, data::Array{Int, N}) where N
        if !all(0 ≤ x ≤ 1 for x in data)
            throw(DomainError("All values x must be 0 ≤ x ≤ 1."))
        end
        for s_I in CartesianIndices(size(data)[1:end-1])
            if !(sum(data[s_I, :]) == 1)
                throw(DomainError("Values should add to one."))
            end
        end
        new{N}(d, data)
    end
end

Base.size(Z::LocalDecisionStrategy) = size(Z.data)
Base.IndexStyle(::Type{<:LocalDecisionStrategy}) = IndexLinear()
Base.getindex(Z::LocalDecisionStrategy, i::Int) = getindex(Z.data, i)
Base.getindex(Z::LocalDecisionStrategy, I::Vararg{Int,N}) where N = getindex(Z.data, I...)

function (Z::LocalDecisionStrategy)(s_I::Path)::State
    findmax(Z[s_I..., :])[2]
end

# --- Decision Strategy ---

"""
    DecisionStrategy

Decision strategy type.
"""
struct DecisionStrategy
    D::Vector{Node}
    I_d::Vector{Vector{Node}}
    Z_d::Vector{LocalDecisionStrategy}
end


text:
"""
    is_leap_year(year)

Return if `year` is a leap year in the gregorian calendar.
"""
isleapyear(y) = (y % 4 == 0) && ((y % 100 != 0) || (y % 400 == 0))


text:
""" `numerical_deriv( x, y)`
Calculate simple estimate of numerical derivative
"""
function numerical_deriv( x::AA1, y::AA2 )   where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}  }
	@assert length(x) == length(y)
	dfluxdlnλ = zeros(size(y))
	dfluxdlnλ[1] = (y[2]-y[1])/(x[2]-x[1])
	dfluxdlnλ[2:end-1] .= (y[3:end].-y[1:end-2])./(x[3:end].-x[1:end-2])
	dfluxdlnλ[end] = (y[end]-y[end-1])/(x[end]-x[end-1])
	return dfluxdlnλ
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    dfdlogλ = Array{T1,1}(undef,length(flux))
    calc_dfluxdlnlambda!(dfdlogλ,flux,λ)
    #dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    #dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    #dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda!(dfdlogλ::AbstractArray{T1,1}, flux::AbstractArray{T2,1}, λ::AbstractArray{T3,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(dfdlogλ) == size(flux)
    @assert length(flux) >= 3
    #dfdlogλ = Array{T1,1}(undef,length(flux))
    dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    logλ = log.(λ)
    d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    #d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2])./(logλ[3:end].-logλ[1:end-2])).^2
    #d2fdlogλ2[end] = d2fdlogλ2[end-1]
    #d2fdlogλ2[1] = d2fdlogλ2[2]
    calc_d2fluxdlnlambda2!(d2fdlogλ2,flux,λ)
    return d2fdlogλ2
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2!(d2fdlogλ2::AbstractArray{T1,1}, flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(d2fdlogλ2) == size(flux)
    @assert length(flux) >= 3
    #logλ = log.(λ)
    #d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2]).^2 ./ (2.0.*log.(λ[3:end]./λ[1:end-2])))
    d2fdlogλ2[end] = d2fdlogλ2[end-1]
    d2fdlogλ2[1] = d2fdlogλ2[2]
    return d2fdlogλ2
end

function calc_weights_for_running_mean(n::Integer; decay_factor::Float64 = 2.0)
	@assert 1<=n<=100
	@assert 1 <= decay_factor < Inf
	w = zeros(n)
	if isodd(n)
		center = convert(Int,1 + (n-1)//2)
		w[center] = 1
		for i in 1:(center-1)
			w[center+i] = w[center-i] = decay_factor^(-i)
		end
	else
		center = convert(Int,n//2)
		for i in 1:center
			w[center+i] = w[center+1-i] = decay_factor^(-(i-0.5))
		end
	end
	w ./= sum(w)
	return w
end

""" Smooth array via weighted running mean of n consecutive points """
function calc_running_mean_weighted(x::AbstractVector{T1}; n::Integer = 5, decay_factor::Float64 = 2.0,
				weights::AbstractVector{T2} = calc_weights_for_running_mean(n, decay_factor=decay_factor) ) where { T1<:Real, T2<:Real }
	@assert 2 <= n < length(x)/2
	@assert isodd(n)   # Could generalize, but haven't
	smoothed = Array{T1,1}(undef,length(x))
	n = length(weights)
	m = floor(Int,n//2)
	for i in 1:m
		start_idx = 2+m-i
		smoothed[i] = sum(view(x,1:i+m).*view(weights,start_idx:n)) /sum(view(weights,start_idx:n))
	end
	for i in (m+1):(length(x)-m-1)
		smoothed[i] = sum(view(x,i-m:i+m) .* weights)
	end
	for i in (length(x)-m):length(x)
		stop_idx = length(x)-(i-m)+1
		smoothed[i] = sum(view(x,i-m:length(x)).*view(weights,1:stop_idx)) /sum(view(weights,1:stop_idx))
	end
	return smoothed
end


text:
CRCBS.extract_feature(solver::NBSSolver,::LowLevelIterationCount, args...) = max_iterations(route_planner(solver))

export
    ObjectPathSummaries

struct ObjectPathSummaries <: FeatureExtractor{Dict{Int,Vector{}}} end
function get_object_path_summaries end
function CRCBS.extract_feature(solver,::ObjectPathSummaries,pc_mapf,env,timer_results)
    summaries = get_object_path_summaries(env)
    to_string_dict(to_string_dict(Dict(k=>to_string_dict(v) for (k,v) in summaries)))
end

extract_object_data(node,t0,tF) = Dict{Symbol,Union{Vector{Int},Int}}()
extract_object_data(node::OBJECT_AT,t0,tF) = Dict{Symbol,Union{Vector{Int},Int}}(
    :object_id=>get_id(get_object_id(node)),
    :start_vtx=>get_id(get_initial_location_id(node)),
    :start_time=>t0,
    )
extract_object_data(node::COLLECT,t0,tF) = Dict{Symbol,Union{Vector{Int},Int}}(
    :object_id=>get_id(get_object_id(node)),
    :start_vtx=>get_id(get_initial_location_id(node)),
    :robot_ids=>map(get_id,get_valid_robot_ids(node)),
    :collect_time=>tF,
    )
extract_object_data(node::CARRY,t0,tF) = Dict{Symbol,Union{Vector{Int},Int}}(
    :object_id=>get_id(get_object_id(node)),
    :start_vtx=>get_id(get_initial_location_id(node)),
    :end_vtx=>get_id(get_destination_location_id(node)),
    :robot_ids=>map(get_id,get_valid_robot_ids(node)),
    :collect_time=>t0,
    )
extract_object_data(node::DEPOSIT,t0,tF) = Dict{Symbol,Union{Vector{Int},Int}}(
    :object_id=>get_id(get_object_id(node)),
    :end_vtx=>get_id(get_destination_location_id(node)),
    :robot_ids=>map(get_id,get_valid_robot_ids(node)),
    :deposit_time=>tF,
    )
function get_object_path_summaries(env::SearchEnv)
    summaries = Dict{Int,Dict{Symbol,Union{Vector{Int},Int}}}()
    for v in vertices(get_schedule(env))
        node = get_node_from_vtx(get_schedule(env),v)
        dict = extract_object_data(node,get_t0(env,v),get_tF(env,v))
        if haskey(dict,:object_id)
            object_id = dict[:object_id]
            merge!(get!(summaries,object_id,valtype(summaries)()),dict)
        end
    end
    summaries
end

function reconstruct_object_paths(robot_paths,object_path_summaries)
    object_paths = Dict{Int,Vector{Int}}()
    object_intervals = Dict{Int,Tuple{Int,Int}}()
    for object_summary in values(object_path_summaries)
        id = object_summary[:object_id]
        t0 = object_summary[:start_time]
        tC = object_summary[:collect_time]
        tF = object_summary[:deposit_time]
        v0 = object_summary[:start_vtx]
        vF = object_summary[:end_vtx]
        robot_ids = object_summary[:robot_ids]
        if robot_ids[1] != -1
            T = length(robot_paths[robot_ids[1]])
            if T >= tC-1
                p1 = map(t->v0, 0:tC-1)
                p2 = robot_paths[robot_ids[1]][tC+1:min(tF+1,T)]
                object_paths[id] = vcat(p1,p2)
                object_intervals[id] = (t0,tF)
            end
        end
    end
    return object_paths, object_intervals
end

export SolutionAdjacencyMatrix
struct SolutionAdjacencyMatrix <: FeatureExtractor{Dict{String,Any}} end
function CRCBS.extract_feature(solver,::SolutionAdjacencyMatrix,pc_tapf,env,timer_results)
    mat = adjacency_matrix(get_graph(get_schedule(env)))
    I,J,_ = findnz(mat)
    return Dict("edge_list" => map(idx->[idx...], zip(I,J)),
        "size_i" => size(mat,1),
        "size_j" => size(mat,2)
        )
end

export
    TaskGraphsProblemLoader,
    add_env!,
    get_env!,
    problem_type

"""
    TaskGraphsProblemLoader{T}

Helper cache for loading problems of type `T`
"""
struct TaskGraphsProblemLoader{T}
    envs::Dict{String,GridFactoryEnvironment}
    prob_specs::Dict{String,ProblemSpec}
    TaskGraphsProblemLoader{T}() where {T} = new{T}(
        Dict{String,GridFactoryEnvironment}(),
        Dict{String,ProblemSpec}()
    )
end
function add_env!(loader::TaskGraphsProblemLoader,env_id::String,
    env=read_env(env_id))
    if haskey(loader.envs,env_id)
    else
        loader.envs[env_id] = env
        # loader.prob_specs[env_id] = ProblemSpec(D=get_dist_matrix(env))
        loader.prob_specs[env_id] = ProblemSpec(D=env)
    end
    loader
end
function get_env!(loader::TaskGraphsProblemLoader,env_id::String)
    if !haskey(loader.envs,env_id)
        loader.envs[env_id] = read_env(env_id)
    end
    return loader.envs[env_id]
end
problem_type(::TaskGraphsProblemLoader{T}) where {T} = T

export
    PCTAPF_Loader,
    PCTA_Loader,
    PCMAPF_Loader,
    ReplanningProblemLoader

const PCTAPF_Loader = TaskGraphsProblemLoader{PC_TAPF}
const PCMAPF_Loader = TaskGraphsProblemLoader{PC_MAPF}
const PCTA_Loader   = TaskGraphsProblemLoader{PC_TA}
const ReplanningProblemLoader = TaskGraphsProblemLoader{RepeatedPC_TAPF}

export write_problem

"""
    write_problem(loader::TaskGraphsProblemLoader,problem_def,prob_path,env_id="")

Write a problem that can later be loaded and solved.
"""
function write_problem(loader::TaskGraphsProblemLoader,problem_def,prob_path,env_id="")
    if isdir(prob_path)
        prob_path = joinpath(prob_path,"problem.toml")
    end
    toml_dict = TOML.parse(problem_def)
    toml_dict["env_id"] = env_id
    open(prob_path, "w") do io
        TOML.print(io, toml_dict)
    end
    prob_path
end

init_random_problem(loader::TaskGraphsProblemLoader,env,config) = random_pctapf_def(env,config)
function write_problems!(loader::TaskGraphsProblemLoader,config::Dict,base_path::String,prob_counter::Int=1)
    env = get_env!(loader,config[:env_id])
    num_trials = get(config,:num_trials,1)
    for i in 1:num_trials
        prob_name = padded_problem_name(prob_counter,"problem","")
        prob_path = joinpath(base_path,prob_name)
        mkpath(prob_path)
        open(joinpath(prob_path,"config.toml"),"w") do io
            TOML.print(io,config)
        end
        prob_def = init_random_problem(loader,env,config)
        write_problem(loader,prob_def,prob_path,config[:env_id])
        prob_counter += 1
    end
    return prob_counter
end
function write_problems!(loader::TaskGraphsProblemLoader,configs::Vector{D},base_path::String,prob_counter::Int=1) where {D<:Dict}
    for config in configs
        prob_counter = write_problems!(loader,config,base_path,prob_counter)
    end
    return prob_counter
end

"""
    CRCBS.load_problem(loader::TaskGraphsProblemLoader,solver_config,prob_path)

Currently only impemented for PC_TAPF and PC_TA
"""
function CRCBS.load_problem(loader::TaskGraphsProblemLoader,solver_config,prob_path)
    if !isfile(prob_path)
        prob_path = joinpath(prob_path,"problem.toml")
    end
    toml_dict = TOML.parsefile(prob_path)
    def = read_problem_def(toml_dict)
    env_id = toml_dict["env_id"]
    env = get_env!(loader,env_id)
    problem_type(loader)(solver_config.solver,def,env)
end
function CRCBS.load_problem(loader::PCTA_Loader,solver_config,prob_path)
    if !isfile(prob_path)
        prob_path = joinpath(prob_path,"problem.toml")
    end
    toml_dict = TOML.parsefile(prob_path)
    def = read_problem_def(toml_dict)
    env_id = toml_dict["env_id"]
    env = get_env!(loader,env_id)
    problem_type(loader)(def,env,solver_config.objective)
end
function is_problem_file(loader::PCTA_Loader,path)
    isfile(joinpath(path,"problem.toml")) && isfile(joinpath(path,"config.toml"))
end
function is_problem_file(loader::PCTAPF_Loader,path)
    isfile(joinpath(path,"problem.toml")) && isfile(joinpath(path,"config.toml"))
end

CRCBS.extract_feature(solver,::SolutionCost,  mapf,solution::OperatingSchedule,timer_results) = best_cost(solver)
"""
    CRCBS.run_profiling(loader::TaskGraphsProblemLoader,solver_config,problem_dir)

Run profiling with a `TaskGraphsProblemLoader`.
"""
function CRCBS.run_profiling(loader::TaskGraphsProblemLoader,solver_config,problem_dir)
    solver = solver_config.solver
    for prob_path in readdir(problem_dir;join=true)
        is_problem_file(loader,prob_path) ? nothing : continue
        prob = load_problem(loader,solver_config,prob_path)
        solution, timer_results = profile_solver!(solver,prob)
        results_dict = compile_results(
            solver,
            solver_config.feats,
            prob,
            solution,
            timer_results
            )
        results_dict["problem_file"] = prob_path
        write_results(loader,solver_config,prob,prob_path,results_dict)
    end
end

export warmup

"""
    warmup(loader::TaskGraphsProblemLoader,solver_config,problem_dir,dummy_path = "dummy_path")

Do a small dry run of `run_profiling(loader,solver_config,problem_dir)` to
ensure that all code is fully compiled before collecting results.
"""
function warmup(loader::TaskGraphsProblemLoader,solver_config,problem_dir,dummy_path = "dummy_path")
    solver = solver_config.solver
    for prob_path in readdir(problem_dir;join=true)
        is_problem_file(loader,prob_path) ? nothing : continue
        prob = load_problem(loader,solver_config,prob_path)
        solution, timer_results = profile_solver!(solver,prob)
        results_dict = compile_results(
            solver,
            solver_config.feats,
            prob,
            solution,
            timer_results
            )
        results_dict["problem_file"] = dummy_path
        write_results(loader,solver_config,prob,prob_path,results_dict)
        break
    end
    run(`rm -rf $dummy_path`)
    return loader
end


text:
@verilog function multioutput(inputwire::Wire{3:0v})
  val1 = inputwire[1:0v]
  val2 = inputwire[msb:2v]
  return (val1, val2)
end

@test multioutput(0b1100) == 0b0011

#next, write two functions.
@verilog function vf(x::Wire{3:0v})
  o1 = x[1:0v]
  o2 = x[msb:2v]
  (o1, o2)
end

@test vf() == """
module vf(
  input [3:0] x,
  output [1:0] o1,
  output [1:0] o2);

  assign o1 = x[1:0];
  assign o2 = x[3:2];
endmodule
"""

@verilog function calls_vf(z::Wire{3:0v})
  c1, c2 = vf(z)
  res = c1 & c2
end

@test calls_vf() == """
module calls_vf(
  input [3:0] z,
  output [1:0] res);

  wire [1:0] c1;
  wire [1:0] c2;

  vf vf_c1_c2(
    .x (z),
    .o1 (c1),
    .o2 (c2));

  assign res = (c1 & c2);
endmodule
"""

@verilog function calls_vf_2(z::Wire{3:0v})
  c2 = Wire(3:0v)
  c1, c2[1:0v] = vf(z)
  res = c1 & c2
end

@test calls_vf_2() == """
module calls_vf_2(
  input [3:0] z,
  output [1:0] res);

  wire [1:0] c1;
  wire [3:0] c2;

  vf vf_c1_c2_1_0(
    .x (z),
    .o1 (c1),
    .o2 (c2[1:0]));

  assign res = (c1 & c2);
endmodule
"""