text:
# ======================================================
# Structures
# ======================================================

"""
    ValGraph{V <: Integer, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C} <: AbstractValGraph

A type representing an undirected simple graph with vertex, edge and graph values.

# Parameters
- `V`: The type of the vertex indices. Should be a concrete subtype of `Integer`.
- `V_VALS`: The type of the vertex values, either a `Tuple` or `NamedTuple`.
- `E_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `G_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `V_VALS_C`: Internal storage parameter, is derived from `V_VALS`
- `E_VALS_C`: Internal storage parameter, is derived from `E_VALS`

The internal parameters `V_VALS_C` and `E_VALS_C` are automatically calculated
by the constructors so that they should usually not be manually specified.
"""
mutable struct ValGraph{  V <: Integer,
                            V_VALS <: AbstractTuple,
                            E_VALS <: AbstractTuple,
                            G_VALS <: AbstractTuple,
                            V_VALS_C <: AbstractTuple,
                            E_VALS_C <: AbstractTuple
                        } <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    vertexvals::V_VALS_C
    edgevals::E_VALS_C
    graphvals::G_VALS
end

"""
    ValOutDiGraph{V <: Integer, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C} <: AbstractValGraph

A type representing a directed simple graph with vertex, edge and graph values.

# Parameters
- `V`: The type of the vertex indices. Should be a concrete subtype of `Integer`.
- `V_VALS`: The type of the vertex values, either a `Tuple` or `NamedTuple`.
- `E_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `G_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `V_VALS_C`: Internal storage parameter, is derived from `V_VALS`
- `E_VALS_C`: Internal storage parameter, is derived from `E_VALS`

The internal parameters `V_VALS_C` and `E_VALS_C` are automatically calculated
by the constructors so that they should usually not be manually specified.

Uses less memory than `ValDiGraph` by only storing outgoing edges.
"""
mutable struct ValOutDiGraph{ V <: Integer,
                                V_VALS <: AbstractTuple,
                                E_VALS <: AbstractTuple,
                                G_VALS <: AbstractTuple,
                                V_VALS_C <: AbstractTuple,
                                E_VALS_C <: AbstractTuple
                            } <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    vertexvals::V_VALS_C
    edgevals::E_VALS_C
    graphvals::G_VALS
end

"""
    ValDiGraph{V <: Integer, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C} <: AbstractValGraph

A type representing a directed simple graph with vertex, edge and graph values.

# Parameters
- `V`: The type of the vertex indices. Should be a concrete subtype of `Integer`.
- `V_VALS`: The type of the vertex values, either a `Tuple` or `NamedTuple`.
- `E_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `G_VALS`: The type of the edge values, either a `Tuple` or `NamedTuple`.
- `V_VALS_C`: Internal storage parameter, is derived from `V_VALS`
- `E_VALS_C`: Internal storage parameter, is derived from `E_VALS`

The internal parameters `V_VALS_C` and `E_VALS_C` are automatically calculated
by the constructors so that they should usually not be manually specified.

Uses more memory than `ValOutDiGraph` by not only storing outgoing edges, but also
back edges.
"""
mutable struct ValDiGraph{    V <: Integer,
                                V_VALS <: AbstractTuple,
                                E_VALS <: AbstractTuple,
                                G_VALS <: AbstractTuple,
                                V_VALS_C <: AbstractTuple,
                                E_VALS_C <: AbstractTuple
                           } <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    badjlist::Adjlist{V}
    vertexvals::V_VALS_C
    edgevals::E_VALS_C
    redgevals::E_VALS_C
    graphvals::G_VALS
end

# ======================================================
# Constructors
# ======================================================

#  ------------------------------------------------------
#  helpers
#  ------------------------------------------------------

"""
The default eltype to use in a graph constructor when no eltype is specified.
"""
const default_eltype = Int32

@generated function create_edgevals(n, ::Type{E_VALS}) where {E_VALS <: Tuple}
    args = Expr[]
    for T in fieldtypes(E_VALS)
        push!(args, :(Adjlist{$T}(n)))
    end
    return Expr(:call, :tuple, args...)
end

@generated function create_edgevals(n, ::Type{E_VALS}) where {E_VALS <: NamedTuple}
    args = Expr[]
    for T in fieldtypes(E_VALS)
        push!(args, :(Adjlist{$T}(n)))
    end
    names = fieldnames(E_VALS)
    return Expr(:call, Expr(:curly, :NamedTuple, names), Expr(:call, :tuple, args...))
end

create_vertexvals(n, V_VALS::Type{<:NTuple{0}}, ::Nothing) = ()
create_vertexvals(n, V_VALS::Type{<:NamedNTuple{0}}, ::Nothing) = (;)

@generated function create_vertexvals(n, ::Type{V_VALS}, ::UndefInitializer) where {V_VALS <: Tuple}
    args = Expr[]
    for T in fieldtypes(V_VALS)
        push!(args, :(Vector{$T}(undef, n)))
    end
    return Expr(:call, :tuple, args...)
end

@generated function create_vertexvals(n, ::Type{V_VALS}, ::UndefInitializer) where {V_VALS <: NamedTuple}
    args = Expr[]
    for T in fieldtypes(V_VALS)
        push!(args, :(Vector{$T}(undef, n)))
    end
    names = fieldnames(V_VALS)
    return Expr(:call, Expr(:curly, :NamedTuple, names), Expr(:call, :tuple, args...))
end

function create_vertexvals(n, V_VALS::Type{<:AbstractTuple}, f::Function)
    vertexvals = create_vertexvals(n, V_VALS, undef)
    for v in 1:n
        t = convert_to_tuple(V_VALS, f(v))
        for i in 1:length(V_VALS.types)
            vertexvals[i][v] = t[i]
        end
    end
    return vertexvals
end

#  ------------------------------------------------------
#  Constructors for empty graphs
#  ------------------------------------------------------

"""
    ValGraph{V = $default_eltype}(n; vertexval_types=(), edgeval_types=(), vertexval_init=nothing, graphvals=())
    ValGraph{V, V_VALS, E_VALS}(n; vertexval_init=nothing, graphvals=())

Construct a `ValGraph` with `n` vertices, zero edges and graph values `graphvals`.  The
vertex value types are either `V_VALS` or `vertexval_types` and the edge value type are
`edgeval_types`.

If omitted, the element type `V` is $(default_eltype).
"""
function ValGraph{V, V_VALS, E_VALS}(n::Integer; vertexval_init=nothing, graphvals=()) where {V <: Integer, V_VALS <: AbstractTuple, E_VALS <: AbstractTuple}

    fadjlist = Adjlist{V}(n)
    vertexvals = create_vertexvals(n, V_VALS, vertexval_init)
    edgevals = create_edgevals(n, E_VALS)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)
    G_VALS = typeof(graphvals)

    return ValGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(0, fadjlist, vertexvals, edgevals, graphvals)
end

"""
    ValOutDiGraph{V = $default_eltype}(n; vertexval_types=(), edgeval_types=(), vertexval_init=nothing, graphvals=())
    ValOutDiGraph{V, V_VALS, E_VALS}(n; vertexval_init=nothing, graphvals=())

Construct a `ValOutDiGraph` with `n` vertices, zero edges and graph values `graphvals`.  The
vertex value types are either `V_VALS` or `vertexval_types` and the edge value type are
`edgeval_types`.

If omitted, the element type `V` is $(default_eltype).
"""
function ValOutDiGraph{V, V_VALS, E_VALS}(n::Integer; vertexval_init=nothing, graphvals=()) where {V<:Integer, V_VALS, E_VALS}

    fadjlist = Adjlist{V}(n)
    vertexvals = create_vertexvals(n, V_VALS, vertexval_init)
    edgevals = create_edgevals(n, E_VALS)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)
    G_VALS = typeof(graphvals)

    return ValOutDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(0, fadjlist, vertexvals, edgevals, graphvals)
end

"""
    ValDiGraph{V = $default_eltype}(n; vertexval_types=(), edgeval_types=(), vertexval_init=nothing, graphvals=())
    ValDiGraph{V, E_VALS}(n, vertexval_init=nothing; graphvals=())

Construct a `ValDiGraph` with `n` vertices, zero edges and graph values `graphvals`.  The
vertex value types are either `V_VALS` or `vertexval_types` and the edge value type are
`edgeval_types`.

If omitted, the element type `V` is $(default_eltype).
"""
function ValDiGraph{V, V_VALS, E_VALS}(n::Integer; vertexval_init=nothing, graphvals=()) where {V<:Integer, V_VALS, E_VALS}

    fadjlist = Adjlist{V}(n)
    badjlist = Adjlist{V}(n)
    vertexvals = create_vertexvals(n, V_VALS, vertexval_init)
    edgevals = create_edgevals(n, E_VALS)
    redgevals = create_edgevals(n, E_VALS)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)
    G_VALS = typeof(graphvals)

    return ValDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(
                0, fadjlist, badjlist, vertexvals, edgevals, redgevals, graphvals)
end

for G in (:ValGraph, :ValOutDiGraph, :ValDiGraph)

    @eval function $G(n::Integer; vertexval_types::AbstractTypeTuple=(), edgeval_types::AbstractTypeTuple=(), vertexval_init=nothing, graphvals=())

        V_VALS = typetuple_to_type(vertexval_types)
        E_VALS = typetuple_to_type(edgeval_types)
        return $G{default_eltype, V_VALS, E_VALS}(n, vertexval_init=vertexval_init, graphvals=graphvals)
    end

    @eval function $G{V}(n::Integer; vertexval_types::AbstractTypeTuple=(), edgeval_types::AbstractTypeTuple=(), vertexval_init=nothing, graphvals=()) where {V <: Integer}

        V_VALS = typetuple_to_type(vertexval_types)
        E_VALS = typetuple_to_type(edgeval_types)
        return $G{V, V_VALS, E_VALS}(V(n), vertexval_init=vertexval_init, graphvals=graphvals)
    end
end

#  ------------------------------------------------------
#  Constructors from other value graphs
#  ------------------------------------------------------

function ValGraph{V}(g::ValGraph) where {V}

    neg = ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, vertexvals, edgevals, g.graphvals)
end

ValGraph(g::ValGraph) = ValGraph{eltype(g)}(g)

function ValDiGraph{V}(g::ValDiGraph) where {V}

    neg = ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    badjlist = deepcopy_adjlist(V, g.badjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)
    redgevals= copy_edgevals(g.redgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, badjlist, vertexvals, edgevals, redgevals, g.graphvals)
end

ValDiGraph(g::ValDiGraph) = ValDiGraph{eltype(g)}(g)

function ValDiGraph{V}(g::ValGraph) where {V}

    neg = ne(g) - num_self_loops(g) + ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    badjlist = deepcopy_adjlist(V, g.fadjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)
    redgevals= copy_edgevals(g.edgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, badjlist, vertexvals, edgevals, redgevals, g.graphvals)
end

ValDiGraph(g::ValGraph) = ValDiGraph{eltype(g)}(g)

# TODO ValOutDiGraph -> ValDiGraph

function ValOutDiGraph{V}(g::ValOutDiGraph) where {V}

    neg = ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValOutDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, vertexvals, edgevals, g.graphvals)
end

ValOutDiGraph(g::ValOutDiGraph) = ValOutDiGraph{eltype(g)}(g)

function ValOutDiGraph{V}(g::ValDiGraph) where {V}

    neg = ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValOutDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, vertexvals, edgevals, g.graphvals)
end

ValOutDiGraph(g::ValDiGraph) = ValOutDiGraph{eltype(g)}(g)

function ValOutDiGraph{V}(g::ValGraph) where {V}

    neg = ne(g) - num_self_loops(g) + ne(g)
    fadjlist = deepcopy_adjlist(V, g.fadjlist)
    vertexvals = copy_vertexvals(g.vertexvals)
    edgevals = copy_edgevals(g.edgevals)

    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    G_VALS = graphvals_type(g)
    V_VALS_C = typeof(vertexvals)
    E_VALS_C = typeof(edgevals)

    return ValOutDiGraph{V, V_VALS, E_VALS, G_VALS, V_VALS_C, E_VALS_C}(neg, fadjlist, vertexvals, edgevals, g.graphvals)
end

ValOutDiGraph(g::ValGraph) = ValOutDiGraph{eltype(g)}(g)

# =========================================================
# Interface
# =========================================================

#  ------------------------------------------------------
#  nv
#  ------------------------------------------------------

Graphs.nv(g::ValGraph) = eltype(g)(length(g.fadjlist))
Graphs.nv(g::ValOutDiGraph) = eltype(g)(length(g.fadjlist))
Graphs.nv(g::ValDiGraph) = eltype(g)(length(g.fadjlist))

#  ------------------------------------------------------
#  ne
#  ------------------------------------------------------

Graphs.ne(g::ValGraph) = g.ne
Graphs.ne(g::ValDiGraph) = g.ne
Graphs.ne(g::ValOutDiGraph) = g.ne

#  ------------------------------------------------------
#  add_edge!
#  ------------------------------------------------------

function Graphs.add_edge!(g::ValGraph, s::Integer, d::Integer, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds

    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace value, but return false
        set_values_for_index!(edgevals, s, index, values)
        s == d && return false # selfloop
        index = searchsortedfirst(g.fadjlist[d], s)
        set_values_for_index!(edgevals, d, index, values)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(edgevals, s, index, values)
    g.ne += 1

    s == d && return true # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    insert_values_for_index!(edgevals, d, index, values)
    return true # edge successfully added
end

function Graphs.add_edge!(g::ValOutDiGraph, s::Integer, d::Integer, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace values, but return false
        set_values_for_index!(edgevals, s, index, values)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(edgevals, s, index, values)
    g.ne += 1
    return true # edge successfully added
end

function Graphs.add_edge!(g::ValDiGraph, s::Integer, d::Integer, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace values, but return false
        set_values_for_index!(g.edgevals, s, index, values)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(g.edgevals, s, index, values)
    g.ne += 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    insert_values_for_index!(g.redgevals, d, index, values)

    return true # edge successfully added
end

#  ------------------------------------------------------
#  rem_edge!
#  ------------------------------------------------------

function Graphs.rem_edge!(g::ValGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(edgevals, s, index)

    g.ne -= 1
    s == d && return true # self-loop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    delete_values_for_index!(edgevals, d, index)

    return true
end

function Graphs.rem_edge!(g::ValOutDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(edgevals, s, index)

    g.ne -= 1
    return true
end

function Graphs.rem_edge!(g::ValDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(g.edgevals, s, index)

    g.ne -= 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    delete_values_for_index!(g.redgevals, d, index)

    return true
end

#  ------------------------------------------------------
#  add_vertex!
#  ------------------------------------------------------

function Graphs.add_vertex!(g::ValGraph, values)

    V = eltype(g)
    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    values = convert_to_tuple(V_VALS, values)

    if V != BigInt && nv(g) == typemax(V)
        # maybe consider capping at nv(g) - 1
        return false # cannot add more vertices
    end

    push!(g.fadjlist, V[])
    for i in 1:length(g.edgevals)
        push!(g.edgevals[i], E_VALS.types[i][])
    end

    for i in 1:length(g.vertexvals)
        push!(g.vertexvals[i], values[i])
    end

    return true
end

function Graphs.add_vertex!(g::ValOutDiGraph, values)

    V = eltype(g)
    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    values = convert_to_tuple(V_VALS, values)

    if V != BigInt && nv(g) == typemax(V)
        # maybe consider capping at nv(g) - 1
        return false # cannot add more vertices
    end

    push!(g.fadjlist, V[])
    for i in 1:length(g.edgevals)
        push!(g.edgevals[i], E_VALS.types[i][])
    end

    for i in 1:length(g.vertexvals)
        push!(g.vertexvals[i], values[i])
    end

    return true
end

function Graphs.add_vertex!(g::ValDiGraph, values)

    V = eltype(g)
    V_VALS = vertexvals_type(g)
    E_VALS = edgevals_type(g)
    values = convert_to_tuple(V_VALS, values)

    if V != BigInt && nv(g) == typemax(V)
        # maybe consider capping at nv(g) - 1
        return false # cannot add more vertices
    end

    push!(g.fadjlist, V[])
    push!(g.badjlist, V[])
    for i in 1:length(g.edgevals)
        push!(g.edgevals[i], E_VALS.types[i][])
        push!(g.redgevals[i], E_VALS.types[i][])
    end

    for i in 1:length(g.vertexvals)
        push!(g.vertexvals[i], values[i])
    end

    return true
end

#  ------------------------------------------------------
#  rem_vertex!
#  ------------------------------------------------------

# TODO

#  ------------------------------------------------------
#  has_edge
#  ------------------------------------------------------

function Graphs.has_edge(g::ValGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        d = s
        list_s = list_d
    end
    return Graphs.insorted(d, list_s)
end

function Graphs.has_edge(g::ValOutDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]

    return Graphs.insorted(d, list_s)
end

function Graphs.has_edge(g::ValDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_fadj = g.fadjlist[s]
    @inbounds list_badj = g.badjlist[d]
    if length(list_fadj) > length(list_badj)
        d = s
        list_fadj = list_badj
    end
    return Graphs.insorted(d, list_fadj)
end

#  ------------------------------------------------------
#  get_edgeval
#  ------------------------------------------------------

@inbounds function get_edgeval(g::ValGraph, s::Integer, d::Integer, key::Integer)

    hasedgekey_or_throw(g, key) # TODO might be sufficient to just check index
    return _get_edgeval(g, s, d, g.edgevals[key])
end

@inbounds function get_edgeval(g::ValGraph, s::Integer, d::Integer, key::Symbol)

    hasedgekey_or_throw(g, key)
    return _get_edgeval(g, s, d, g.edgevals[key])
end

function _get_edgeval(g::ValGraph, s::Integer, d::Integer, adjlist::Adjlist)

    verts = vertices(g)
    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end

    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    # TODO more specific error
    error("No such edge")
end

@inline function get_edgeval(g::ValOutDiGraph, s::Integer, d::Integer, key::Integer)

    hasedgekey_or_throw(g, key)
    return _get_edgeval(g, s, d, g.edgevals[key])
end

@inbounds function get_edgeval(g::ValOutDiGraph, s::Integer, d::Integer, key::Symbol)

    hasedgekey_or_throw(g, key)
    return _get_edgeval(g, s, d, g.edgevals[key])
end

function _get_edgeval(g::ValOutDiGraph, s::Integer, d::Integer, adjlist::Adjlist)

    verts = vertices(g)
    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    error("No such edge")
end

@inline function get_edgeval(g::ValDiGraph, s::Integer, d::Integer, key::Integer)

    hasedgekey_or_throw(g, key)

    return _get_edgeval(g, s, d, g.edgevals[key])
end

@inbounds function get_edgeval(g::ValDiGraph, s::Integer, d::Integer, key::Symbol)

    hasedgekey_or_throw(g, key)
    return _get_edgeval(g, s, d, g.edgevals[key])
end

function _get_edgeval(g::ValDiGraph, s::Integer, d::Integer, adjlist::Adjlist)

    verts = vertices(g)
    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    error("No such edge")
end

#  ------------------------------------------------------
#  get_edgeval_or
#  ------------------------------------------------------

@inline function get_edgeval_or(g::ValGraph, s::Integer, d::Integer, key::Integer, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

@inline function get_edgeval_or(g::ValGraph, s::Integer, d::Integer, key::Symbol, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

function _get_edgeval_or(g::ValGraph, s::Integer, d::Integer, adjlist::Adjlist, alternative)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end

    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    return alternative
end

@inline function get_edgeval_or(g::ValOutDiGraph, s::Integer, d::Integer, key::Integer, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

@inline function get_edgeval_or(g::ValOutDiGraph, s::Integer, d::Integer, key::Symbol, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

function _get_edgeval_or(g::ValOutDiGraph, s::Integer, d::Integer, adjlist::Adjlist, alternative)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    return alternative
end

@inline function get_edgeval_or(g::ValDiGraph, s::Integer, d::Integer, key::Integer, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

@inline function get_edgeval_or(g::ValDiGraph, s::Integer, d::Integer, key::Symbol, alternative)

    hasedgekey_or_throw(g, key)
    return _get_edgeval_or(g, s, d, g.edgevals[key], alternative)
end

function _get_edgeval_or(g::ValDiGraph, s::Integer, d::Integer, adjlist::Adjlist, alternative)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return adjlist[s][index]
    end
    return alternative
end

#  ------------------------------------------------------
#  get_edgeval(g, s, d, :)
#  ------------------------------------------------------

function get_edgeval(g::ValGraph, s::Integer, d::Integer, ::Colon)

    verts = vertices(g)
    (s in verts && d in verts) || error("Values not found")
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    error("Values not found")
end

function get_edgeval(g::ValOutDiGraph, s::Integer, d::Integer, ::Colon)

    verts = vertices(g)
    (s in verts && d in verts) ||  error("Values not found")

    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    error("Values not found")
end

# TODO could probably be made faster by checking the shorter list
function get_edgeval(g::ValDiGraph, s::Integer, d::Integer, ::Colon)

    verts = vertices(g)
    (s in verts && d in verts) || error("Values not found")

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    error("Values not found")
end

#  -----------------------------------------------------
#  get_edgeval_or(g, s, d, :)
#  -----------------------------------------------------

function get_edgeval_or(g::ValGraph, s::Integer, d::Integer, ::Colon, alternative)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    return alternative
end

function get_edgeval_or(g::ValOutDiGraph, s::Integer, d::Integer, ::Colon, alternative)

    verts = vertices(g)
    (s in verts && d in verts) ||  return alternative

    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    return alternative
end

# TODO could probably be made faster by checking the shorter list
function get_edgeval_or(g::ValDiGraph, s::Integer, d::Integer, ::Colon, alternative)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return values_for_index(g.edgevals, edgevals_type(g), s, index)
    end
    return alternative
end

#  ------------------------------------------------------
#  set_edgeval!
#  ------------------------------------------------------

function set_edgeval!(g::ValGraph, s::Integer, d::Integer, key::Integer, value)

    hasedgekey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    @inbounds edgevals[d][index] = value

    return true
end

function set_edgeval!(g::ValOutDiGraph, s::Integer, d::Integer, key::Integer, value)

    hasedgekey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    return true

end

function set_edgeval!(g::ValDiGraph, s::Integer, d::Integer, key::Integer, value)

    hasedgekey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    redgevals = g.redgevals[key]

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    @inbounds redgevals[d][index] = value

    return true
end

function set_edgeval!(g::ValGraph, s::Integer, d::Integer, ::Colon, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)

    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(edgevals, s, index, values)

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    set_values_for_index!(edgevals, d, index, values)

    return true
end

function set_edgeval!(g::ValOutDiGraph, s::Integer, d::Integer, ::Colon, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)

    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(edgevals, s, index, values)

    return true
end

function set_edgeval!(g::ValDiGraph, s::Integer, d::Integer, ::Colon, values)

    E_VALS = edgevals_type(g)
    values = convert_to_tuple(E_VALS, values)

    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(g.edgevals, s, index, values)

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    set_values_for_index!(g.redgevals, d, index, values)

    return true
end

#  ------------------------------------------------------
#  get_vertexval
#  ------------------------------------------------------

# TODO maybe implement also for Colon

@inline function get_vertexval(g::ValGraph, v::Integer, key::Integer)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

@inline function get_vertexval(g::ValGraph, v::Integer, key::Symbol)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

_get_vertexval(g::ValGraph, v::Integer, list::Vector) = list[v]

@inline function get_vertexval(g::ValOutDiGraph, v::Integer, key::Integer)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

@inline function get_vertexval(g::ValOutDiGraph, v::Integer, key::Symbol)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

_get_vertexval(g::ValOutDiGraph, v::Integer, list::Vector) = list[v]

@inline function get_vertexval(g::ValDiGraph, v::Integer, key::Integer)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

@inline function get_vertexval(g::ValDiGraph, v::Integer, key::Symbol)

    # TODO verify key, vertex
    return _get_vertexval(g, v, g.vertexvals[key])
end

_get_vertexval(g::ValDiGraph, v::Integer, list::Vector) = list[v]

#  ------------------------------------------------------
#  set_vertexval!
#  ------------------------------------------------------

function set_vertexval!(g::ValGraph, v::Integer, key::Integer, value)

    has_vertex(g, v) || return false
    g.vertexvals[key][v] = value

    return true
end

function set_vertexval!(g::ValOutDiGraph, v::Integer, key::Integer, value)

    has_vertex(g, v) || return false
    g.vertexvals[key][v] = value

    return true
end

function set_vertexval!(g::ValDiGraph, v::Integer, key::Integer, value)

    has_vertex(g, v) || return false
    g.vertexvals[key][v] = value

    return true
end

#  ------------------------------------------------------
#  get_graphval
#  ------------------------------------------------------

function get_graphval(g::ValGraph, key::Integer)

    # TODO verify key

    return g.graphvals[key]
end

function get_graphval(g::ValOutDiGraph, key::Integer)

    # TODO verify key

    return g.graphvals[key]
end

function get_graphval(g::ValDiGraph, key::Integer)

    # TODO verify key

    return g.graphvals[key]
end

function get_graphval(g::ValGraph, ::Colon)

    return g.graphvals
end

function get_graphval(g::ValOutDiGraph, ::Colon)

    return g.graphvals
end

function get_graphval(g::ValDiGraph, ::Colon)

    return g.graphvals
end

#  ------------------------------------------------------
#  set_graphval!
#  ------------------------------------------------------

function set_graphval!(g::ValGraph, key::Integer, value)

    g.graphvals = replace_in_tuple(g.graphvals, key, value)
    return true
end

function set_graphval!(g::ValOutDiGraph, key::Integer, value)

    g.graphvals = replace_in_tuple(g.graphvals, key, value)
    return true
end

function set_graphval!(g::ValDiGraph, key::Integer, value)

    g.graphvals = replace_in_tuple(g.graphvals, key, value)
    return true
end

#  ------------------------------------------------------
#  is_directed
#  ------------------------------------------------------

Graphs.is_directed(::Type{<:ValGraph}) = false
Graphs.is_directed(::Type{<:ValOutDiGraph}) = true
Graphs.is_directed(::Type{<:ValDiGraph}) = true

#  ------------------------------------------------------
#  outneighbors
#  ------------------------------------------------------

Graphs.outneighbors(g::ValGraph, v::Integer) = g.fadjlist[v]
Graphs.outneighbors(g::ValDiGraph, v::Integer) = g.fadjlist[v]
Graphs.outneighbors(g::ValOutDiGraph, v::Integer) = g.fadjlist[v]

#  ------------------------------------------------------
#  inneighbors
#  ------------------------------------------------------

Graphs.inneighbors(g::ValGraph, v::Integer) = outneighbors(g, v)
Graphs.inneighbors(g::ValDiGraph, v::Integer) = g.badjlist[v]

#  ------------------------------------------------------
#  outedgevals
#  ------------------------------------------------------

outedgevals(g::ValGraph, v::Integer, key::Integer) = g.edgevals[key][v]
outedgevals(g::ValDiGraph, v::Integer, key::Integer) = g.edgevals[key][v]
outedgevals(g::ValOutDiGraph, v::Integer, key::Integer) = g.edgevals[key][v]

#  ------------------------------------------------------
#  inedgevals
#  ------------------------------------------------------

inedgevals(g::ValGraph, v::Integer, key::Integer) = outedgevals(g, v, key)
inedgevals(g::ValDiGraph, v::Integer, key::Integer) = g.redgevals[key][v]

# ====================================================================
# Iterators
# ====================================================================

@inline function Base.iterate(eit::ValEdgeIter{<:ValGraph, key}, state=(one(eltype(eit.graph)), 1)) where {key}

    g = eit.graph
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    V = eltype(g)
    n = V(nv(g))
    u, i = state

    @inbounds while u < n
        list_u = fadjlist[u]
        if i > length(list_u)
            u += one(u)
            i = searchsortedfirst(fadjlist[u], u)
            continue
        end
        e = if key == nothing
                ValEdge(u, list_u[i], ())
            else
                ValEdge(u, list_u[i], values_for_index(edgevals, edgevals_type(g), u, i))
            end

        state = (u, i + 1)
        return e, state
    end

    # i > length(fadjlist[end]) || fadlist[end][i] == n

    @inbounds (n == 0 || i > length(fadjlist[n])) && return nothing

    e = ValEdge(n, n, values_for_index(edgevals, edgevals_type(g), u, i))
    state = (u, i + 1)
    return e, state
end

function Base.iterate(
            iter::ValEdgeIter{<:Union{ValOutDiGraph, ValDiGraph}, key},
            state=(one(eltype(iter.graph)), 1)) where {key}

    g = iter.graph
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    V = eltype(g)
    n = V(nv(g))
    u, i = state

    @inbounds while u <= n
        if i > length(fadjlist[u])
            u == n && return nothing

            u += one(u)
            i = 1
            continue
        end
        e = if key == nothing
                ValDiEdge(u, fadjlist[u][i], ())
            else
                ValDiEdge(u, fadjlist[u][i], values_for_index(edgevals, edgevals_type(g), u, i))
            end
        return e, (u, i + 1)
    end

    return nothing
end


text:
# module TestOperations

using Test, Dates, Temporal, Random

@testset "Operations" begin
    @testset "Utilities" begin
        Random.seed!(1)
        x = ts(cumsum(randn(N))) + N
        dx = diff(log(x), pad=true, padval=NaN)
        @test findall(trues(x)) == [CartesianIndex(i,1) for i in 1:size(x,1)]
        @test findfirst(trues(x)) == CartesianIndex(1,1)
        idx_nans = isnan(dx)
        @test findall(idx_nans) == [CartesianIndex(1,1)]
        @test findfirst(idx_nans) == CartesianIndex(1,1)
        @test sign(dx[findall(!idx_nans)]).values[:] == sign.(dx.values[.!isnan.(dx.values)])
        x = ts(rand(N))
        @test round(Int,x).values == round.(Int,round(x).values)
        @test ones(x).values == ones(eltype(x), size(x))
        @test zeros(x).values == zeros(eltype(x), size(x))
        @test length(rand(TS,N)) == N
        @test length(rand(TS,N,K)) == N*K
        @test length(rand(TS,(N,K))) == N*K
        @test length(randn(TS,N)) == N
        @test length(randn(TS,N,K)) == N*K
        @test length(randn(TS,(N,K))) == N*K
        x = ts(-1.0 .* rand(N))
        y = +x
        @test all(y.values .<= 0.0)
        y = -x
        @test all(y.values .>= 0.0)
    end
    @testset "Logical" begin
        x1 = x2 = TS(rand(100))
        @test x1 == x2
        @test all(trues(x1))
        @test all(!falses(x2))
        x1 += 1.0
        @test x1 != x2
        @test all(x1 > x2)
        @test all(x2 < x1)
        @test all(x1 >= x2)
        @test all(x2 <= x1)
    end
    @testset "Arithmetic" begin
        x1 = TS(rand(100)) + 100
        x2 = TS(rand(100)) + 100
        z = x1 + x2
        @test z.values == x1.values + x2.values
        z = x1 - x2
        @test z.values == x1.values - x2.values
        z = x1 * x2
        @test z.values == x1.values .* x2.values
        z = x1 / x2
        @test z.values == x1.values ./ x2.values
        z = x1 * 2.0
        @test z.values == x1.values * 2.0
        z = x1 / 2.0
        @test z.values == x1.values / 2.0
        z = x1 % 2.0
        @test z.values == x1.values .% 2.0
        @test x1 + x1.values == x1 * 2
        @test x2 - x2.values == zeros(x2)
    end
    @testset "Statistics" begin
        x = TS(rand(100)) + 100
        x[1] = 0.0
        @test prod(x) == prod(round(x)) == cumprod(x).values[end] == prod(x.values) == 0.0
        @test cumsum(x).values[end] ≈ sum(x) ≈ sum(x.values)
        @test mean(x) == mean(x.values)
        dx = diff(x; pad=true, padval=NaN)
        @test isnan(dx.values[1])
        @test size(dx,1) == size(x,1)
        @test size(dropnan(dx),1) == size(diff(x,pad=false),1)
        dx = pct_change(x)
        @test size(dx,1) == size(x,1)-1
        dx = pct_change(x, continuous=false, pad=false)
        @test size(dx,1) == size(x,1)-1
        dx = pct_change(x, continuous=false, pad=true)
        @test size(dx,1) == size(x,1)
        bx = lag(x, pad=true, padval=NaN)
        @test isnan(bx.values[1])
        @test size(bx,1) == size(x,1)
        @test size(dropnan(bx),1) == size(lag(x,pad=false),1)
        @test dropnan(bx) == dropnan(shift(x, pad=true, padval=NaN))
        @test acf(x)[1] == 1.0
    end
end

# end

ext: jl


text:
@testset "Q_update_test" begin

#sparse matrix
  A=sprandn(100,100,0.1)
  B=Vector{SparseMatrixCSC{Float64,Int64}}(undef, 2)
  B[1]=sprandn(100,100,0.1)
  B[2]=sprandn(100,100,0.1)
  #make sure all nonzero diagonals in B are also in A, otherwise CDS doesnt work
  A=A+B[1]+B[2]
  A2=deepcopy(A)
  B2=deepcopy(B)
  A3=deepcopy(A)
  B3=deepcopy(B)
  maxit=100;
  p=2
  log_PARSDMM = log_type_PARSDMM(zeros(maxit,p-1),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());
  rho=Vector{Float64}(undef, 2)
  rho[1]=1.0
  rho[2]=1.0
  log_PARSDMM.rho[1,1]=2.0
  log_PARSDMM.rho[1,2]=3.0

  #add A = A+ 1.0*B[1] + 2.0*B[2]
  i=1
  ind_updated = findall(rho .!= log_PARSDMM.rho[i,:])  #::Vector{Integer}# locate changed rho index
  ind_updated = convert(Array{Int64,1},ind_updated)
  #use explicit solution:
    if isempty(ind_updated) == false
    for ii=1:length(ind_updated)
        A = A + ( B[ind_updated[ii]] )*( rho[ind_updated[ii]]-log_PARSDMM.rho[i,ind_updated[ii]] );
    end
    end
  #use Q_update:
    TI=Int64
    TF=Float64
    TD_Prop=set_properties(zeros(Bool, 99), zeros(Bool, 99), zeros(Bool, 99),
                           fill!(Vector{Tuple{TI,TI}}(undef, 99), (0, 0)),
                           fill!(Vector{Tuple{String,String, String, String}}(undef, 99), ("", "", "", "")),
                           zeros(Bool, 99),
                           Vector{Vector{TI}}(undef,99))
    A2=Q_update!(A2,B2,TD_Prop,rho,ind_updated,log_PARSDMM,i,[])

    x=randn(size(A,2))
    @test isapprox(A2*x,A*x,rtol=10*eps())

#sparse matrix in CDS format
  (R_A,offset_A)=mat2CDS(A3)
  (R_B1,offset_B1)=mat2CDS(B[1])
  (R_B2,offset_B2)=mat2CDS(B[2])
  R_B=Vector{Array{Float64,2}}(undef, 2)
  R_B[1]=R_B1
  R_B[2]=R_B2
  TD_Prop.AtA_offsets[1]=offset_B1
  TD_Prop.AtA_offsets[2]=offset_B2
  Q_update!(R_A,R_B,TD_Prop,rho,ind_updated,log_PARSDMM,i,offset_A)

  y=zeros(length(x))
  y=CDS_MVp(length(x),size(R_A,2),R_A,offset_A,x,y)
  @test isapprox(y,A*x,rtol=10*eps())

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# Algorithms extending

"""
# Usage

    solve(prob::MultiTermsFODEProblem, ClosedForm())

Use Closed-Form solution to obtain numerical solution at zero initial condition.

## Reference:

Dingyu Xue, Northeastern University, China
ISBN:9787030543981
"""
struct ClosedForm <: FractionalDiffEqAlgorithm end

function solve(prob::MultiTermsFODEProblem, h, ::ClosedForm)
    @unpack parameters, orders, rightfun, rparameters, rorders, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    D = sum(parameters./(h.^orders))
    t = collect(t0:h:T)
    nT = length(t)
    u = rightfun.(t)

    #TODO: For small orders? e.g. if rparameters and rorders is Number???
    isa(rparameters, Number) && typeof(rorders) <: Number ? D1 = rparameters/h^rorders : D1 = rparameters[:]./h.^rorders[:]

    nA = length(parameters)
    vect = hcat(orders, rorders)

    y1 = zeros(nT)
    W = ones(nT, length(vect))

    @fastmath @inbounds @simd for j=2:nT
        W[j, :] = W[j-1, :].*(1 .-(vect'.+1)/(j.-1))
    end

    @fastmath @inbounds @simd for i=2:nT
        A = y1[i-1:-1:1]'*W[2:i, 1:nA]
        y1[i] = (u[i]-sum(A.*parameters./(h.^orders)))/D
    end

    y = zeros(nT) # Prelocate the final result

    @fastmath @inbounds @simd for i=2:nT
        y[i] = (W[1:i, (nA+1):end]*D1)'*y1[i:-1:1]
    end
    return FODESolution(t, y)
end


text:
"""
`Base.length(tab::ODERKTableau)`

Defines the length of a Runge-Kutta method to be the number of stages.
"""
Base.length(tab::ODERKTableau) = tab.stages

"""
`stability_region(z,tab::ODERKTableau)`

Calculates the stability function from the tableau at `z`. Stable if <1.

```math
r(z) = \\frac{\\det(I-zA+zeb^T)}{\\det(I-zA)}
```
"""
stability_region(z,tab::ODERKTableau) = det(eye(tab.stages)- z*tab.A + z*ones(tab.stages)*tab.α')/det(eye(tab.stages)-z*tab.A)

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
include(joinpath(nlpmodels_path, "nls_consistency.jl"))

function consistent_nls()
  @testset "Consistency of Linear problem" begin
    m, n = 50, 20
    A = Matrix(1.0I, m, n) .+ 1
    b = collect(1:m)
    lvar = -ones(n)
    uvar = ones(n)
    lls_model = LLSModel(A, b, lvar=lvar, uvar=uvar)
    autodiff_model = ADNLSModel(x->A*x-b, zeros(n), m, lvar=lvar, uvar=uvar)
    nlp = ADNLPModel(x->0, zeros(n), lvar=lvar, uvar=uvar, c=x->A*x-b,
                     lcon=zeros(m), ucon=zeros(m))
    feas_res_model = FeasibilityResidual(nlp)
    model = Model()
    @variable(model, x[1:n], start=0.0)
    @NLexpression(model, F[i=1:m], sum(A[i,j] * x[j] for j = 1:n) - b[i])
    mpnls = MathProgNLSModel(model, F)
    nlss = [lls_model, autodiff_model, feas_res_model, mpnls]
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_nls_functions(nlss)
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    for nls in nlss
      reset!(nls)
    end

    f(x) = begin
      r = A*x - b
      return 0.5*dot(r, r)
    end
    nlps = [nlss; ADNLPModel(f, zeros(n))]
    consistent_functions(nlps, exclude=[hess_coord])
  end

  @testset "Consistency of Linear problem with linear constraints" begin
    m, n = 50, 20
    A = Matrix(1.0I, m, n) .+ 1
    b = collect(1:m)
    lvar = -ones(n)
    uvar = ones(n)
    nc = 10
    C = [ones(nc, n); 2 * ones(nc, n); -ones(nc, n); -Matrix(1.0I, nc, n)]
    lcon = [   zeros(nc); -ones(nc); fill(-Inf,nc); zeros(nc)]
    ucon = [fill(Inf,nc);  ones(nc);     zeros(nc); zeros(nc)]
    K = ((1:4:4nc) .+ (0:3)')[:]
    lcon, ucon = lcon[K], ucon[K]
    lls_model = LLSModel(A, b, lvar=lvar, uvar=uvar, C=C, lcon=lcon,
                         ucon=ucon)
    autodiff_model = ADNLSModel(x->A*x-b, zeros(n), m, lvar=lvar,
                                uvar=uvar, c=x->C*x, lcon=lcon,
                                ucon=ucon)
    model = Model()
    @variable(model, x[1:n], start=0.0)
    @NLexpression(model, F[i=1:m], sum(A[i,j] * x[j] for j = 1:n) - b[i])
    @NLconstraint(model, [i=1:4nc], lcon[i] ≤ sum(C[i,j] * x[j] for j = 1:n) ≤ ucon[i])
    mpnls = MathProgNLSModel(model, F)
    nlss = [lls_model, autodiff_model, mpnls]
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_nls_functions(nlss)
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_functions(nlss, exclude=[hess_coord])
  end

  @testset "Consistency of Nonlinear problem" begin
    m, n = 10, 2
    lvar = -ones(n)
    uvar =  ones(n)
    F(x) = [2 + 2i - exp(i*x[1]) - exp(i*x[2]) for i = 1:m]
    x0 = [0.3; 0.4]
    autodiff_model = ADNLSModel(F, x0, m, lvar=lvar, uvar=uvar)
    nlp = ADNLPModel(x->0, x0, lvar=lvar, uvar=uvar, c=F, lcon=zeros(m), ucon=zeros(m))
    feas_res_model = FeasibilityResidual(nlp)
    model = Model()
    @variable(model, x[i=1:n], start=x0[i])
    @NLexpression(model, F1[i=1:m], 2 + 2i - exp(i*x[1]) - exp(i*x[2]))
    mpnls = MathProgNLSModel(model, F1)
    nlss = [autodiff_model, feas_res_model, mpnls]
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_nls_functions(nlss)
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    for nls in nlss
      reset!(nls)
    end

    f(x) = begin
      r = F(x)
      return 0.5*dot(r, r)
    end
    nlps = [nlss; ADNLPModel(f, zeros(n))]
    consistent_functions(nlps, exclude=[hess_coord])
  end

  @testset "Consistency of Nonlinear problem with constraints" begin
    m, n = 10, 2
    lvar = -ones(n)
    uvar =  ones(n)
    F(x) = [2 + 2i - exp(i*x[1]) - exp(i*x[2]) for i = 1:m]
    x0 = [0.3; 0.4]
    c(x) = [x[1]^2 - x[2]^2; 2 * x[1] * x[2]; x[1] + x[2]]
    lcon = [0.0; -1.0; -Inf]
    ucon = [Inf;  1.0;  0.0]

    autodiff_model = ADNLSModel(F, x0, m, lvar=lvar, uvar=uvar,
                                lcon=lcon, ucon=ucon, c=c)
    model = Model()
    @variable(model, x[i=1:n], start=x0[i])
    @NLconstraint(model, x[1]^2 - x[2]^2 ≥ 0.0)
    @NLconstraint(model, -1.0 ≤ 2 * x[1] * x[2] ≤ 1.0)
    @NLconstraint(model, x[1] + x[2] ≤ 0.0)
    @NLexpression(model, F1[i=1:m], 2 + 2i - exp(i*x[1]) - exp(i*x[2]))
    mpnls = MathProgNLSModel(model, F1)

    nlss = [autodiff_model, mpnls]
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_nls_functions(nlss)
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_functions(nlss, exclude=[hess_coord])
  end

  @testset "Consistency of slack variant" begin
    model = Model()
    x0 = [0.3; 0.4; 0.5]
    @variable(model, x[i=1:3], start=x0[i])
    @NLexpression(model, F1, x[1] - 1.0)
    @NLexpression(model, F2, x[2] - x[1]^2)
    @NLexpression(model, F3, sin(x[1] * x[2]) * x[3])
    @NLconstraint(model, x[1]^2 - x[1]^2 ≥ 0.0)
    @NLconstraint(model, x[1] + x[2] ≤ 0.0)
    @NLconstraint(model, cos(x[1]) - x[2] == 0.0)
    nls = MathProgNLSModel(model, [F1; F2; F3])
    nls_auto_slack = SlackNLSModel(nls)

    smodel = Model()
    x0 = [x0; zeros(2)]
    ℓ = [-Inf; -Inf; -Inf; 0.0; -Inf]
    u = [ Inf;  Inf;  Inf; Inf;  0.0]
    @variable(smodel, ℓ[i] ≤ x[i=1:5] ≤ u[i], start=x0[i])
    @NLexpression(smodel, F1s, x[1] - 1.0)
    @NLexpression(smodel, F2s, x[2] - x[1]^2)
    @NLexpression(smodel, F3s, sin(x[1] * x[2]) * x[3])
    @NLconstraint(smodel, x[1]^2 - x[1]^2 - x[4] == 0.0)
    @NLconstraint(smodel, x[1] + x[2] - x[5] == 0.0)
    @NLconstraint(smodel, cos(x[1]) - x[2] == 0.0)
    nls_manual_slack = MathProgNLSModel(smodel, [F1s; F2s; F3s])

    nlss = [nls_manual_slack, nls_auto_slack]
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_nls_functions(nlss)
    consistent_nls_counters(nlss)
    consistent_counters(nlss)
    consistent_functions(nlss)
  end
end

consistent_nls()


text:
using RickTracy
using Base.Test

sub2ind_nonsquare(row_lens, k) = begin
    sum = 0
    if length(k) > 1
        for i in 1:(k[1]-1)
            sum += row_lens[i]
        end
    end
    sum += k[end]
end

dsum = 0
Nxs = rand(1:7)
Nys = [rand(1:10) for nx in 1:Nxs]
Nzs = [rand(1:10) for nx in 1:Nxs for ny in 1:Nys[nx]]
# Nxs = 3
# Nys = [4,3,4]
# Nzs = [5,4,5,3,7,4,4,10,8,9,3]
@show Nxs Nys Nzs

@testset "iter extrema" begin

    @resetallsnaps
    @watch x y z dsum
    i = 0
    for x in 1:Nxs
        for y in 1:Nys[x]
            i += 1
            for z in 1:Nzs[i]
                @snapall loc=inner
                dsum += 1
            end
        end
    end

    ies = RickTracy.get_iter_extrema([:x, :y, :z])
    maxNs = [[Nxs], Nys, Nzs]
    for i in 1:3
        for (k,v) in ies[i]
            j = i==1 ? 1 : sub2ind_nonsquare(maxNs[i-1], k)
            maxv = maxNs[i][j]
            @test (k,v) == (k,(1,maxv))
        end
    end
end

max_stars_repo_name: JobJob/RickTracy


text:
module powertrain_control
using ReachabilityAnalysis, ModelingToolkit
#using Polyhedra, CDDLib # ??

vars = @variables p, lam, pe, i, t

@taylorize function startup!(dx, x, params, t)
    p, lam, pe, i = x
    dx[1] = 0.41328*(2*247*(-2.3421*p*p+2.7799*p-0.3273)-0.9*(-3.66+0.08979*104.71975511*p-0.0337*104.71975511*p*p+0.0001*104.71975511*104.71975511*p))
    dx[2] = 4*(13.893-35.2518*1*((1/14.7)*(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))+20.7364*1*1*((1/14.7)*(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))*((1/14.7)*(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))+2.6287*(0.9*(-3.66+0.08979*104.71975511*p-0.0337*104.71975511*p*p+0.0001*p*104.71975511*104.71975511))-1.592*(0.9*(-3.66+0.08979*104.71975511*p-0.0337*104.71975511*p*p+0.0001*p*104.71975511*104.71975511))*1*((1/14.7)*(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))-lam)
    dx[3] = 0.41328*(2*1*(247)*(-2.3421*p*p+2.7799*p-0.3273)-(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))
    dx[4] = zero(i)
    return dx
end

@taylorize function normal!(dx, x, params, t)
    p, lam, pe, i = x
    dx[1] = 0.41328*(2*1*(247)*(-2.3421*p*p+2.7799*p-0.3273)-(-3.66+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))
    dx[2] = 4*(13.893-35.2518*1*((1/14.7)*(1+ivalue+0.04*(1*lam-14.7))*(-0.366+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))+20.7364*1*1*((1/14.7)*(1+ivalue+0.04*(1*lam-14.7))*(-0.366+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))*((1/14.7)*(1+ivalue+0.04*(1*lam-14.7))*(-0.366+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))+2.6287*(0.9*(-0.366+0.08979*104.71975511*p-0.0337*104.71975511*p*p+0.0001*p*104.71975511*104.71975511))-1.592*(0.9*(-0.366+0.08979*104.71975511*p-0.0337*104.71975511*p*p+0.0001*p*104.71975511*104.71975511)) *1*((1/14.7)*(1+ivalue+0.04*(1*lam-14.7))*(-0.366+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))-lam)
    dx[3] = 0.41328*(2*1*(247.0)*(-2.3421*p*p+2.7799*p-0.3273)-(-0.366+0.08979*104.71975511*pe-0.0337*104.71975511*pe*pe+0.0001*pe*104.71975511*104.71975511))
    dx[4] = 0.14*(1*lam-14.7)
    return dx
end

function powertrain_control_hybrid()
    n = 4 + 1 # variables

    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)

    ## mode 1 "startup"
    invariant = HPolyhedron([HalfSpace(t >= 9.5, vars), HalfSpace(t <= 20, vars)])
    mode1 = @system(x' = startup!(x), dim:5, x ∈ invariant)
    ## mode 1 "normal"
    invariant = Universe(n)
    mode2 = @system(x' = normal!(x), dim:5, x ∈ invariant)
    modes = [mode1, mode2]

    reset = Dict(n => 0.)

    ## transition startup -> normal
    guard = HalfSpace(t >= 9.5, vars)
    trans1 = ConstrainedResetMap(n, guard, reset)
    resetmaps = [trans1]

    return HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])
end

function model(X0)
    H = powertrain_control_hybrid()
    return IVP(H, [(1, X0)])
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using HorizonSideRobots

function reverse_side(side::HorizonSide)
    return HorizonSide((Int(side) + 2) % 4)
end

function draw_cross(robot::Robot)
    sides = [Nord, West, Sud, Ost]
    for side in sides
        while !isborder(robot, side)
            move!(robot, side)
            putmarker!(robot)
        end

        rside = reverse_side(side)

        while ismarker(robot)
            move!(robot, rside)
        end
    end
    putmarker!(robot)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
module MoistAir

    using Unitful
    using ThermoState
    using Roots

    using StaticArrays
    using ThermoState.Types
    using ThermoState.QuickStates

    import ThermoState: pressure,temperature,mass,moles,molar_mass
    import ThermoState: mass_volume, mol_volume, total_volume
    import ThermoState: mass_density, mol_density
    import ThermoState: mass_enthalpy, mol_enthalpy,total_enthalpy
    import ThermoState: mass_gibbs, mol_gibbs, total_gibbs
    import ThermoState: mass_helmholtz, mol_helmholtz, total_helmholtz
    import ThermoState: mass_internal_energy, mol_internal_energy, total_internal_energy
    import ThermoState: mass_entropy, mol_entropy, total_entropy
    import ThermoState: mass_fraction, mol_fraction
    import ThermoState: mass_number, mol_number

    include("types.jl")
#=
,WebBulbPT  => :hum_wetbulb
,HumidityRatioPT => :hum_ratio
,MolHumidityPT => :hum_molfrac
,MassHumidityPT => :hum_massfrac
,RelativeHumidityPT => :rel_hum
,DewPointPT => :hum_dewpoint

=#
    export ASHRAEMoistAir, ASHRAEDryAir, ASHRAEWaterVapor
    export MoistAirModel

    #dry quantities
    export dry_volume,dry_enthalpy,dry_entropy,molecular_weight

    #specific humidity properties
    export hum_wetbulb,hum_dewpoint
    export hum_molfrac,hum_massfrac
    export rel_hum,hum_ratio

end


text:
function TV.as(d::ProductMeasure{F,A}) where {F,A<:AbstractArray}
    d1 = marginals(d).f(first(marginals(d).data))
    as(Array, as(d1), size(marginals(d))...)
end

###############################################################################
# I <: Base.Generator

function TV.as(d::ProductMeasure{F,I}) where {F, I<:Base.Generator}
    d1 = marginals(d).f(first(marginals(d).iter))
    as(Array, as(d1), size(marginals(d))...)
end

function Base.rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure, d1::Dists.Distribution) where {T}
    mar = marginals(d)

    # Distributions doens't (yet) have the three-argument form
    elT = typeof(rand(rng, first(mar)))

    sz = size(mar)
    x = Array{elT, length(sz)}(undef, sz)
    for (j,parj) in enumerate(d.pars)
        x[j] = rand(rng, d.f(parj))
    end
    x
end

# e.g. set(Normal(μ=2)^5, params, randn(5))
function Accessors.set(d::ProductMeasure{F,A}, ::typeof(params), p::AbstractArray) where {F,A<:AbstractArray}
    set.(marginals(d), params, p)
end

function Accessors.set(d::ProductMeasure{F,A}, ::typeof(params), p) where {F,A<:AbstractArray}
    par = typeof(d.pars[1])(p)
    ProductMeasure(d.f, Fill(par, size(d.pars)))
end


text:
# ---
# title: 1161. Maximum Level Sum of a Binary Tree
# id: problem1161
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Breadth-first Search
# link: <https://leetcode.com/problems/maximum-level-sum-of-a-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, the level of its root is `1`, the level of
# its children is `2`, and so on.
#
# Return the **smallest** level `X` such that the sum of all the values of nodes
# at level `X` is **maximal**.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/03/capture.JPG)
#
#
#
#     Input: root = [1,7,0,7,-8,null,null]
#     Output: 2
#     Explanation:
#     Level 1 sum = 1.
#     Level 2 sum = 7 + 0 = 7.
#     Level 3 sum = 7 + -8 = -1.
#     So we return the level with the maximum sum which is level 2.
#
#
# **Example 2:**
#
#
#
#     Input: root = [989,null,10250,98693,-89388,null,null,null,-32127]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 104]`.
#   * `-105 <= Node.val <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# # L2-projection
#
# ![](heat_square_fluxes.png)

# ## Introduction
# This example continues from the Heat equation example, where the temperature field was determined on a square domain. In this example, we first compute the heat flux in each integration point (based on the solved temperature field) and then we do an L2-projection of the fluxes to the nodes of the mesh. By doing this, we can more easily visualize integration points quantities.

# ## Implementation
# Start by simply running the Heat equation example
include("heat_equation.jl");

# Next we define a function that computes the heat flux for each integration point in the domain. Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit conductivity $\lambda = 1 ⇒ q = - \nabla T$.
function compute_heat_fluxes(cellvalues::CellScalarValues{dim,T}, dh::DofHandler, a) where {dim,T}

    n = getnbasefunctions(cellvalues)
    cell_dofs = zeros(Int, n)
    nqp = getnquadpoints(cellvalues)

    ## Allocate storage for the fluxes to store
    q = [Vec{2,T}[] for _ in 1:getncells(dh.grid)]

    @inbounds for (cell_num, cell) in enumerate(CellIterator(dh))
        q_cell = q[cell_num]
        celldofs!(cell_dofs, dh, cell_num)
        aᵉ = a[cell_dofs]
        reinit!(cellvalues, cell)

        for q_point in 1:nqp
            q_qp = -function_gradient(cellvalues, q_point, aᵉ)
            push!(q_cell, q_qp)
        end
    end
    return q
end
#md nothing # hide

# Now call the function to get all the fluxes.
q = compute_heat_fluxes(cellvalues, dh, u);

# Next, create an L2-projector using the same interpolation as was used to approximate the temperature field. On instantiation, the projector assembles the coefficient matrix `M` and computes the Cholesky factorization of it. By doing so, the projector can be reused without having to invert `M` every time.
projector = L2Projector(cellvalues, ip, grid);

# Project the integration point values to the nodal values
q_nodes = project(q, projector);

# ## Exporting to VTK
# To visualize the heat flux, we export the projected field `q`
# to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
vtk_grid("heat_equation_flux", grid) do vtk
    vtk_point_data(vtk, q_nodes, "q")
end;

#md # ## [Plain Program](@id l2_projection-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [l2_projection.jl](l2_projection.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
export ft_setdeadmantimeout

function ft_setdeadmantimeout(ft_handle::Culong, deadmantimeout::Integer = 5000)
  ft_status = ccall((:FT_SetDeadmanTimeout, d2xx),
                     Cuint,
                     (Culong, Cuint),
                     ft_handle, deadmantimeout)
  checkstatus(ft_status)
  return nothing
end


text:
@testset "snake game env" begin
    # env = SnakeGameEnv()
    # RLBase.test_runnable!(env)
end


text:
import Base: length, size, tail, iterate, eltype, IteratorSize, IteratorEltype, haslength, SizeUnknown, @propagate_inbounds, HasEltype

# progress(minimize(f, repeat(data,10)))
# A stream (iterator) based implementation: minimize works like map
# taking a stream of args and generating a stream of func values
# except applying gradient based updates to params at each step

#"Example: `minimize(f,repeat(data,10))`"
minimize(f,d::I,a=Adam(); params=nothing) where {I} = Minimize{I}(d,f,a,params)
minimize!(x...; o...) = for x in minimize(x...; o...); end

struct Minimize{I}; data::I; func; algo; params; end

IteratorSize(::Type{Minimize{I}}) where {I} = IteratorSize(I)
IteratorEltype(::Type{<:Minimize}) = Base.EltypeUnknown()
length(m::Minimize) = length(m.data)
size(m::Minimize) = size(m.data)

@propagate_inbounds function iterate(m::Minimize, s...)
    next = iterate(m.data, s...)
    next === nothing && return nothing
    (args, s) = next
    y = @diff m.func(args...)
    for x in (m.params === nothing ? params(y) : m.params)
        if x.opt === nothing
            x.opt = clone(m.algo)
        end
        update!(x, grad(y,x))
    end
    #Returning the tape is risky: prevents gc and causes >2x slowdown at memory limit
    #return (y,s)
    return (value(y),s)
end

"""
    converge(itr; alpha=0.1)

Return an iterator which acts exactly like `itr`, but quits when values from `itr` stop
decreasing. `itr` should produce numeric values.

It can be used to train a model with the data cycled:

    progress!(converge(minimize(model,cycle(data))))

`alpha` controls the exponential average of values to detect convergence. Here is how
convergence is decided:

    p = x - avgx
    avgx = c.alpha * x + (1-c.alpha) * avgx
    avgp = c.alpha * p + (1-c.alpha) * avgp
    avgp > 0.0 && return nothing

`converge!(...)` is equivalent to `(for x in converge(...) end)`, i.e.  iterates over the
object created by `converge(...)` and returns `nothing`.

"""
converge(iter::I; alpha=0.1) where {I} = Converge{I}(iter, alpha)
converge!(x...; o...) = for x in converge(x...; o...); end

struct Converge{I}; iter::I; alpha::Float64; end

# Converge is large Filter, does not have known size
# length(c::Converge) = length(c.iter)
# size(c::Converge) = size(c.iter)
eltype(c::Converge) = eltype(c.iter)
IteratorEltype(::Type{Converge{I}}) where {I} = IteratorEltype(I)
IteratorSize(::Type{<:Converge}) = SizeUnknown()

@propagate_inbounds function iterate(c::Converge, s=(0.0,Inf))
    avgp,avgx,state = s[1],s[2],tail(tail(s))
    next = iterate(c.iter, state...)
    next === nothing && return nothing
    (item, state) = next
    x = value(item)
    if avgx == Inf; avgx = x; end
    p = x - avgx
    avgx = c.alpha * x + (1-c.alpha) * avgx
    avgp = c.alpha * p + (1-c.alpha) * avgp
    avgp > 0.0 && return nothing
    (item, (avgp, avgx, state))
end

"""
    param(array; atype)
    param(dims...; init, atype)
    param0(dims...; atype)

The first form returns `Param(atype(array))` where `atype=identity` is the default.

The second form Returns a randomly initialized `Param(atype(init(dims...)))`.  By default,
`init` is `xavier_uniform` and `atype` is `KnetArray{Float32}` if `gpu() >= 0`,
`Array{Float32}` otherwise.

The third form `param0` is an alias for `param(dims...; init=zeros)`.
"""
param,param0

# TODO: Knet.Param <: AutoGrad.Tracked as a separate type?
param(x::Union{Array,KnetArray}; atype=identity) = Param(atype(x))
param(d...; init=xavier_uniform, atype=atype())=Param(atype(init(d...)))
param0(d...; atype=atype())=param(d...; init=zeros, atype=atype)
atype()=(gpu() >= 0 ? KnetArray{Float32} : Array{Float32})

### DEPRECATED:

"""
    train!(model, data; loss, optimizer, callback, o...)

Train a model with given data. This function is deprecated, please use `sgd`, `adam` etc.

* `model`: A callable object. `model(x; o...)` should return a prediction. `params(model)`
   will automatically iterate over model parameters.
* `data`: An iterator. `for (x,y) in data` should iterate over input-output pairs.
* `loss=nll`: A loss function, called with `J = @diff loss(model,x,y; o...)`.
* `optimizer=Adam()`: An optimizer object that will be copied for each parameter and used by
  `[update!]`(@ref).
* `callback`: To facilitate reporting and termination, a callback function is called before
   every update with `callback(J)`. Training continues if the return value is true, terminates
   if it is false. By default training will end after one pass over the data.
* Other keyword arguments `(o...)` will be passed to `loss` and possibly by `loss` to `model`.
"""
function train!(model, data; loss=nll, optimizer=Adam(), callback=epochs(data,1), o...)
    @warn "train! is deprecated, use sgd!, adam! etc. instead." maxlog=1
    ps = params(model)
    for param in ps
        if param.opt === nothing
            param.opt = clone(optimizer)
        end
    end
    while true
        for (x,y) in data
            J = @diff loss(model,x,y; o...)
            if !callback(J)
                return
            end
            for param in ps
                g = grad(J,param)
                update!(value(param),g,param.opt)
            end
        end
    end
end

# """
# Pre-defined callback function constructors:

# * converge(): Trains until convergence
# * updates(n): Stops after n updates
# * epochs(data,n): Trains for n epochs, equivalent to updates(n*length(data))
# """
# converge, updates, epochs

function converge(alpha::Number = 0.001)
    avgx = Inf
    avgp = 0.0
    # prog = Progress()
    function callback(x)
        x = value(x)
        if avgx == Inf; avgx = x; end
        p = x - avgx
        avgx = alpha * x + (1-alpha) * avgx
        avgp = alpha * p + (1-alpha) * avgp
        # display_progress!(prog, x)
        return avgp <= 0.0
    end
    return callback
end

function updates(n)
    # p = Progress(n)
    function callback(x)
        # display_progress!(p, value(x))
        n -= 1
        return n > 0
    end
end

epochs(d,n)=updates(n*length(d))

# # Iterator version:
# "Example: `progress!(train(f,repeat(data,10)))`"
# train(pred, data::I; loss=nll, optimizer=Adam(), callback=nothing, params=nothing, kw...) where {I} = Train{I}(data,pred,loss,optimizer,callback,params,kw,Any)

# # Let's not overwrite old train! for backward compatibility
# #train!(x...; o...) = for x in train(x...; o...); end

# struct Train{I}; data::I; pred; loss; optimizer; callback; params; kw; eltype; end

# length(c::Train) = length(c.data)
# size(c::Train) = size(c.data)
# eltype(c::Train) = (c.eltype === Any ? (c.eltype=typeof(@diff c.loss(c.pred,first(c.data)...;c.kw...))) : c.eltype)
# IteratorSize(::Type{Train{I}}) where {I} = IteratorSize(I)
# IteratorEltype(::Type{<:Train}) = Base.HasEltype()

# @propagate_inbounds function iterate(m::Train, s...)
#     next = iterate(m.data, s...)
#     next === nothing && return nothing
#     (args, s) = next
#     y = @diff m.loss(m.pred, args...; m.kw...)
#     m.callback !== nothing && !m.callback(y) && return nothing
#     for x in (m.params === nothing ? params(y) : m.params)
#         if x.opt === nothing
#             x.opt = clone(m.optimizer)
#         end
#         update!(x, grad(y,x))
#     end
#     return (value(y),s)
# end

### DEAD CODE:

    ## This may be slightly faster but risky if active params change
    # if m.params === nothing
    #     m.params = params(y, m.algo)
    # end
    # for x in m.params
    #     update!(x, grad(y,x))
    # end

# function AutoGrad.params(y::AutoGrad.Tape, optimizer=nothing)
#     p = Param[]
#     for node in y.list
#         x = node.Value
#         if isa(x, Param)
#             if x.opt === nothing && optimizer !== nothing
#                 x.opt = clone(optimizer)
#             end
#             push!(p, x)
#         end
#     end
#     return p
# end

# # Simpler and more flexible alternative to train!
# # Does not care where model ends loss begins or where params are
# # data may consist of tuples of any number of args
# # Epochs can be set by data iterator (convergence cannot)
# function minimize!(func, data, optimizer=Adam())
#     for args in data
#         y = @diff func(args...)
#         for node in y.list      # breaks abstraction
#             x = node.Value
#             if isa(x, Param)
#                 g = grad(y,x)
#                 if x.opt === nothing; x.opt = clone(optimizer); end
#                 update!(x.value, g, x.opt)
#             end
#         end
#     end
# end

# "Returns an iterator over Params on Tape."
# struct Params; tape::AutoGrad.Tape; end

# eltype(::Type{Params}) = Param
# IteratorEltype(::Type{Params}) = HasEltype()
# IteratorSize(::Type{Params}) = SizeUnknown()

# @propagate_inbounds function iterate(p::Params, s::Int=1)
#     next = iterate(p.tape.list, s)
#     while next !== nothing
#         (n,s) = next
#         if isa(n.Value,Param)
#             return (n.Value,s)
#         end
#         next = iterate(p.tape.list, s)
#     end
#     nothing
# end

# # Alternative simpler definition:
# params(t::Tape) = (n.Value for n in t.list if n.Value isa Param)

### DEAD CODE

### Issues:
# + What if we call train multiple times, and don't want to use the optimizers?
# x Do we want parameter initialization as well? init and opt init should happen once.
# - Recording losses with different loss functions.
# x What info does the callback need?
# - Are we doing anything other than pushing kwargs from train to Train?
# - What if we want convergence in trnloss or convergence in devloss? Return earlier (best) model?
# + How do we easily measure epochs?
# + ProgressMeter both in time mode and converge mode.
# + Printing loss with ProgressMeter seems difficult.
# + Frequency of progress updates and loss calculations?

# + Keyword argument problem:
# - optimizer, loss, model can all take keyword args; how do we specify them through train?
# + We can give a constructed optimizer and clone it for each param.
# ? We don't call model directly, only through loss (because it may need model params for regularization).
# ? So we pass all unrecognized kwargs to loss and let it sort out.

# x What to pass to the callback:
# x model, data, loss, optimizer and (o...) are all available to the caller. No need to pass to callback.
# x The only things that are not available are J,x,y. I can't think of a use for x,y.
# x That leaves J. I considered passing value(J), however that prevents the callback from looking at gradients.
# + (e.g. for reporting the gradient norms), so I decided to pass back J as is.

# x We assume a model is just a callable object (https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects-1)
# x model(x) will give us a prediction, and params(model) will iterate over the parameters.

# + 20190105: Do we even need to assume this? train! can simply look at the Tape to find the
# + parameters! In that case optimizers would need to be set elsewhere.

# x use HasLength after data
# x converge may not have length?
# + first efficiency of iterating y.list
# x separate Param in Knet?

# + write train(model,data) iterator style
# + fix update between display_progress and progress
# x progress should handle HasLength
# + use tape iter in train
# + write params tape as iterator
# - check regularization:
#     do we need opt args?
#     regularizer as parametric fn?
#     regularizer as part of optimizer?
# - write docs
# x use throttle?

# + use cycle for repeat
# + use take for updates: take(cycle(data),n)
# + shuffling during repeats?
# x filter for params(tape) and converge?
# + make params an optional argument


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

@doc raw"""AllowedHostPath defines the host volume conditions that will be enabled by a policy for pods to use. It requires the path prefix to be defined. Deprecated: use AllowedHostPath from policy API Group instead.

    IoK8sApiExtensionsV1beta1AllowedHostPath(;
        pathPrefix=nothing,
        readOnly=nothing,
    )

    - pathPrefix::String : pathPrefix is the path prefix that the host volume must match. It does not support `*`. Trailing slashes are trimmed when validating the path prefix with a host path.  Examples: `/foo` would allow `/foo`, `/foo/` and `/foo/bar` `/foo` would not allow `/food` or `/etc/foo`
    - readOnly::Bool : when set to true, will allow host volumes matching the pathPrefix only if all volume mounts are readOnly.
"""
mutable struct IoK8sApiExtensionsV1beta1AllowedHostPath <: SwaggerModel
    pathPrefix::Any # spec type: Union{ Nothing, String } # spec name: pathPrefix
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiExtensionsV1beta1AllowedHostPath(;pathPrefix=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1AllowedHostPath, Symbol("pathPrefix"), pathPrefix)
        setfield!(o, Symbol("pathPrefix"), pathPrefix)
        validate_property(IoK8sApiExtensionsV1beta1AllowedHostPath, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiExtensionsV1beta1AllowedHostPath

const _property_map_IoK8sApiExtensionsV1beta1AllowedHostPath = Dict{Symbol,Symbol}(Symbol("pathPrefix")=>Symbol("pathPrefix"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiExtensionsV1beta1AllowedHostPath = Dict{Symbol,String}(Symbol("pathPrefix")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1AllowedHostPath }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1AllowedHostPath))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1AllowedHostPath }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1AllowedHostPath[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1AllowedHostPath }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1AllowedHostPath[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1AllowedHostPath)
    true
end

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


text:
abstract type RandomSearch end
include("purerandomsearch.jl")
include("simulatedannealing.jl")


text:
# default root and comm for this package
"Default MPI root for `RMPI`."
const mpi_root = zero(Int32)

"Default MPI communicator for `RMPI`."
mpi_comm() = MPI.COMM_WORLD

"""
    mpi_size(comm = mpi_comm())
Size of MPI communicator.
"""
mpi_size(comm = mpi_comm()) = MPI.Comm_size(comm)

"""
    next_mpiID()
Produce a new ID number for MPI distributed objects. Uses an internal counter.
"""
next_mpiID

let mpiID::Int = 0
    global next_mpiID
    next_mpiID() = (mpiID +=1; mpiID)
end

"""
    mpi_rank(comm = mpi_comm())
Return the current MPI rank.
"""
mpi_rank(comm = mpi_comm()) = MPI.Comm_rank(comm)

"""
    is_mpi_root(root = mpi_root)
Returns `true` if called from the root rank
"""
is_mpi_root(root = mpi_root) = mpi_rank() == root

"""
    @mpi_root expr
Evaluate expression only on the root rank.
Extra care needs to be taken as `expr` *must not* contain any code that involves
syncronising MPI operations, i.e. actions that would require syncronous action
of all MPI ranks.

Example:
```julia
wn = walkernumber(dv)   # an MPI syncronising function call that gathers
                        # information from all MPI ranks
@mpi_root @info "The current walker number is" wn # print info message on root only
```
"""
macro mpi_root(args...)
    :(is_mpi_root() && $(esc(args...)))
end

macro mpi_on(p, args...)
    :(mpi_rank() == $p && $(esc(args...)))
end

"""
    mpi_barrier(comm = mpi_comm())
The MPI barrier with optional argument. MPI syncronizing.
"""
mpi_barrier(comm = mpi_comm()) = MPI.Barrier(comm)

"""
    targetrank(key, np)
Compute the rank where the `key` belongs.
"""
targetrank(key::Union{Integer,AbstractFockAddress}, np) = hash(key, hash(1)) % np
targetrank(pair::Pair, np) = targetrank(pair[1], np)

"""
    mpi_combine_walkers!(target, source, [strategy])
Distribute the entries of `source` to the `target` data structure such that all
entries in the `target` are on the process with the correct mpi rank
as controlled by [`targetrank()`](@ref).
MPI syncronizing.

Note: the [`storage`](@ref) of the `source` is communicated rather than the `source` itself.
"""
function mpi_combine_walkers!(dtarget::MPIData, source::AbstractDVec)
    ltarget = localpart(dtarget)
    empty!(ltarget) # start with empty slate
    strategy = dtarget.s
    mpi_combine_walkers!(ltarget, storage(source), strategy)
end

# This function is just a wrapper that makes allreduce treat a SVector as a scalar
function Rimu.sort_into_targets!(dtarget::MPIData, ws::NTuple{NT,W}, statss) where {NT,W}
    # multi-threaded MPI version
    # should only ever run on thread 1
    @assert Threads.threadid() == 1 "`sort_into_targets!()` is running on `threadid()` == $(Threads.threadid()) instead of 1!"
    lwm = ws[1]
    for i in 2:NT # combine new walkers generated from different threads
        add!(lwm, ws[i])
    end
    mpi_combine_walkers!(dtarget,lwm) # combine walkers from different MPI ranks
    stats = Rimu.combine_stats(statss) # combine stats from all threads
    res_stats = MPI.Allreduce(Rimu.MultiScalar(stats), +, dtarget.comm)
    return dtarget, ws, res_stats
end
function Rimu.sort_into_targets!(dtarget::MPIData, w::AbstractDVec, stats)
    # single threaded MPI version
    mpi_combine_walkers!(dtarget,w) # combine walkers from different MPI ranks
    res_stats = MPI.Allreduce(Rimu.MultiScalar(stats), +, dtarget.comm)
    return dtarget, w, res_stats
end

"""
    sync_cRandn(md::MPIData)
Generate one random number with [`cRandn()`](@ref) in a synchronous way such
that all MPI ranks have the same random number.
The argument is ignored unless it is of type `MPIData`, in which case a random
number from the root rank is broadcasted to all MPI ranks. MPI syncronizing.
"""
function ConsistentRNG.sync_cRandn(md::MPIData)
    MPI.bcast(cRandn(), md.root, md.comm)
end

"""
    ConsistentRNGs.check_crng_independence(dv::MPIData)
Does a sanity check to detect dependence of random number generators across
all MPI ranks. Returns the size of the combined RNG state,
i.e. `mpi_size()*Threads.nthreads()*fieldcount(ConsistentRNG.CRNG)`.
MPI syncronizing.
"""
ConsistentRNG.check_crng_independence(dv::MPIData) = _check_crng_independence(dv.comm)

function _check_crng_independence(comm::MPI.Comm) # MPI syncronizing
    # get vector of threaded RNGs on this rank
    crngs = ConsistentRNG.CRNGs[]
    # extract the numbers that make up the state of the RNGs and package into
    # an MPI-suitable buffer
    statebuffer = [getfield(rng,i) for rng in crngs for i in 1:fieldcount(typeof(rng))]
    # gather from all ranks
    combined_statebuffer = MPI.Allgather(statebuffer, comm)  # MPI syncronizing
    # check independence
    @mpi_root @assert union(combined_statebuffer) == combined_statebuffer "Dependency in parallel rngs detected"

    @assert length(ConsistentRNG.CRNGs[]) == Threads.nthreads() "Number of CNRGs should be equal to nthreads()"
    return length(combined_statebuffer)
end

"""
    mpi_seed_CRNGs!(seed = rand(Random.RandomDevice(), UInt))
Re-seed the random number generators in an MPI-safe way. If seed is provided,
the random numbers from [`cRand()`](@ref) will follow a deterministic sequence.

Independence of the random number generators on different MPI ranks is achieved
by adding `hash(mpi_rank())` to `seed`.
"""
function mpi_seed_CRNGs!(seed = rand(Random.RandomDevice(), UInt))
    rngs = seedCRNG!(seed + hash(mpi_rank()))
    _check_crng_independence(mpi_comm())
    return rngs
end

"""
    mpi_allprintln(args...)
Print a message to `stdout` from each rank separately, in order. MPI synchronizing.
"""
function mpi_allprintln(args...)
    for i in 0:(mpi_size() - 1)
        if mpi_rank() == i
            println("[ rank ", lpad(i, length(string(mpi_size() - 1))), ": ", args...)
            flush(stdout)
        end
        mpi_barrier()
    end
end


text:
"""
    BagCount{T <: AbstractAggregation}

A wrapper type that when called applies the [`AbstractAggregation`](@ref) stored in it,
and appends one more element containing bag size after ``x ↦ \\log(x + 1)`` transformation to the result.

Used as a functor:

    (bc::BagCount)(x, bags[, w])

where `x` is either `Missing`, `AbstractMatrix` or [`ArrayNode`](@ref),
`bags` is [`AbstractBags`](@ref) structure and optionally `w` is an `AbstractVector` of weights.

# Examples
```jldoctest
julia> x = Float32[0 1 2; 3 4 5]
2×3 Matrix{Float32}:
 0.0  1.0  2.0
 3.0  4.0  5.0

julia> b = bags([1:1, 2:3])
AlignedBags{Int64}(UnitRange{Int64}[1:1, 2:3])

julia> a = vcat(SegmentedMean(2), SegmentedMax(2))
AggregationStack:
 SegmentedMean(ψ = Float32[0.0, 0.0])
 SegmentedMax(ψ = Float32[0.0, 0.0])

julia> a(x, b)
4×2 Matrix{Float32}:
 0.0  1.5
 3.0  4.5
 0.0  2.0
 3.0  5.0

julia> BagCount(a)(x, b)
5×2 Matrix{Float32}:
 0.0       1.5
 3.0       4.5
 0.0       2.0
 3.0       5.0
 0.693147  1.09861
```

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`SegmentedSum`](@ref),
    [`SegmentedMax`](@ref), [`SegmentedMean`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
struct BagCount{T <: AbstractAggregation}
    a::T
end

Flux.@functor BagCount

_bagcount(T, bags) = permutedims(log.(one(T) .+ length.(bags)))
ChainRulesCore.@non_differentiable _bagcount(T, bags)

function (bc::BagCount)(x::Maybe{AbstractArray}, bags::AbstractBags, args...)
    o1 = bc.a(x, bags, args...)
    o2 = _bagcount(eltype(o1), bags)
    vcat(o1, o2)
end

function Base.show(io::IO, m::MIME"text/plain", @nospecialize(bc::BagCount{T})) where T
    print(io, "BagCount(", repr(m, bc.a))
    print(io, T <: AggregationStack ? "\n)" : ")")
end

function Base.show(io::IO, @nospecialize(bc::BagCount))
    print(io, "BagCount")
    if !get(io, :compact, false)
        print(io, "(", bc.a, ")")
    end
end


text:
"""
    Compressor

Representation of compressor supprted by blosc
"""
struct Compressor
    name ::Symbol
    libname ::String
    code ::Int8
    version ::String
end

function fill_available_compressors!(dict::Dict{Symbol, Compressor})
    names = Lib.blosc_list_compressors()
    for n in names
        code, complib, version = Lib.blosc_get_complib_info(n)
        code < 0 && continue
        dict[Symbol(n)] = Compressor(Symbol(n),complib, code, version)
    end
end

function _get_default_compressor_name()
    name = Lib.blosc_compcode_to_compname(0)
    return Symbol(name)
end

"""
    available_compressors()::Dict{Symbol, Compressor}

Get dictionary of available compressors
"""
available_compressors()::Dict{Symbol, Compressor} = _available_compressors

"""
    available_compressors_names()::Vector{Symbol}

Get vector of available compressors names
"""
available_compressors_names() = collect(keys(_available_compressors))

"""
    has_compressor(name::Symbol)

Check if compressor `name` available
"""
has_compressor(name::Symbol) = haskey(_available_compressors, name)

"""
    compressor(name::Symbol)

Get compressor by name
"""
function compressor_by_name(name::Symbol)
    !has_compressor(name) && throw(KeyError(name))
    _available_compressors[name]
end

"""
    default_compressor_name()

Name of default compressor
"""
default_compressor_name() = _default_compressor_name

"""
    default_compressor()

Default compressor
"""
default_compressor() = _available_compressors[_default_compressor_name]


text:
abstract type AbstractSonum <: AbstractFloat end

primitive type Sonum8 <: AbstractSonum 8 end
primitive type Sonum16 <: AbstractSonum 16 end


text:
"Kuramoto model parameter structure."
struct kp
    G::Float64
    σ::Float64
    w::Array{Float64,2}
    l::Array{Float64,2}
end

"Kuramoto order parameter"
kop(sol) = [abs(mean(exp.(1.0im .* u_t))) for u_t in sol.u];

"Kuramoto dfun."
function kf(du,u,h,p,t)
    G = p.G/length(u)
    du .= 2*2*π/1000
    for j=1:length(u)
        for i=1:length(u)
            if i!=j
                uⱼ = h(p, t - p.l[i,j]; idxs=j)[1]
                du[i] += G * p.w[i,j] * sin(uⱼ - u[i])
            end
        end
    end
end

"Kuramoto noise function."
function kg(du,u,h,p,t)
    du .= sqrt(2*p.σ)
end

"Kuramoto history function."
function kh(p,t;idxs=nothing)
    nn = size(p.w,1)
    x = ((1:nn)/nn*2*pi .+ t) .* 2 .* pi
    collect(idxs==nothing ? x : x[idxs])
end


text:
function besselj_zero_approx1(ν::Integer, m::Integer)
    μ = 4.0*ν^2
    a = (m+ν/2-1/4)*π
    a8 = 1/(8a)
    a82 = a8^2
    j = (a -(μ-1)*a8*(1 + a82*(4*(7μ-31)/3 +a82*( 32*(83*μ^2-982*μ+3779)/15 +
         a82*(64*(6949*μ^3-153855*μ^2+1585743*μ-6277237)/105))))
        )
end

function besseljprime_zero_approx1(ν::Integer, m::Integer)
    #if (ν==0 && m==1)
    #    return 0.0
    #end
    if ν==0
        m += 1
    end
    μ = 4*ν^2
    b = (m+ν/2-3/4)*π
    b8 = 1/(8*b)
    b82 = b8^2
    j = (b - b8*(μ+3 + b82*(4*(7*μ^2+82*μ-9)/3 +b82*(32*(83*μ^3+2075*μ^2-3039*μ+3537)/15 +
             b82*(64*(6949*μ^4+296492*μ^3-1248002*μ^2+7414380*μ-5853627)/105)))))
end

function airy_zero_approx(m::Integer)
    if m<=10
        z =[-2.338107410459767,
            -4.087949444130970,
            -5.520559828095552,
            -6.786708090071759,
            -7.944133587120853,
            -9.022650853340980,
           -10.04017434155809,
           -11.00852430373326,
           -11.93601556323626,
           -12.82877675286576][m]
    else
        t = 3/8*π*(4*m-1)
        t2 = t^(-2)
        z = -t^(2/3)*(1+t2*(5/48+t2*(-5/36+t2*(77125/82944+t2*(
            -108056875/6967296+ t2*162375596875/334430208)))))
    end
    z
end

function airyprime_zero_approx(m::Integer)
    if m<=10
        z = [-1.0187929716474710890,
             -3.2481975821798365379,
             -4.8200992111787356394,
             -6.1633073556394865476,
             -7.3721772550477701771,
             -8.4884867340197221329,
             -9.5354490524335474707,
            -10.527660396957407282,
            -11.475056633480245295,
            -12.384788371845747325][m]
    else
        t = 3/8*π*(4*m-3)
        t2 = t^(-2)
        z = -t^(2/3)*(1+t2*(-7/48+t2*(35/288+t2*(-181223/207360+t2*(
            18683371/1244160 - t2*9114588436/191102976)))))
    end
    z
end

function besselj_zero_approx2(ν::Integer, m::Integer; prime::Bool=false)
    if prime
        ζ = ν^(-2/3)*airyprime_zero_approx(m)
    else
         ζ = ν^(-2/3)*airy_zero_approx(m)
    end
    y = 2/3*(-ζ)^(3/2)
    if y>100000
        x = π/2
    elseif y<1
        p = (3*y)^(1/3)
        p2 = p^2
        x = p*(1+p2*(-2/15+p2*(3/175+p2*(-2/1575))))
    else
        p = 1/(y+π/2)
        p2 = p^2
        x = π/2 - p*(1+p2*(2/3+p2*(13/15+p2*(146/105+p2*(781/315+p2*16328/3465)))))
    end
    x2 = (y+x)^2
    r = (x-atan(x+y))/x2
    x = x - (1+x2)*r*(1+r/(x+y))

    z = 1/cos(x)
    h = sqrt(ζ*(1-z^2))
    if prime
        g1 = z/ζ*h*(7/(48*ζ) + h*(7/(z^2-1)+9
        )/24)
        j = ν*z + g1/ν
    else
        f1 = -z/ζ*h*(5/(48*ζ) + h*(5/(z^2-1)+3)/24)
        j = ν*z + f1/ν
    end
    j
end

function besselj_zero_approx(ν::Integer, m::Integer)
  if m>=ν
     j = besselj_zero_approx1(ν, m)
  else
     j = besselj_zero_approx2(ν, m)
  end
  j
end

function besseljprime_zero_approx(ν::Integer, m::Integer)
  if m>=ν
        j = besseljprime_zero_approx1(ν, m)
  else
        j = besselj_zero_approx2(ν, m, prime=true)
  end
  j
end

function besselj_zero_iter(ν::Integer, z::AbstractFloat)
    T = typeof(z)
    ɛ = eps(T)
    for i = 1:200
        J = besselj(ν, z)
        if abs(J) < 1000ɛ
            break
        end
        Jprime = besselj(ν-1, z) - ν*J/z
        z -= J/Jprime
        if i==200
            println("200 iterations, res=",abs(J))
        end
    end
    return z
end

function besseljprime_zero_iter(ν::Integer, z::AbstractFloat)
    T = typeof(z)
    ɛ = eps(T)
    for i = 1:200
        J = besselj(ν, z)
        J1 = besselj(ν-1, z) - ν*J/z
        if abs(J1) < 1000ɛ
            break
        end
        J2 = -J1/z + ((ν/z)^2-1)*J
        z -= J1/J2
        if i==200
            println("200 iterations, res=",abs(J1))
        end
    end
    return z
end

besselj_zero(ν, m; T=Float64) = besselj_zero_iter(ν, convert(T, besselj_zero_approx(ν, m)))

besseljprime_zero(ν::Integer, m::Integer; T=Float64) =
    besseljprime_zero_iter(ν, convert(T, besseljprime_zero_approx(ν, m)))


text:
using Distributed; addprocs(1)
@everywhere using RemoteObjects
using Sockets
using Test

@everywhere begin
    struct MyStruct
        x
    end
    mycall(m::MyStruct) = m.x
end
@everywhere mimic(MyStruct; force=true)

@testset "Basics" begin

r = RemoteObject(MyStruct, 1)
@test r isa RemoteObject{MyStruct}
x = mycall(r)
@test x isa RemoteObject{Int}
@test r.uuid != x.uuid
@test Base.fetch(x) == 1
iob = IOBuffer()
Base.show(iob, r)
str = String(take!(iob))
@test str == "RemoteObject (worker 2):\nMyStruct(1)"

r = @remote MyStruct(2)
@test r isa RemoteObject{MyStruct}
x = Base.fetch(r)
@test x.x == 2

end

@testset "Remote Server" begin

port = rand(20000:30000)
tsk = @async run_server(port)
yield()
conn = connect_remote(Sockets.localhost, port)

# FIXME: Set timedwait in case server hangs

r = @remote conn rand(Float32)
@test r isa RemoteObject{Float32}
x = Base.fetch(r)
@test x isa Float32

@test_throws UndefVarError @remote conn lol()

sleep(1)

end


text:
using Documenter
using ReadStatTables
using CategoricalArrays

DocMeta.setdocmeta!(ReadStatTables, :DocTestSetup, :(using ReadStatTables, CategoricalArrays))

makedocs(
    modules = [ReadStatTables],
    format = Documenter.HTML(
        canonical = "https://junyuan-chen.github.io/ReadStatTables.jl/stable/",
        prettyurls = get(ENV, "CI", nothing) == "true",
        edit_link = "main"
    ),
    sitename = "ReadStatTables.jl",
    authors = "Junyuan Chen",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting Started" => "man/getting-started.md",
            "Table Interface" => "man/table-interface.md",
            "Value Labels" => "man/value-labels.md",
            "Date and Time Values" => "man/date-and-time-values.md"
        ],
        "About" => [
            "License" => "license.md"
        ]
    ],
    workdir = joinpath(@__DIR__, "..")
)

deploydocs(
    repo = "github.com/junyuan-chen/ReadStatTables.jl.git",
    devbranch = "main"
)