text:
#  ======================================================
#  Constructors
#  ======================================================

#=
mutable struct SimpleValueDiGraph{V <: Integer, E_VAL, E_VAL_C <: Union{Nothing, EdgeValContainer{E_VAL}}, RE_VAL_C <: Union{Nothing, EdgeValContainer{E_VAL}} } <: AbstractSimpleValueGraph{V, E_VAL}
    ne::Int
    fadjlist::Adjlist{V}
    badjlist::Adjlist{V}
    edgevals::E_VAL_C
    redgevals::RE_VAL_C
end
=#

mutable struct SimpleValueDiGraph{V,
                                  E_VAL,
                                  E_VAL_C <: Union{Nothing, EdgeValContainer{E_VAL}},
                                  RE_VAL_C <: Union{Nothing, EdgeValContainer{E_VAL}},
                                 } <: AbstractSimpleValueGraph{V, E_VAL}
    ne::Int
    fadjlist::Adjlist{V}
    badjlist::Adjlist{V}
    edgevals::E_VAL_C
    redgevals::RE_VAL_C
end

const ConstValueDiGraph{V, E_VAL} = SimpleValueDiGraph{V, E_VAL, Nothing, Nothing}
const OutValueDiGraph{V, E_VAL} = SimpleValueDiGraph{V, E_VAL, <: EdgeValContainer{E_VAL}, Nothing}
const InOutValueDiGraph{V, E_VAL} = SimpleValueDiGraph{V, E_VAL, <: EdgeValContainer{E_VAL}, <: EdgeValContainer{E_VAL}}

function SimpleValueDiGraph(nv::V,
                            E_VAL::Type=default_edgeval_type;
                            const_edgevals::Bool=false,
                            reverse_edgevals::Bool=false
                           ) where {V <: Integer}
    fadjlist = Adjlist{V}(nv)
    badjlist = Adjlist{V}(nv)
    if const_edgevals
        return ConstValueDiGraph(0, fadjlist, badjlist, nothing, nothing)
    end
    edgevals = create_edgeval_list(nv, E_VAL)
    if !reverse_edgevals
        return OutValueDiGraph{V, E_VAL, typeof(edgevals)}(0,
                                                          fadjlist, badjlist,
                                                          edgevals, nothing)
    end
    redgevals = create_edgeval_list(nv, E_VAL)
    return InOutValueDiGraph{V, E_VAL, typeof(edgevals), typeof(redgevals)}(0,
                                                                            fadjlist,
                                                                            badjlist,
                                                                            edgevals,
                                                                            redgevals)
end

SimpleValueDiGraph{V}(nv::Integer, E_VAL::Type=default_edgeval_type; args...) where {V} = SimpleValueDiGraph(V(nv), E_VAL, args...)
SimpleValueDiGraph{V, E_VAL}(nv::Integer, args...) where {V, E_VAL} = SimpleValueDiGraph(V(nv), E_VAL, args...)

# TODO rewrite for tuples and named tuples
# TODO const edgevals & reverse edgevals
function SimpleValueDiGraph(g::SimpleDiGraph{V},
                          E_VAL::Type=default_edgeval_type;
                          const_edgevals::Bool=false,
                          reverse_edgevals::Bool=false,
                          edgeval_initializer = (s, d) -> default_edgeval(E_VAL)) where {V}
    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)
    badjlist = deepcopy_adjlist(g.badjlist)
    if const_edgevals
        return ConstValueDiGraph{V, E_VAL}(ne(g), fadjlist, nothing)
    end
    edgevals = Vector{Vector{E_VAL}}(undef, n)
    E_VAL_C = typeof(edgevals)
    for s in Base.OneTo(n)
        edgevals[s] = [edgeval_initializer(V(s), d) for d in fadjlist[s]]
    end
    if !reverse_edgevals
        return OutValueDiGraph{V, E_VAL, E_VAL_C}(ne(g),
                                                  fadjlist,
                                                  badjlist,
                                                  edgevals,
                                                  nothing)
    end

    redgevals = E_VAL_C( T(undef, n) for T in E_VAL_C.types )
    # TODO initialize redgevals

    InOutValueDiGraph{V, E_VAL, E_VAL_C,E_VAL_C}(ne(g),
                                                 fadjlist,
                                                 badjlist,
                                                 edgevals,
                                                 redgevals)
end

# TODO this function has some issues with typesafety
# TODO weights are not symmetric
# TODO const edgevals & reverse edgevals
function SimpleValueDiGraph(g::SimpleDiGraph{V},
                          E_VAL::Type{<:TupleOrNamedTuple};
                          const_edgevals::Bool=false,
                          reverse_edgevals::Bool=false,
                          edgeval_initializer = (s, d) -> default_edgeval(E_VAL)) where {V}
    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist) # TODO deepcopy seems not be typesave
    badjlist = deepcopy_adjlist(g.fadjlist)
    if const_edgevals
        return ConstValueDiGraph{V, E_VAL}(ne(g), fadjlist, nothing)
    end

    E_VAL_C = edgevals_container_type(Val(E_VAL))
    edgevals = E_VAL_C( T(undef, n) for T in E_VAL_C.types )
    for s in Base.OneTo(n)
        for (i, T) in enumerate(E_VAL.types)
            edgevals[i][s] = Vector{T}(undef, length(fadjlist[s]))
        end
    end
    for s in Base.OneTo(n)
        for (j, d) in enumerate(fadjlist[s])
            w = edgeval_initializer(V(s), d)
            for (i, T) in enumerate(E_VAL.types)
                edgevals[i][s][j] = w[i]
            end
        end
    end
    if !reverse_edgevals
        return OutValueDiGraph{V, E_VAL, E_VAL_C}(ne(g),
                                                  fadjlist,
                                                  badjlist,
                                                  edgevals,
                                                  nothing)
    end

    redgevals = E_VAL_C( T(undef, n) for T in E_VAL_C.types )
    # TODO initialize redgevals

    return InOutValueDiGraph{V, E_VAL, E_VAL_C, E_VAL_C}(ne(g),
                                                        fadjlist,
                                                        badjlist,
                                                        edgevals,
                                                        redgevals)
end

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

function add_edge!(g::InOutValueDiGraph{V, E_VAL},
                   s::Integer,
                   d::Integer,
                   value=default_edgeval(E_VAL)) where {V, E_VAL}
    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 value, but return false
        set_value_for_index!(g.edgevals, s, index, value)
        return false
    end

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

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    set_value_for_index!(g.redgevals, d, index, value)
    return true # edge successfully added
end

function add_edge!(g::OutValueDiGraph{V, E_VAL},
                   s::Integer,
                   d::Integer,
                   value=default_edgeval(E_VAL)) where {V, E_VAL}
    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 value, but return false
        set_value_for_index!(g.edgevals, s, index, value)
        return false
    end

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

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

function add_edge!(g::ConstValueDiGraph,
                   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 if index <= length(list) && list[index] == d
        # edge already there, return false
        return false
    end

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

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

add_edge!(g::SimpleValueDiGraph, e::SimpleEdge)      = add_edge!(g, src(e), dst(e))
add_edge!(g::Union{OutValueDiGraph, InOutValueDiGraph}, e::SimpleEdge, u)   = add_edge!(g, src(e), dst(e), u)
add_edge!(g::Union{OutValueDiGraph, InOutValueDiGraph}, e::SimpleValueEdge) = add_edge!(g, src(e), dst(e), val(e))

function rem_edge!(g::InOutValueDiGraph, 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_value_for_index!(g.edgevals, s, index)

    g.ne -= 1

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

    return true
end

function rem_edge!(g::OutValueDiGraph, 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_value_for_index!(g.edgevals, s, index)

    g.ne -= 1

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

    return true
end

function rem_edge!(g::ConstValueDiGraph, 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)

    g.ne -= 1

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

    return true
end

rem_edge!(g::SimpleValueDiGraph, e::SimpleEdge) = rem_edge!(g, src(e), dst(e))
rem_edge!(g::Union{OutValueDiGraph, InOutValueDiGraph}, e::SimpleValueEdge) = rem_edge!(g, src(e), dst(e))

# TODO rem_vertex!, rem_vertices!

function has_edge(g::SimpleValueDiGraph, 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 LightGraphs.insorted(d, list_fadj)
end

function has_edge(g::SimpleValueDiGraph, s::Integer, d::Integer, value)
    verts = vertices(g)
    E_VAL = edgeval_type(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds return (index <= length(list) && list[index] == d && value_for_index(g.edgevals, E_VAL, s, index) == value)
end

has_edge(g::SimpleValueDiGraph, e::SimpleEdge)      = has_edge(g, src(e), dst(e))
has_edge(g::SimpleValueDiGraph, e::SimpleEdge, u)   = has_edge(g, src(e), dst(e), u)
has_edge(g::SimpleValueDiGraph, e::SimpleValueEdge) = has_edge(g, src(e), dst(e), val(e))

# TODO could probably be made faster by checking the shorter list
function get_edgeval(g::InOutValueDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    E_VAL = edgeval_type(g)
    (s in verts && d in verts) || return default # TODO may raise bounds error?
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index)
    end
    return nothing
end

# TODO rest methods for get_value
# TODO lots of duplicated code
function get_edgeval(g::OutValueDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    E_VAL = edgeval_type(g)
    (s in verts && d in verts) || return default # TODO may raise bounds error?
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index)
    end
    return nothing
end

# TODO could probably be made faster by checking the shorter list
function get_edgeval(g::InOutValueDiGraph{V, <: TupleOrNamedTuple}, s::Integer, d::Integer, key) where {V <: Integer}
    verts = vertices(g)
    E_VAL = edgeval_type(g)
    (s in verts && d in verts) || return default # TODO may raise bounds error?
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index, key)
    end
    return nothing
end

function get_edgeval(g::OutValueDiGraph{V, <: TupleOrNamedTuple}, s::Integer, d::Integer, key) where {V <: Integer}
    verts = vertices(g)
    E_VAL = edgeval_type(g)
    (s in verts && d in verts) || return default # TODO may raise bounds error?
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index, key)
    end
    return nothing
end

function get_edgeval(g::ConstValueDiGraph{V, E_VAL}, s::Integer, d::Integer) where {V, E_VAL}
    if has_edge(g, s, d)
        return default_edgeval(E_VAL)
    end
    return nothing
end

function get_edgeval(g::ConstValueDiGraph{V, E_VAL}, s::Integer, d::Integer, key) where {V, E_VAL}
    if has_edge(g, s, d)
        return default_edgeval(E_VAL)[key]
    end
    return nothing
end

get_edgeval!(g::SimpleValueDiGraph, e::SimpleEdge)   = get_edgeval(g, src(e), dst(e))
get_edgeval!(g::SimpleValueDiGraph, e::SimpleEdge, key)   = get_edgeval(g, src(e), dst(e), key)

function set_edgeval!(g::OutValueDiGraph, s::Integer, d::Integer, value)
    verts = vertices(g)
    E_VAL = edgeval_type(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
    set_value_for_index!(g.edgevals, s, index, value)

    return true
end

function set_edgeval!(g::OutValueDiGraph{V, <: TupleOrNamedTuple}, s::Integer, d::Integer, key, value) where {V}
     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
    set_value_for_index!(g.edgevals, s, index, key, value)

    return true
end

function set_edgeval!(g::InOutValueDiGraph, s::Integer, d::Integer, value)
    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_value_for_index!(g.edgevals, s, index, value)

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

    return true
end

# TODO create generated methods for NamedTuple (and also Tuple?)
function set_edgeval!(g::InOutValueDiGraph{V, <: TupleOrNamedTuple}, s::Integer, d::Integer, key, value) where {V}
    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
    set_value_for_index!(g.edgevals, s, index, key, value)

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

    return true
end

set_edgeval!(g::Union{InOutValueDiGraph, OutValueDiGraph}, e::SimpleEdge, u)   = set_value!(g, src(e), dst(e), u)
set_edgeval!(g::Union{InOutValueDiGraph, OutValueDiGraph}, e::SimpleValueEdge) = set_value!(g, src(e), dst(e), edgeval(e))

is_directed(::Type{<:SimpleValueDiGraph}) = true
is_directed(g::SimpleValueDiGraph) where {T, U} = true

outneighbors(g::SimpleValueDiGraph, v::Integer) = g.fadjlist[v]
inneighbors(g::SimpleValueDiGraph,  v::Integer) = g.badjlist[v]
# TODO neightbors, all_neighbors

# TODO implement these with iterators instead of generators
outedgevals(g::Union{InOutValueDiGraph{V, E_VAL, <: Adjlist},
                     OutValueDiGraph{V, E_VAL, <: Adjlist}}, v::Integer) where {V, E_VAL} =
    g.edgevals[v]
outedegevals(g::ConstValueDiGraph{V, E_VAL}, v::Integer) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL), length(g.fadjlist[v]))

inedgevals(g::InOutValueDiGraph{V, E_VAL , <: Adjlist}, v::Integer) where {V, E_VAL} = r.redgevals[v]
inedegevals(g::ConstValueDiGraph{V, E_VAL}, v::Integer) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL), length(g.fadjlist[v]))
#TODO all_edgevals

# TODO implement these with iterators instead of generators
outedgevals(g::Union{InOutValueDiGraph{V, E_VAL, <: Tuple},
                     OutValueDiGraph{V, E_VAL, <: Tuple}}, v::Integer) where {V, E_VAL} =
    ( Tuple( adjlist[v][i] for adjlist in g.edgevals ) for i in OneTo(length(g.fadjlist[v])) )
outedgevals(g::Union{InOutValueDiGraph{V, E_VAL, T},
                     OutValueDiGraph{V, E_VAL, T}}, v::Integer) where {V, E_VAL, T <: NamedTuple} =
    ( NamedTuple{Tuple(T.names)}( adjlist[v][i] for adjlist in g.edgevals ) for i in OneTo(length(g.fadjlist[v])) )

inedgevals(g::InOutValueDiGraph{V, E_VAL, <: Tuple}, v::Integer) where {V, E_VAL} =
    ( Tuple( adjlist[v][i] for adjlist in g.redgevals ) for i in OneTo(length(g.badjlist[v])) )
inedgevals(g::InOutValueDiGraph{V, E_VAL, T}, v::Integer) where {V, E_VAL, T <: NamedTuple} =
    ( NamedTuple{Tuple(T.names)}( adjlist[v][i] for adjlist in g.redgevals ) for i in OneTo(length(g.badjlist[v])) )

outedgevals(g::Union{InOutValueDiGraph{V, E_VAL, <: TupleOrNamedTuple},
                    OutValueDiGraph{V, E_VAL, <: TupleOrNamedTuple}},
           v::Integer, key) where {V, E_VAL} = g.edgevals[key][v]
outedegevals(g::ConstValueDiGraph{V, E_VAL}, v::Integer, key) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL)[key], length(g.fadjlist[v]))

inedgevals(g::InOutValueDiGraph{V, E_VAL, <: TupleOrNamedTuple},
            v::Integer,
            key) where {V, E_VAL} = g.redgevals[key][v]
inedegevals(g::ConstValueDiGraph{V, E_VAL}, v::Integer, key) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL)[key], length(g.badjlist[v]))

# TODO maybe add a sizehint kwarg
function add_vertex!(g::InOutValueDiGraph{V, E_VAL_C, E_VAL_C}) where {V, E_VAL, E_VAL_C}
    # TODO There are overflow checks in Julia Base, use these
    (nv(g) + one(V) <= nv(g)) && return false # overflow
    push!(g.fadjlist, V[])
    push!(g.badjlist, V[])
    if E_VAL_C <: TupleOrNamedTuple
        for (i, T) in enumerate(E_VAL.types)
            push!(g.edgevals[i], T[])
            push!(g.redgevals[i], T[])
        end
    else
        push!(g.edgevals, E_VAL[])
        push!(g.redgevals, E_VAL[])
    end
    return true
end

# TODO maybe add a sizehint kwarg
function add_vertex!(g::OutValueDiGraph{V, E_VAL, E_VAL_C}) where {V, E_VAL, E_VAL_C}
    # TODO There are overflow checks in Julia Base, use these
    (nv(g) + one(V) <= nv(g)) && return false # overflow
    push!(g.fadjlist, V[])
    push!(g.badjlist, V[])
    if E_VAL_C <: TupleOrNamedTuple
        for (i, T) in enumerate(E_VAL.types)
            push!(g.edgevals[i], T[])
        end
    else
        push!(g.edgevals, E_VAL[])
    end
    return true
end

# TODO maybe add a sizehint kwarg
function add_vertex!(g::ConstValueDiGraph{V}) where {V}
    # TODO There are overflow checks in Julia Base, use these
    (nv(g) + one(V) <= nv(g)) && return false # overflow
    push!(g.fadjlist, V[])
    push!(g.badjlist, V[])
    return true
end

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

function iterate(iter::SimpleValueEdgeIter{<:Union{InOutValueDiGraph, OutValueDiGraph}},
                 state=(one(eltype(iter.g)), 1) )
    g = iter.g
    fadjlist = g.fadjlist
    V = eltype(g)
    E_VAL = edgeval_type(g)
    n::V = nv(g)
    u::V, i = state

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

            u += 1
            i = 1
            continue
        end
        e = SimpleValueEdge(u, fadjlist[u][i], value_for_index(g.edgevals, E_VAL, u, i))
        return e, (u, i + 1)
    end

    return nothing
end

function iterate(iter::SimpleValueEdgeIter{<:ConstValueDiGraph},
                 state=(one(eltype(iter.g)), 1) )
    g = iter.g
    fadjlist = g.fadjlist
    V = eltype(g)
    n::V = nv(g)
    u::V, i = state

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

            u += 1
            i = 1
            continue
        end
        e = SimpleValueEdge(u, fadjlist[u][i])
        return e, (u, i + 1)
    end

    return nothing
end


text:
"""
    Workflow(actionname, files_to_create, [files_to_delete])

Create or modify a GitHub Action workflow. `actionname` is the name given to the action when it runs.
`files_to_create` should consist of one or more `filename => filecontents` pairs specifying the workflow file(s)
you want to create or update. The optional `files_to_delete` lists one or more workflow files you want to delete.
"""
struct Workflow
    name::String
    files_to_create::Dict{String, String}
    files_to_delete::Set{String}

    # Block the generic fallback
    Workflow(name::String, files_to_create::Dict{String, String}, files_to_delete::Set{String}) = new(name, files_to_create, files_to_delete)
end

Workflow(name::AbstractString, files_to_create, files_to_delete) = Workflow(String(name), normdict(files_to_create), normss(files_to_delete))
Workflow(name::AbstractString, files_to_create) = Workflow(String(name), normdict(files_to_create), Set{String}())

normdict(files_to_create::AbstractDict) = convert(Dict{String,String}, files_to_create)
normdict(files_to_create::Pair) = Dict{String,String}(files_to_create)
normdict(files_to_create::AbstractVector{<:Pair}) = Dict{String,String}(files_to_create)

normss(files_to_delete::AbstractSet) = convert(Set{String}, files_to_delete)
normss(files_to_delete) = Set{String}(files_to_delete)
normss(files_to_delete::AbstractString) = normss([files_to_delete])

ext: jl


text:
module Delsarte

using AbstractAlgebra, ClusteredLowRankSolver, BasesAndSamples, QuadGK

export delsarte, Nspherical_cap_packing

function delsarte(n, d,costheta, precision=512; all_free = false, kwargs...)
    #set up the polynomial field
    setprecision(precision)
    FF = RealField
    P, (u, ) = PolynomialRing(FF, ["u"])

    #compute both the gegenbauer polynomials and the sos basis
    gbasis = basis_gegenbauer(2d, n, u)
    samples = sample_points_chebyshev(2d,-1,costheta)
    sosbasis, samples = approximatefekete([u^k for k=0:2d], samples)

    #construct the constraint ∑_k a_k P^n_k(u) + SOS + (u+1)(cos(θ)-u)*SOS = - 1
    c = Dict()
    for k=0:2d
        c[Block(k)] = LowRankMatPol([gbasis[k+1]], [[P(1)]])
    end
    c[Block(:A)] = LowRankMatPol([P(1)], [sosbasis[1:d+1]])
    c[Block(:B)] = LowRankMatPol([(u+1)*(costheta-u)], [sosbasis[1:d]])
    constraint = Constraint(-P(1), c, Dict(), samples)

    #Construct the objective 1 + ∑_k a_k
    objective = Objective(1, Dict(Block(k) => hcat([FF(1)]) for k=0:2d), Dict())

    #Construct the SOS problem: minimize the objective s.t. the constraint
    sos = LowRankPolProblem(false, objective, [constraint])

    #Construct the SDP with or without using free variables for the a_k
    if all_free
        sdp = ClusteredLowRankSDP(sos,[k for k=0:2d])
    else
        sdp = ClusteredLowRankSDP(sos)
    end

    #Solve the SDP and return results
    return status, sol,time, errorcode = solvesdp(sdp; kwargs...)
end

function weights(theta, n)
    #approximate the integral up to a certain precision. For odd n, this can be done exactly (for even n, there is a term √(1-y^2) in the integral, so thats more difficult)
    # exact answer uses hypergeometric 2F1 functions. For n odd>=3, these terminate. Not sure about convergence at x=1.
    # exact answer of integral: [x 2F1(1/2, (3 - n)/2, 3/2, x)]_cos(theta)^1
    # Wolframalpha:  at 1 it equals sqrt(pi) Gamma(n/2-1/2) Gamma(n/2) = 1/the factor in front
    if n==3
        return 1//2*(1-cos(theta)) #odd n give an easy integral (a polynomial in the integration variable)
    elseif n==4
        return (theta-cos(theta)*abs(sin(theta)))/pi #converted to the type of theta
    elseif n==5
        return 1//2-3//4*cos(theta)+1//4 * cos(theta)^3
    end
    #for the others we just compute the integral numerically
    integral =
        sqrt(BigFloat(pi))^(-1) * gamma(BigFloat(n) / 2) / gamma(BigFloat(n - 1) / 2) *
        quadgk(y -> (1 - y^2)^((n - 3) / 2), cos(theta), BigFloat(1), rtol = 10^(-40))[1] # -30 goes relatively fast, -50 is slow
    return integral
end#

function Nspherical_cap_packing(n,d,thetas,N = length(thetas),precision=precision(BigFloat);all_free = false,solver_kwargs...)
    #Construct the polynomial ring and the gegenbauer polynomials
    setprecision(precision)
    FF = RealField
    w = weights.(thetas,n)

    P, (u, ) = PolynomialRing(FF, ["u"])
    basis = basis_gegenbauer(2d, n, u)

    constraints = []
    for i=1:N
        for j=1:i
            # sum_k a_k,i,j P_k^n(x) <= - sqrt(w(theta_i)*w(theta_j)) for x in [-1, cos(theta_i+theta_j)]
            # With sum of squares polynomials this equals:
            # sum_k a_k,i,j P_k^n(x) + SOS + (x+1)*(cos(theta_i+theta_j)-x)* SOS = - sqrt(w(theta_i)*w(theta_j))
            c = Dict()
            for k=0:2d
                if i != j
                    # We set both sides of the matrix
                    c[Block(k,i,j)] = LowRankMatPol([1//2*basis[k+1]], [[P(1)]])
                    c[Block(k,j,i)] = LowRankMatPol([1//2*basis[k+1]], [[P(1)]])
                else
                    c[Block(k,i,i)] = LowRankMatPol([basis[k+1]], [[P(1)]])
                end
            end
            c[Block((:SOS1,i,j))] = LowRankMatPol([P(1)], [[u^k for k=0:d]])
            c[Block((:SOS2,i,j))] = LowRankMatPol([(u+1)*(cos(thetas[i]+thetas[j])-u)], [[u^k for k=0:d-1]])

            samples = sample_points_chebyshev(2d,-1,cos(thetas[i]+thetas[j])) # TODO possibly take more samples
            push!(constraints, Constraint(-P(sqrt(w[i]*w[j])), c, Dict(), samples))
        end

        # The objective is M, but constrained by ∑_k a_kii ≦  M - w_i, i.e. ∑_k a_kii + pos.slack_i - M = -w_i
        obj = Dict(Block(k,i,i) => LowRankMatPol([P(1)],[[P(1)]]) for k=0:2d)
        #an easy mistake would be to use the same slack variable for every i by forgetting the i in the tuple (:SOS_obj, i).
        obj[Block((:SOS_obj,i))] = LowRankMatPol([P(1)],[[P(1)]])
        objective_constraint = Constraint(-P(w[i]),obj,Dict(:M=>P(-1)),[[0]])
        push!(constraints,objective_constraint)
    end

    #The objective is M
    objective = Objective(FF(0), Dict(), Dict(:M => FF(1)))

    #Construct the SOS - problem, minimizing the objective
    sos = LowRankPolProblem(false, objective, constraints)

    #construct and solve the SDP, with or without using free variables for the A_k
    if all_free
        sdp = ClusteredLowRankSDP(sos,[k for k=0:2d])
    else
        sdp = ClusteredLowRankSDP(sos)
    end
    status, sol, time, errorcode = solvesdp(sdp;solver_kwargs...)
end

end # end module


text:
using Hapi
using Documenter

makedocs(;
    modules=[Hapi],
    authors="Kevin Bonham <kevbonham@gmail.com>",
    repo="https://github.com/kescobo/Hapi.jl/blob/{commit}{path}#L{line}",
    sitename="Hapi.jl",
    strict=:doctest,
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        edit_link="main",
        canonical="https://kescobo.github.io/Hapi.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kescobo/Hapi.jl",
    push_preview=true,
    devbranch="main"
)


text:
module GP

import NLopt
import Base: show, size, length

include("cov.jl")
include("means.jl")
include("core.jl")
include("hyperparameters.jl")
include("search.jl")
include("ensemble.jl")
include("util.jl")

end


text:
"""
    OrderedBijector()

A bijector mapping ordered vectors in ℝᵈ to unordered vectors in ℝᵈ.

## See also
- [Stan's documentation](https://mc-stan.org/docs/2_27/reference-manual/ordered-vector.html)
  - Note that this transformation and its inverse are the _opposite_ of in this reference.
"""
struct OrderedBijector <: Bijector{1} end

"""
    ordered(d::Distribution)

Return a `Distribution` whose support are ordered vectors, i.e., vectors with increasingly ordered elements.
"""
ordered(d::ContinuousMultivariateDistribution) = Bijectors.transformed(d, OrderedBijector())

(b::OrderedBijector)(y::AbstractVecOrMat) = _transform_ordered(y)

function _transform_ordered(y::AbstractVector)
    x = similar(y)
    @assert !isempty(y)

    @inbounds x[1] = y[1]
    @inbounds for i = 2:length(x)
        x[i] = x[i - 1] + exp(y[i])
    end

    return x
end

function _transform_ordered(y::AbstractMatrix)
    x = similar(y)
    @assert !isempty(y)

    @inbounds for j = 1:size(x, 2), i = 1:size(x, 1)
        if i == 1
            x[i, j] = y[i, j]
        else
            x[i, j] = x[i - 1, j] + exp(y[i, j])
        end
    end

    return x
end

(ib::Inverse{<:OrderedBijector})(x::AbstractVecOrMat) = _transform_inverse_ordered(x)

function _transform_inverse_ordered(x::AbstractVector)
    y = similar(x)
    @assert !isempty(y)

    @inbounds y[1] = x[1]
    @inbounds for i = 2:length(y)
        y[i] = log(x[i] - x[i - 1])
    end

    return y
end

function _transform_inverse_ordered(x::AbstractMatrix)
    y = similar(x)
    @assert !isempty(y)

    @inbounds for j = 1:size(y, 2), i = 1:size(y, 1)
        if i == 1
            y[i, j] = x[i, j]
        else
            y[i, j] = log(x[i, j] - x[i - 1, j])
        end
    end

    return y
end

logabsdetjac(b::OrderedBijector, x::AbstractVector) = sum(@view(x[2:end]))
logabsdetjac(b::OrderedBijector, x::AbstractMatrix) = vec(sum(@view(x[2:end, :]); dims = 1))


text:
#
# Correctness Tests
#
using Test

using XMLDict
using BioMedQuery.Processes
using BioMedQuery.PubMed
using BioMedQuery.DBUtils
using BioMedQuery
using MySQL
using SQLite
using DataStreams

my_tests = [
            ("dbutils_sqlite.jl",   "       Testing: DBUtils SQLite"),
            ("dbutils_mysql.jl",    "       Testing: DBUtils MySQL"),
            ("pubmed.jl",           "       Testing: Eutils/PubMed"),
            ("export_citations.jl", "       Testing: Export Citations"),
            ("pubmed_parse.jl",     "       Testing: Entrez Parsing"),
            ("ct.jl",               "       Testing: CLINICAL TRIALS"),
            ("processes_mysql.jl",  "       Testing: Processes MySQL"),
            ("processes_sqlite.jl", "       Testing: Processes SQLite"),
            ("processes_df.jl",     "       Testing: Processes DataFrame")
            ]

println("Running tests:")

for (my_test, test_string) in my_tests
    println("-----------------------------------------")
    println("-----------------------------------------")
    println(test_string)
    println("-----------------------------------------")
    println("-----------------------------------------")

    include(my_test)
end


text:
function scatterplot(;x::AbstractVector = Int[],
                      y::AbstractVector = Int[],
                      group::AbstractVector = Int[],
                      pointShape::AbstractString = "circle",
                      pointSize::Real = 50)

    v = VegaVisualization()

    add_data!(v, x = x, y = y, group = group)

    #Get unique table name
    table = v.data[1].name

    default_scales!(v)
    default_axes!(v)

    #if non-null array for group, add legend

    if group != Int[]
        legend!(v)
    end

    #Old add_points! code
    res = VegaMark(_type = "symbol",
                from = VegaMarkFrom(data=table),
                properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                            y = VegaValueRef(scale = "y", field = "y"))
                                                )
                                )

    v.marks == nothing? v.marks = [res] : push!(v.marks, res)
    #End add_points! code

    v.marks[1].properties.enter.fill = VegaValueRef(scale = "group", field = "group")
    v.marks[1].properties.enter.size = VegaValueRef(value = pointSize)
    v.marks[1].properties.enter.shape = VegaValueRef(value = pointShape)

    return v

end

max_stars_repo_name: jpthor/vega


text:
#=
Window Discriminator Spike Sorting
=#

function canvas_release_win(widget::Ptr,param_tuple,user_data::Tuple{Gui_Handles,RHD2000})

    event = unsafe_load(param_tuple)
    han, rhd = user_data
    if event.button==1

        (x1,x2,y1,y2)=coordinate_transform(han,event)

        #If this is distributed, it is going to be *really* slow
        #because it will be sending over the entire DArray from whatever processor
        #Change this so it is just communicating with the cluster part; JNeuron does something similar
        if (han.clus==0)||(han.var2[han.spike,2]==0) #do nothing if zeroth cluster or window
        elseif (length(rhd.s[han.spike].c.win) < han.clus) #new cluster
            push!(rhd.s[han.spike].c.win,[SpikeSorting.mywin(x1,x2,y1,y2)])
            push!(rhd.s[han.spike].c.hits,0)
            han.total_clus[han.spike] += 1
            han.var2[han.spike,1]+=1
            han.spike_win=rhd.s[han.spike].c.win[end]
        elseif length(rhd.s[han.spike].c.win[han.clus]) < han.var2[han.spike,2] #new window
            push!(rhd.s[han.spike].c.win[han.clus],SpikeSorting.mywin(x1,x2,y1,y2))
            han.var2[han.spike,1]+=1
            han.spike_win=rhd.s[han.spike].c.win[han.clus]
        else #replace old window
            rhd.s[han.spike].c.win[han.clus][han.var2[han.spike,2]]=SpikeSorting.mywin(x1,x2,y1,y2)
            han.spike_win=rhd.s[han.spike].c.win[han.clus]
        end

        if ((han.clus>0)&(han.var2[han.spike,2]>0))&((han.buf_count>0)&(han.pause))

            window_cluster(han,han.clus)
            plot_new_color(getgc(han.c2),han,han.clus)
        end
    end

    nothing
end

#Delete clusters
function b1_cb_win(widgetptr::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    if (han.clus==0)||(han.clus>han.total_clus[han.spike]) #do nothing if zeroth cluster selected
    else
        deleteat!(rhd.s[han.spike].c.win,han.clus)
        pop!(rhd.s[han.spike].c.hits)
        han.total_clus[han.spike] -= 1
        han.clus = 0
        han.var2[han.spike,1] = 0
        han.var2[han.spike,2] = 0
        setproperty!(han.tb1,:label,string("Cluster: ",han.var1[han.spike,2]))
        setproperty!(han.tb2,:label,"Window: 0")
    end
    nothing
end

#Delete Window
function b2_cb_win(widgetptr::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    if (han.var2[han.spike,2]==0)||(han.var2[han.spike,2]>han.var2[han.spike,1]) #do nothing if zeroth window
    else
        deleteat!(rhd.s[han.spike].c.win[han.clus],han.var2[han.spike,2])
        han.var2[han.spike,1]-= 1
        han.var2[han.spike,2] = 0
        setproperty!(han.tb2,:label,"Window: 0")
    end

    nothing
end

#Display Windows
function b3_cb_win(widgetptr::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data
    ctx = getgc(han.c2)

    increment=div(500,han.wave_points)

    #loop over clusters
    for i=1:length(rhd.s[han.spike].c.win)
        #loop over windows
        for j=1:length(rhd.s[han.spike].c.win[i])

            x1=(rhd.s[han.spike].c.win[i][j].x1-1)*increment+1
            x2=(rhd.s[han.spike].c.win[i][j].x2-1)*increment+1
            y1=rhd.s[han.spike].c.win[i][j].y1
            y2=rhd.s[han.spike].c.win[i][j].y2
            move_to(ctx,x1,(y1-han.offset[han.spike])*han.scale[han.spike,1]+300)
            line_to(ctx,x2,(y2-han.offset[han.spike])*han.scale[han.spike,1]+300)
            set_line_width(ctx,5.0)
            select_color(ctx,i+1)
            stroke(ctx)
        end
    end

    reveal(han.c2)

    nothing
end

#Select Cluster
function b4_cb_win(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    #go to next cluster
    clus=han.clus+1

    if clus==han.total_clus[han.spike]+2
        han.clus=0
        han.var2[han.spike,1]=0
    elseif clus==han.total_clus[han.spike]+1
        #create new cluster
        han.clus=clus
        han.var2[han.spike,1]=0
    else
        han.clus=clus
        han.var2[han.spike,1]=length(rhd.s[han.spike].c.win)
    end

    #reset currently selected window to zero
    han.var2[han.spike,2]=0

    setproperty!(han.tb1,:label,string("Cluster: ",han.clus))
    setproperty!(han.tb2,:label,"Window: 0")

    nothing
end

#Select Window
function b5_cb_win(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    #go to next window
    win=han.var2[han.spike,2]+1

    if win==han.var2[han.spike,1]+2
        han.var2[han.spike,2]=0
    elseif win==han.var1[han.spike,1]+1
        #create new window
        han.var2[han.spike,2]=win
    else
        han.var2[han.spike,2]=win
    end

    setproperty!(han.tb2,:label,string("Window: ",han.var2[han.spike,2]))

    nothing
end

#Finds all spikes in buffer that are selected by windows
function window_cluster(han::Gui_Handles,clus::Int64)

    @inbounds for i=1:han.buf_ind

        hits=0
        for j=1:length(han.spike_win) #Loop over all windows
            a1=han.spike_win[j].x1
            a2=han.spike_win[j].x2
            b1=han.spike_win[j].y1
            b2=han.spike_win[j].y2
            for k=(han.spike_win[j].x1-1):(han.spike_win[j].x2+1)
                if SpikeSorting.intersect(a1,a2,k,k+1,b1,b2,han.spike_buf[k,i],han.spike_buf[k+1,i])
                    hits+=1
                    break
                end
            end
        end
        if hits==length(han.spike_win)
            han.buf_clus[i]=clus
        end
    end

    nothing
end


text:
using Base.Threads

export spawnloop

function spawnloop(f, n)
    fence[] = 0
    @threads for i = 1:n
        f(Val(i))
    end
end

const fence = Atomic{Int}(0)

function fun(::Val{N}) where N
    b = a = 0
    n = 0
    t = 0
    while a == b
        # t0 = time_ns()
        a = fence[]
        # t1 = time_ns()
        #t += Int(t1 - t0)
        n += 1
    end
    println(N, "(", threadid(), ") :", n, " ", t)
end

function fun(::Val{1})
    sleep(0.001)
    t0 = time_ns()
    fence[] = 1
    t1 = time_ns()
    println(1, "(", threadid(), ") time = ", Int(t1-t0))
end

function releaseto(a::Atomic{T}, v::T) where T
    a[] = v
    v
end
function releaseto(::Any, v::T) where T
    v
end

function waitfor(my::Atomic{T}, v::T) where T
    av = my[]
    while v > av
        av = my[]
    end
    av
end

function waitfor(delta::T, v::T) where T
    v + delta
end

function stepforward(my::Union{T,Atomic}, next::Union{T,Atomic}, irel::T, ineed::T) where T
    releaseto(next, irel)
    waitfor(my, ineed)
end

N = 10^6

function funsweep(n::Integer)
    A = [Atomic{Int}(0) for i = 1:n]
    @threads for t = 1:n
        println("starting thread $t")
        my = t == 1 ? 1000 : A[t]
        nx = t == n ? 0 : A[t+1]
        i = 1
        v = 1
        while i < N
            # println("a i = $i, tid = $t, v = $v")
            v = stepforward(my, nx, i, v)
            # println("b i = $i, tid = $t, v = $v")
            while i <= min(v, N)
                dojob(i, t)
                i += 1
            end
        end
        stepforward(0, nx, i, 0)
    end
end

const B = zeros(Int, N)
function dojob(i, t)
    B[i] |= 1<<(t-1)
    # println("i = $i , tid = $t")
end


text:
macro noop(ex)
    esc(ex)
end

const var"@inlinecall" = if VERSION < v"1.8"
    var"@noop"
else
    var"@inline"
end

@static if VERSION < v"1.8.0-DEV.410"
    using Base: @_inline_meta
else
    const var"@_inline_meta" = Base.var"@inline"
end


text:
isgpu(x) = false

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
        @info "DeconvOptim.jl: CUDA.jl is loaded, so include GPU functionality"

        gpu_or_cpu(x) = CUDA.CuArray
        isgpu(x::CUDA.CuArray) = true
        # prevent slow scalar indexing on GPU
        CUDA.allowscalar(false);

        # we need to fix some operations so that they are fast o GPUs
#        # Reference: https://discourse.julialang.org/t/cuarray-and-optim/14053
#        LinearAlgebra.norm1(x::CUDA.CuArray{T,N}) where {T,N} = sum(abs, x); # specializes the one-norm
#        LinearAlgebra.normInf(x::CUDA.CuArray{T,N}) where {T,N} = maximum(abs, x); # specializes the one-norm
#        Optim.maxdiff(x::CUDA.CuArray{T,N},y::CUDA.CuArray{T,N}) where {T,N} = maximum(abs.(x-y));

    end
end


text:
"""
    mutable struct RTEMatpowerGenerator <: AbstractNodeLabel

Mutable structure descendant of AbstractNodeLabel

# Fields
- `power_min::Complex128`
- `power_max::Complex128`
- `cost_degree::Int64`
- `cost_coeffs::Array{Float64}`
- `gen_status::Bool`: generator on/off

"""
mutable struct RTEMatpowerGenerator <: AbstractNodeLabel
  power_min::Complex #Smin = Pmin + i Qmin
  power_max::Complex #Smax = Pmax + i Qmax
  active_power_input::Float64
  gen_status::Bool
  flag::String
  voltage_input::Complex
end

"""
    create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: RTEMatpowerGenerator

Create power variables n for generator `elemid` at `bus` in `bus_vars` if `elem_formulation == :NewVar`\n
Return nothing

"""
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: RTEMatpowerGenerator
    bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)
    sign = element.flag
    bus_vars["lambda_$(sign)"] = Variable("lambda_$(sign)", Real)
    bus_vars["Volt"] = Variable(variable_name("VOLT", bus, "", scenario), Complex)
  return
end

"""
[cstrname, polynom, lb, ub] = Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: RTEMatpowerGenerator

Return the polynomial contribution in power of generator `elemid` at `bus`(name, value, lower bound, upper bound). Will be used to construct power balance constraints in polynomial problem.\n
If `elem_formulation == :NbMinVar`, return generator bounds `["UNIT", Polynomial() ,Smin, Smax]`\n
If `elem_formulation == :NewVar`, return -Sgen `["UNIT", -Sgen, 0, 0]`\n
Return no contribution `["", Polynomial(), 0, 0]` otherwise
"""
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: RTEMatpowerGenerator
  cstrname = "UNIT"
  return [cstrname, -bus_vars[elemid], 0, 0]

end

"""
    constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: RTEMatpowerGenerator

Return all the constraints defined by generator `elemid` at `bus`. Will be used to construct constraints in polynomial problem.\n
If `elem_formulation == :NewVar`, return generator bounds : "Genbounds" => Smin <= Sgen <= Smax\n
Return empty dictionary otherwise
"""
function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: RTEMatpowerGenerator
  #bounds on the reactive power
    Qmin = imag(element.power_min)
    Qmax = imag(element.power_max)
    reactive_power = 0.5*(-im)*(bus_vars[elemid] - conj(bus_vars[elemid]))
    Pmin = real(element.power_min)
    Pmax = real(element.power_max)
    Pinput = element.active_power_input
    if !element.gen_status
      Pmin = Pmax = 0
      Qmin = Qmax = 0
    end
    active_power = 0.5*(bus_vars[elemid]+conj(bus_vars[elemid]))
    sign = element.flag
    λ = bus_vars["lambda_$(sign)"]
    if sign == "minus"
      return SortedDict{String, Constraint}( "Reactivepower" => Qmin << reactive_power << Qmax,
      "Activepower" => -active_power + Pmin+2*(Pinput-Pmin)*λ == 0,
      "ctr_lambda_$(sign)" => 0 << λ << 0.5)
    elseif sign == "plus"
      return SortedDict{String, Constraint}( "Reactivepower" => Qmin << reactive_power << Qmax,
      "Activepower" => -active_power + 2*Pinput-Pmax+2*(Pmax-Pinput)*λ == 0,
      "ctr_lambda_$(sign)"=> 0.5 << λ << 1)
    else
      exit("Problem with flag")
    end
end

"""
    cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_vars::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: RTEMatpowerGenerator

Return the polynomial contribution in objective generator `elemid` at `bus`.
"""
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_var::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: RTEMatpowerGenerator
  #autre option : minimisation ecart tensions
  v = bus_elems_var["Volt"]
  v_input = element.voltage_input
  return abs2(v-v_input)
end


text:
include("ring.jl")
using .rings
using Plots

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
module ConcreteAbstractions

export
    @base,
    @extend

const _base_types = Dict{Symbol, Tuple}()

macro base(typeexpr::Expr)
    # dump(typeexpr)
    # name = :AbstractFoo # TODO: extract type name
    @assert typeexpr.head == :type
    mutable, nameblock, args = typeexpr.args

    # extract name and parameters
    name, params = if isa(nameblock, Symbol)
        nameblock, []
    elseif nameblock.head == :curly
        nameblock.args[1], nameblock.args[2:end]
    end

    # add this to our list of base types, then return the call to abstract
    _base_types[name] = params, args
    # @show name
    esc(quote
        abstract $name
    end)
end

macro extend(typeexpr::Expr)
    # dump(typeexpr)
    @assert typeexpr.head == :type
    mutable, nameblock, args = typeexpr.args
    # @show args

    # create a curly expression for the new type and grab the base name
    @assert isa(nameblock, Expr) && nameblock.head == :(<:)
    curly, basename = nameblock.args
    if isa(curly, Symbol)
        curly = :($curly{})
    end

    # check that we defined this base type before extracting the base params/args
    if !haskey(_base_types, basename)
        error("You must define a base type $basename")
    end
    params, baseargs = _base_types[basename]

    # add base parameters to def
    curly.args = vcat(curly.args[1], params, curly.args[2:end])
    nameblock.args[1] = curly

    # add base fields to def
    prepend!(args.args, baseargs.args)

    # return the new type expression
    esc(typeexpr)
end

end # module


text:
using BioMedQuery
using MySQL

@testset "MySQL BDUtils" begin

    dbname="pubmed_test"

    conn = DBUtils.init_mysql_database(MYSQL_HOST, MYSQL_USER, MYSQL_PASSWORD, dbname)
    PubMed.create_tables!(conn)

    println(conn)

    #check collection of tables
    tables_query = DBUtils.select_all_tables(conn)
    tables = ["basic","author_ref","mesh_desc","mesh_heading","mesh_qual","pub_type","abstract_structured","abstract_full","file_meta"]
    @test sort(tables) == sort(tables_query)

    #check minimum insert
    DBUtils.insert_row!(conn, "basic", Dict(:pmid => 1234, :title=>"Test Article", :pub_year=>nothing))

    #check duplicate error insert
    duplicate_id = DBUtils.insert_row!(conn, "basic", Dict(:pmid => 1234,
        :title=>"Test Article", :pub_year=>nothing))

    @test duplicate_id < 1
    sel = DBUtils.db_select(conn, ["pmid"], "basic", Dict(:title=>"Test Article"))

    @test length(sel[1]) == 1
    @test sel[1][1] == 1234

    # #clean up
    MySQL.execute!(conn, "DROP DATABASE IF EXISTS $dbname;")
end


text:
using BinaryBuilder

name = "EDFlib"
version = v"1.16.0"

sources = [
    ArchiveSource("https://www.teuniz.net/edflib/edflib_116.tar.gz",
                  "cc9f9cc63869fa5742a7dd7e1aa3ff69fedcd4547f2c56ada43d4a4bfa4c6a4e";
                  unpack_target="edflib")
]

script = raw"""
cd ${WORKSPACE}/srcdir/edflib
mkdir ${libdir}
${CC} edflib.c -shared -fPIC -D_LARGEFILE64_SOURCE -D_LARGEFILE_SOURCE -o ${libdir}/libedflib.${dlext}
# No separate license file, it just lives in the README and in the source files
install_license README.md
"""

platforms = supported_platforms()

products = [
    LibraryProduct("libedflib", :libedflib)
]

dependencies = [
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function ThermalGround(; name)
    @named hp = HeatPort()
    eqs = [hp.T ~ 0]
    ODESystem(eqs, t, systems=[hp], name=name)
end

function HeatCapacitor(; name, C=1.0)
    c_th = C

    @named hp = HeatPort()
    @parameters C
    @variables T(t) dt(t)

    D = Differential(t)
    eqs = [
        T ~ hp.T
        dt ~ D(T)
        D(T) ~ hp.Q_flow / C
        ]
    ODESystem(eqs, t, [T, dt], [C], systems=[hp], defaults=Dict(C => c_th), name=name)
end

function ThermalConductor(; name, G=1.0)
    g_th = G

    @named hp1 = HeatPort()
    @named hp2 = HeatPort()
    @parameters G
    @variables Q_flow(t) T(t)

    eqs = [
        T ~ hp1.T - hp2.T
        Q_flow ~ G*T
        Q_flow ~ hp1.Q_flow
        -Q_flow ~ hp2.Q_flow
    ]
    ODESystem(eqs, t, [Q_flow, T], [G], systems=[hp1, hp2], defaults=Dict(G => g_th), name=name)
end

function ThermalResistor(; name, R=1.0)
    r_th = R
    @named hp1 = HeatPort()
    @named hp2 = HeatPort()
    @parameters R
    @variables Q_flow(t) T(t)

    eqs = [
        T ~ R*Q_flow
        T ~ hp1.T - hp2.T
        hp1.Q_flow ~ Q_flow
        hp2.Q_flow ~ -Q_flow
    ]
    ODESystem(eqs, t, [Q_flow, T], [R], systems=[hp1, hp2], defaults=Dict(R => r_th), name=name)
end

function ConvectiveConductor(; name, G=1.0)
    g_c = G

    @named solidport = HeatPort()
    @named fluidport = HeatPort()
    @parameters G # Convective thermal conductance
    @variables Q_flow(t) dT(t)

    eqs = [
        dT ~ solidport.T - fluidport.T
        solidport.Q_flow ~ Q_flow
        fluidport.Q_flow ~ -Q_flow
        dT ~ G*Q_flow
    ]
    ODESystem(eqs, t, [Q_flow, dT], [G], systems=[solidport, fluidport], defaults=Dict(G => g_c), name=name)
end

function ConvectiveResistor(; name, R=1.0)
    r_c = R

    @named solidport = HeatPort()
    @named fluidport = HeatPort()
    @parameters R # Convective thermal resistance
    @variables Q_flow(t) dT(t)

    eqs = [
        dT ~ solidport.T - fluidport.T
        solidport.Q_flow ~ Q_flow
        fluidport.Q_flow ~ -Q_flow
        dT ~ R*Q_flow
    ]
    ODESystem(eqs, t, [Q_flow, dT], [R], systems=[solidport, fluidport], defaults=Dict(R => r_c), name=name)
end

function BodyRadiation(; name, G=1.0)
    g_r = G
    σ   = 5.6703744191844294e-8 # Stefan-Boltzmann constant

    @named hp1 = HeatPort()
    @named hp2 = HeatPort()
    @parameters G # Net radiation conductance between two surfaces
    @variables Q_flow(t)

    eqs = [
        Q_flow ~ G*σ*(hp1.T^4 - hp2.T^4)
    ]
    ODESystem(eqs, t, [Q_flow], [G], systems=[hp1, hp2], defaults=Dict(G => g_r), name=name)
end

function ThermalCollector(; name, N=1)
    hp = []
    for i in 1:N
        _hp = HeatPort(name=Symbol(:hp, i))
        push!(hp, _hp)
    end
    @named collector_port = HeatPort()

    eqs = [
        collector_port.Q_flow + sum(k -> k.Q_flow, hp) ~ 0
        collector_port.T ~ hp[1].T
    ]
    for i in 1:N-1
        push!(eqs, hp[i].T ~ hp[i+1].T)
    end
    ODESystem(eqs, t, [], [], systems=[hp..., collector_port], name=name)
end


text:
const PACKETS = Dict{PacketType, Type}(
CONNACK => Connack,
PUBLISH => Publish,
PUBACK => Puback,
PUBREC => Pubrec,
PUBREL => Pubrel,
PUBCOMP => Pubcomp,
SUBACK => Suback,
UNSUBACK => Unsuback,
PINGRESP => Pingresp)

function read_packet(s::IO)
    header = read(s, UInt8)
    len = read_len(s)
    # read variable header and payload into buffer
    buffer = PipeBuffer(read(s, len))
    packet_type = convert(PacketType,header & 0xF0)
    flags = header & 0x0F
    packet = read(buffer, flags, PACKETS[packet_type])
    return packet
end

function write_packet(s::IO, packet::Packet)
    buffer = PipeBuffer()
    write(buffer, packet)
    data = take!(buffer)
    write(s, packet.header)
    write_len(s, length(data))
    write(s, data)
end

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
# Autogenerated wrapper script for CHiPPS_BLIS_jll for x86_64-w64-mingw32-cxx11
export blis, libblis

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using ALPS_jll
using BiCePS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CHiPPS_BLIS")
JLLWrappers.@declare_library_product(libblis, "libBlis-0.dll")
JLLWrappers.@declare_executable_product(blis)
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, ALPS_jll, BiCePS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libblis,
        "bin\\libBlis-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        blis,
        "bin\\blis.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# -----------------------------------------
# Functions related to "traditional" decomposition
# -----------------------------------------
" Find the transformation matrix `H` to decompose the vector `s` into its parts and stacks them into `sbar`"
function find_decomposition_matrix!(ws)

  # allocate H and new decomposed cones
  n = COSMO.find_H_col_dimension(ws.p.C.sets, ws.ci.sp_arr)
  H_I = zeros(Int64, n)

  # find number of decomposed and total sets and allocate structure for new compositve convex set
  num_total, num_new_psd_cones = COSMO.num_cone_decomposition(ws)
  # decomposed_psd_cones = Array{COSMO.PsdCone}(undef, 0)
  C_new = Array{COSMO.AbstractConvexSet{Float64}}(undef, num_total)
  C_new[1] = COSMO.ZeroSet{Float64}(ws.ci.originalM)

  # loop over all convex sets and fill H and composite convex set accordingly
  row = 1
  entry = 1
  sp_ind = 1
  set_ind = 2
  for (kkk, C) in enumerate(ws.p.C.sets)
    set_ind, sp_ind, entry = COSMO.decompose!(H_I, C_new, set_ind, C, entry, row, ws.ci.sp_arr, sp_ind)
    row += C.dim
  end
  ws.p.C = COSMO.CompositeConvexSet(C_new)
  ws.ci.H = sparse(H_I, collect(1:n), ones(n))
end

"Determine the total number of sets `num_total` after decomposition and the number of new psd cones `num_new_psd_cones`."
function num_cone_decomposition(ws)
  num_sets = length(ws.p.C.sets)
  num_old_psd_cones = length(ws.ci.psd_cones_ind)
  num_new_psd_cones = 0
  for iii = 1:num_old_psd_cones
    sp = ws.ci.sp_arr[iii]
    num_new_psd_cones += COSMO.num_cliques(sp.sntree)
  end
  ws.ci.num_decom_psd_cones = ws.ci.num_psd_cones - ws.ci.num_decomposable + num_new_psd_cones
  # the total number is the number of original non-psd cones + number of decomposed psd cones + 1 zeroset for the augmented system
  num_total = num_sets - num_old_psd_cones + num_new_psd_cones + 1
  return num_total, num_new_psd_cones
end

function decompose!(H_I::Vector{Int64}, C_new, set_ind::Int64, C::COSMO.AbstractConvexSet, entry::Int64, row::Int64, sp_arr::Array{SparsityPattern}, sp_ind::Int64)
  for i = 1:C.dim
    H_I[entry] = row + i - 1
    entry += 1
  end
  C_new[set_ind] = C

  return set_ind + 1, sp_ind, entry
end

"For a PSD cone `C` of dimension `nBlk`, return the number of entries/rows in vectorized form."
get_blk_rows(nBlk::Int64, C::COSMO.PsdCone{T}) where {T} = Base.power_by_squaring(nBlk, 2)
get_blk_rows(nBlk::Int64, C::COSMO.PsdConeTriangle{T}) where {T} = div((nBlk + 1) * nBlk, 2)

function decompose!(H_I::Vector{Int64}, C_new, set_ind::Int64,  C::DecomposableCones{ <: Real}, entry::Int64, row_start::Int64, sp_arr::Array{SparsityPattern}, sp_ind::Int64)
  sparsity_pattern = sp_arr[sp_ind]
  sntree = sparsity_pattern.sntree
  # original matrix size
  original_size = C.sqrt_dim

  for iii = 1:num_cliques(sntree)

    c = COSMO.get_clique(sntree, iii)
    # the graph and tree algorithms determined the cliques vertices of an AMD-permuted matrix. Since the location of the data hasn't changed in reality, we have to map the clique vertices back
    c = map(v -> sparsity_pattern.ordering[v], c)
    sort!(c)
    entry = COSMO.add_subblock_map!(H_I, c, original_size, entry, row_start, C)

    # create and add new cone for subblock
    num_rows = get_blk_rows(get_nBlk(sntree, iii), C)
    C_new[set_ind] = typeof(C)(num_rows)
    set_ind += 1
  end
  return set_ind, sp_ind + 1, entry
end

# fills the corresponding entries of H for clique c
function add_subblock_map!(H_I::Vector{Int64}, clique_vertices::Array{Int64}, original_size::Int64, entry::Int64, row_start::Int64, ::PsdCone{<: Real})
  for vj in clique_vertices
    for vi in clique_vertices
      row = mat_to_vec_ind(vi, vj, original_size)
      H_I[entry] = row_start + row - 1
      entry += 1
    end
  end
  return entry::Int64
end

function add_subblock_map!(H_I::Vector{Int64}, clique_vertices::Array{Int64}, original_size::Int64, entry::Int64,  row_start::Int64, ::PsdConeTriangle{<: Real})
  for vj in clique_vertices
    for vi in clique_vertices
      if vi <= vj
        row = mat_to_svec_ind(vi, vj)
        H_I[entry] = row_start + row - 1
        entry += 1
      end
    end
  end
  return entry::Int64
end

function find_H_col_dimension(sets, sp_arr)
  sum_cols = 0
  sp_arr_ind = 1
  for C in sets
    dim, overlaps, sp_arr_ind = decomposed_dim(C, sp_arr, sp_arr_ind)
    sum_cols += dim
  end
  return sum_cols::Int64
end

function augment_system!(ws)
  _augment!(ws.p, ws.ci.H)
  nothing
end

function _augment!(problem, H::SparseMatrixCSC)
  mH, nH = size(H)
  m, n = size(problem.A)
  problem.P = blockdiag(problem.P, spzeros(nH, nH))
  problem.q = vec([problem.q; zeros(nH)])
  problem.A = [problem.A H; spzeros(nH, n) -sparse(1.0I, nH, nH)]
  problem.b = vec([problem.b; zeros(nH)])
  problem.model_size[1] = size(problem.A, 1)
  problem.model_size[2] = size(problem.A, 2)
  nothing
end

# -----------------------------------------
# Functions related to clique tree based transformation (compact decomposition)
# see: Kim: Exploiting sparsity in linear and nonlinear matrix inequalities via positive semidefinite matrix completion (2011), p.53
# -----------------------------------------

"""
    augment_clique_based!(ws::COSMO.Workspace{T})

If `settings.compact_decomposition = true`, decompose the problem based on the clique tree.
"""
function augment_clique_based!(ws::COSMO.Workspace{T}) where {T}
  A = ws.p.A
  b = ws.p.b
  q = ws.p.q
  P = ws.p.P
  cones = ws.p.C.sets
  sp_arr = ws.ci.sp_arr

  # determine number of final augmented matrix and number of overlapping entries
  mA, nA, num_overlapping_entries = find_A_dimension(ws.p.model_size[2], ws.p.C.sets, ws.ci.sp_arr)
  # find number of decomposed and total sets and allocate structure for new compositve convex set
  num_total, num_new_psd_cones = COSMO.num_cone_decomposition(ws)

  C_new = Array{COSMO.AbstractConvexSet{Float64}}(undef, num_total - 1)

  # allocate memory for Aa_I, Aa_J, Aa_V, b_I
  nz = nnz(A)
  Aa_I = zeros(Int64, nz + 2 * num_overlapping_entries)
  Aa_J = extra_columns(nz + 2 * num_overlapping_entries, nz + 1, ws.p.model_size[2] + 1)
  Aa_V =  alternating_sequence(nz + 2 * num_overlapping_entries, nz + 1)
  findnz!(Aa_I, Aa_J, Aa_V, A)
  bs = sparse(b)
  b_I = zeros(Int64, length(bs.nzval))
  b_V = bs.nzval

  row_ranges = COSMO.get_set_indices(cones) # the row ranges of each cone in the original problem
  row_ptr = 1 # a continuously updated pointer to the start of the first entry of a cone in A_I
  overlap_ptr = nz + 1 # a continuously updated pointer to the next free entry to enter the rows for the +1, -1 overlap entries

  sp_ind = 1 # every decomposable cone is linked to a sparsity_pattern with index sp_ind
  set_ind = 1 # the set index counter set_ind is used to create a map of the set number in the original problem to the set number in the decomposed problem

  # Main Loop: Loop over the cones of the original problem and add entries to the row vectors Aa_I and b_I
  for (k, C) in enumerate(cones)
    row_range = row_ranges[k]
    row_ptr, overlap_ptr, set_ind, sp_ind = COSMO.add_entries!(Aa_I, b_I, C_new, row_ptr, A, bs, row_range, overlap_ptr, set_ind, sp_ind, sp_arr, C, k, ws.ci.cone_map)
  end

  ws.p.A = allocate_sparse_matrix(Aa_I, Aa_J, Aa_V, mA, nA)
  ws.p.b = Vector(SparseArrays._sparsevector!(b_I, b_V, mA))
  ws.p.P = blockdiag(P, spzeros(num_overlapping_entries, num_overlapping_entries))
  ws.p.q = vec([q; zeros(num_overlapping_entries)])
  ws.p.model_size[1] = size(ws.p.A, 1)
  ws.p.model_size[2] = size(ws.p.A, 2)
  ws.p.C = COSMO.CompositeConvexSet(C_new)
  return nothing
end

"Given the row, column, and nzval vectors and dimensions, assemble the sparse matrix `Aa` of the decomposed problem in a slightly more memory efficient way."
function allocate_sparse_matrix(Aa_I::Array{Int64, 1}, Aa_J::Array{Int64, 1}, Aa_V::Array{Float64, 1}, mA::Int64, nA::Int64)
  csrrowptr = zeros(Int64, mA + 1)
  csrcolval = zeros(Int64, length(Aa_I))
  csrnzval = zeros(Float64, length(Aa_I))
  klasttouch = zeros(Int64, nA + 1)
  csccolptr = zeros(Int64, nA + 1)
  # sort_col_wise!(Aa_I, Aa_V, A.colptr, size(A, 2))
  #Aa = SparseMatrixCSC{Float64, Int64}(mA, nA, Aa_J, Aa_I, Aa_V)
  Aa = SparseArrays.sparse!(Aa_I, Aa_J, Aa_V, mA, nA, +, klasttouch, csrrowptr, csrcolval, csrnzval, csccolptr, Aa_I, Aa_V )
end

"Return the dimension of the problem after a clique tree based decomposition, given the sparsity patterns in `sp_arr`."
function find_A_dimension(n_original::Int64, sets, sp_arr)
  num_cols = n_original
  num_overlapping_entries = 0
  num_rows = 0
  sp_arr_ind = 1
  for C in sets
    dim, overlaps, sp_arr_ind = decomposed_dim(C, sp_arr, sp_arr_ind)
    num_rows += dim
    num_overlapping_entries += overlaps
  end
  return num_rows::Int64, (num_cols + num_overlapping_entries)::Int64, num_overlapping_entries
end

# This method handles all non-decomposable sets C. For these non-decomposable sets you just have to shift the row indices of all entries by an offset
function add_entries!(A_I::Array{Int64, 1}, b_I::Array{Int64, 1}, C_new::Array{COSMO.AbstractConvexSet{Float64}, 1}, row_ptr::Int64, A0::SparseMatrixCSC, b0::SparseVector, row_range::UnitRange{Int64}, overlap_ptr::Int64, set_ind::Int64, sp_ind::Int64,
  sp_arr::Array{SparsityPattern, 1}, C::AbstractConvexSet, k::Int64, cone_map::Dict{Int64, Int64})

  m, n = size(A0)

  # let's handle vector b first
  offset = row_ptr - row_range.start

  row_range_col = COSMO.get_rows(b0, row_range)
  if row_range_col != nothing
    for k in row_range_col
      b_I[k] = b0.nzind[k] + offset
    end
  end

  # also create new set once
  C_new[set_ind] = C
  cone_map[set_ind] = k

  for col = 1:n
    ptr = copy(row_ptr)
    # indices which store the rows in column for C in A0
    row_range_col = COSMO.get_rows(A0, col, row_range)
    for k in row_range_col
      A_I[k] = A0.rowval[k] + offset
    end
  end

  return row_ptr + C.dim, overlap_ptr, set_ind + 1, sp_ind
end

# This method handles decomposable cones, e.g. PsdConeTriangle. The row vectors A_I and b_I have to be edited in such a way that entries of one clique are appear continously
function add_entries!(A_I::Array{Int64, 1}, b_I::Array{Int64, 1}, C_new::Array{COSMO.AbstractConvexSet{Float64}, 1}, row_ptr::Int64, A0::SparseMatrixCSC, b0::SparseVector, row_range::UnitRange{Int64}, overlap_ptr::Int64, set_ind::Int64, sp_ind::Int64,
  sp_arr::Array{SparsityPattern, 1},  C::PsdConeTriangle{Float64}, k::Int64,  cone_map::Dict{Int64, Int64})

  sp = sp_arr[sp_ind] # The SparsityPattern correspondig to this cone C
  sntree = sp.sntree # The Supernodal Elimination Tree that stores information about the cliques
  ordering = sp.ordering # A reordering that was applied by the LDL routine
  N_v = length(ordering)
  m, n = size(A0) # Dimensions of the original problem

  # determine the row ranges for each of the subblocks
  clique_to_rows = COSMO.clique_rows_map(row_ptr, sntree, C)

  # loop over cliques in descending topological order
  for iii = num_cliques(sntree):-1:1

    # get supernodes and seperators and undo the reordering
    sep = map(v -> ordering[v], get_sep(sntree, iii))
    isa(sep, Array{Any, 1}) && (sep = Int64[])
    snd = map(v -> ordering[v], get_snd(sntree, iii))
    # compute sorted block indices (i, j, flag) for this clique with an information flag whether an entry (i, j) is an overlap
    block_indices = COSMO.get_block_indices(snd, sep, N_v)

    # If we encounter an overlap with a parent clique we have to be able to find the location of the overlapping entry
    # Therefore load and reorder the parent clique
    if iii == num_cliques(sntree)
      par_clique = Int64[]
      par_rows = 0:0
    else
      par_ind = COSMO.get_clique_par(sntree, iii)
      par_rows = clique_to_rows[par_ind]
      par_clique = map(v -> ordering[v], get_clique_by_ind(sntree, par_ind))
      sort!(par_clique)
    end

    # Loop over all the columns and shift the rows in A_I and b_I according to the clique strucutre
    for col = 1:n
      row_range_col = COSMO.get_rows(A0, col, row_range)
      row_range_b = col == 1 ? COSMO.get_rows(b0, row_range) : 0:0
      overlap_ptr = add_clique_entries!(A_I, b_I, A0.rowval, b0.nzind, block_indices, par_clique, par_rows, col, C.sqrt_dim, row_ptr, overlap_ptr, row_range, row_range_col, row_range_b)
    end

    # create and add new cone for subblock
    num_rows = get_blk_rows(get_nBlk(sntree, iii), C)
    cone_map[set_ind] = k
    C_new[set_ind] = typeof(C)(num_rows, sp_ind, iii)
    row_ptr += num_rows
    set_ind += 1
  end
  return row_ptr, overlap_ptr, set_ind, sp_ind + 1
end

" Loop over all entries (i, j) in the clique and either set the correct row in `A_I` and `b_I` if (i, j) is not an overlap or add an overlap column with (-1 and +1) in the correct positions."
function add_clique_entries!(A_I::Array{Int64, 1}, b_I::Array{Int64, 1}, A_rowval::Array{Int64}, b_nzind::Array{Int64, 1}, block_indices::Array{Tuple{Int64, Int64, Int64},1},  par_clique::Array{Int64, 1}, par_rows::UnitRange{Int64}, col::Int64,  C_sqrt_dim::Int64, row_ptr::Int64, overlap_ptr::Int64, row_range::UnitRange{Int64}, row_range_col::UnitRange{Int64}, row_range_b::UnitRange{Int64})
  counter = 0
  for block_idx in block_indices
    new_row_val = row_ptr + counter
    # a block index that corresponds to an overlap
    if block_idx[3] == 0
      if col == 1
        i = block_idx[1]
        j = block_idx[2]
        A_I[overlap_ptr] = new_row_val # this creates the +1 entry
        A_I[overlap_ptr + 1] = par_rows.start + COSMO.parent_block_indices(par_clique, i, j) - 1 # this creates the -1 entry
        overlap_ptr += 2
      end
    else
      # (i, j) of the clique
      i = block_idx[1]
      j = block_idx[2]
      # k = svec(i, j)
      k = COSMO.mat_to_svec_ind(i, j)
      modify_clique_rows!(A_I, k, A_rowval, C_sqrt_dim, new_row_val, row_range, row_range_col)
      col == 1 && modify_clique_rows!(b_I, k, b_nzind, C_sqrt_dim, new_row_val, row_range, row_range_b)
    end
  counter += 1
  end
  return overlap_ptr
end

" Given the nominal entry position `k = svec(i, j)` find and modify with `new_row_val` the actual location of that entry in the global row vector `rowval`."
function modify_clique_rows!(v::Array{Int64, 1}, k::Int64, rowval::Array{Int64, 1}, C_sqrt_dim::Int64, new_row_val::Int64, row_range::UnitRange{Int64}, row_range_col::UnitRange{Int64})
  row_0 = COSMO.get_row_index(k, rowval, C_sqrt_dim, row_range, row_range_col)
  # row_0 happens when (i, j) references an edge that was added by merging cliques, the corresponding value will be zero
  # and can be disregarded
  if row_0 != 0
    v[row_0] = new_row_val
  end
  return nothing
end

"Given the svec index `k` and an offset `row_range_col.start`, return the location of the (i, j)th entry in the row vector `rowval`."
function get_row_index(k::Int64, rowval::Array{Int64, 1}, sqrt_dim::Int64, row_range::UnitRange{Int64}, row_range_col::UnitRange{Int64})

  k_shift = row_range.start + k - 1

  # determine upper set boundary of where the row could be
  u = min(row_range_col.stop, row_range_col.start + k_shift - 1)

  # find index of first entry >= k, starting in the interval [l, u]
  # if no, entry is >= k, returns u + 1
  r = searchsortedfirst(rowval, k_shift, row_range_col.start, u, Base.Order.Forward)
  # if no r s.t. rowval[r] = k_shift was found that means that the (i, j)th entry represents an edded edge (zero) from clique merging
  if r > u || rowval[r] != k_shift
    return 0
  else
    return r
  end
end

# -----------------------------------------
# Utility Functions
# -----------------------------------------

" Find the index of k=svec(i, j) in the parent clique `par_clique`."
function parent_block_indices(par_clique::Array{Int64, 1}, i::Int64, j::Int64)
  ir = searchsortedfirst(par_clique, i)
  jr = searchsortedfirst(par_clique, j)
  return COSMO.mat_to_svec_ind(ir, jr)
end

"""
    get_block_indices(snd::Array{Int64}, sep::Array{Int64})

For a clique consisting of supernodes `snd` and seperators `sep`, compute all the indices (i, j) of the corresponding matrix block
in the format (i, j, flag) where flag is equal to 0 if entry (i, j) corresponds to an overlap of the clique and 1 otherwise.

`Nv` is the number of vertices in the graph that we are trying to decompose.
"""
function get_block_indices(snd::Array{Int64}, sep::Array{Int64}, Nv::Int64)
  N = length(sep) + length(snd)
  d = div(N * (N + 1), 2)

  block_indices = Array{Tuple{Int64, Int64, Int64}, 1}(undef, d)
  ind = 1

  for j in sep, i in sep
    if i <= j
      block_indices[ind] = (i, j, 0)
      ind += 1
    end
  end

  for j in snd, i in snd
    if i <= j
      block_indices[ind] = (i, j, 1)
      ind += 1
    end
  end

  for i in snd
    for j in sep
      block_indices[ind] = (min(i, j), max(i, j), 1)
      ind += 1
    end
  end

  sort!(block_indices, by = x -> x[2] * Nv + x[1] )
  return block_indices
end

"For a cone `C` determine the sum of vector-dimensions of its cliques after decomposition (`dim`) and the total number of overlaps between them (`overlaps`)."
decomposed_dim(C::AbstractConvexSet, sp_arr::Array{SparsityPattern}, sp_arr_ind::Int64) = (C.dim, 0, sp_arr_ind)
function decomposed_dim(C::DecomposableCones{ <: Real}, sp_arr::Array{SparsityPattern}, sp_arr_ind::Int64)
  sntree = sp_arr[sp_arr_ind].sntree
  dim, overlaps = get_decomposed_dim(sntree, C)
  return dim::Int64, overlaps::Int64, (sp_arr_ind + 1)::Int64
end

"Return the row ranges of each clique after the decomposition of `C` shifted by `row_start`."
function clique_rows_map(row_start::Int64, sntree::SuperNodeTree, C::DecomposableCones{<:Real})
  Nc = num_cliques(sntree)
  rows = Array{UnitRange{Int64}}(undef,  Nc)
  ind = zeros(Int64, Nc)
  for iii = Nc:-1:1
    num_rows = COSMO.get_blk_rows(COSMO.get_nBlk(sntree, iii), C)
    rows[iii] = row_start:row_start+num_rows-1
    ind[iii] = sntree.snd_post[iii]
    row_start += num_rows
  end
  return Dict(ind .=> rows)
end

function get_rows(b::SparseVector, row_range::UnitRange{Int64})
  rows = b.nzind
  if length(rows) > 0
    s = searchsortedfirst(rows, row_range.start)
    if rows[s] > row_range.stop || s == 0
        return nothing, 0:0
    else
      e = searchsortedlast(rows, row_range.stop)
      return s:e
    end
  else
    return nothing
  end

end

function get_rows(A::SparseMatrixCSC, col::Int64, row_range::UnitRange{Int64})
  colrange = A.colptr[col]:(A.colptr[col + 1]-1)

  # if the column has entries
  if colrange.start <= colrange.stop
    # create a view into the row values of column col
    rows = view(A.rowval, colrange)
    # find the rows within row_start:row_start+C.dim-1
    # s: index of first entry in rows >= row_start
    s = searchsortedfirst(rows, row_range.start)
    if rows[s] > row_range.stop || s == 0
      return 0:0
    else
      # e: index of last value in rows <= row_start + C.dim - 1
      e = searchsortedlast(rows, row_range.stop)
      return colrange[s]:colrange[e]
    end
  else
    return 0:0
  end
end

"Returns the appropriate amount of memory for `A.nzval`, including, starting from `n_start`, the (+1 -1) entries for the overlaps."
function alternating_sequence(total_length::Int64, n_start::Int64)
  v = ones(Float64, total_length)
  for i=n_start + 1:2:length(v)
    v[i] = -1
  end
  return v
end

"Returns the appropriate amount of memory for the columns of the augmented problem matrix `A`, including, starting from `n_start`, the columns for the (+1 -1) entries for the overlaps."
function extra_columns(total_length::Int64, n_start::Int64, start_val::Int64)
  v = zeros(Int64, total_length)
  for i = n_start:2:length(v)-1
    v[i] = start_val
    v[i + 1] = start_val
    start_val += 1
  end
  return v
end

function augmented_col_ptr(colptr::Array{Int64}, num_overlapping_entries::Int64)
  l = length(colptr)
  v = zeros(Int64, l + num_overlapping_entries)
  col = colptr[end] + 2
  for i = l+1:length(v)
    v[i] = col
    col += 2
  end
  return v
end

"Given a sparse matrix `S`, write the columns and non-zero values into the first `numnz` entries of `J` and `V`."
function findnz!(I::Vector{Ti}, J::Vector{Ti}, V::Vector{Tv}, S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    numnz = nnz(S)
    count = 1
    @inbounds for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1)
       # I[count] = S.rowval[k]
        J[count] = col
        V[count] = S.nzval[k]
        count += 1
    end
end

function add_nz_entries!(J::Vector{Ti}, V::Vector{Tv}, S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
  for i = 1:length(S.colptr)
    J[i] = S.colptr[i]
  end
  for i = 1:length(S.nzval)
    V[i] = S.nzval[i]
  end
  return nothing
end

"Given the row vector and colptr of a sparse matrix, sort the  entries within `rowvals` and `nzvals` withing the first `n0` columns."
function sort_col_wise!(rowval::Array{Int64, 1}, nzval::Array{Float64, 1}, colptr::Array{Int64, 1}, n0::Int64)
  for col = 1:n0
    r = colptr[col]:colptr[col + 1]-1
    row_view = view(rowval, r)
    nzval_view = view(nzval, r)
    p = sortperm(row_view)
    if !issorted(p)
      permute!(row_view, p)
      permute!(nzval_view, p)
    end
  end
  return nothing
end

function vec_to_mat_ind(ind::Int64, n::Int64)
  ind > n^2 && error("Index ind out of range.")
  ind == 1 && (return 1, 1)

  r = ind % n

  if r == 0
    j = Int(ind / n)
    i = n
  else
    j = Int(floor(ind / n) + 1)
    i = r
  end
  return i, j
end

"Convert the matrix indices `(i, j)` of `A ∈ R^{m x n}` into the corresponding linear index of `v = vec(A)`."
function mat_to_vec_ind(i::Int64, j::Int64, m::Int64)
  (i > m || i <= 0 || j <= 0) && throw(BoundsError("Indices outside matrix bounds."))
  return (j - 1) * m + i
end

"Convert the matrix indices `(i, j)` of `A ∈ R^{m x n}` into the corresponding linear index of `v = svec(A, ::UpperTriangle)`."
function mat_to_svec_ind(i::Int64, j::Int64)
  if i <= j
    return div((j - 1) * j, 2) + i
  else
    return div((i - 1) * i, 2) + j
  end
end

vectorized_ind(i::Int64, j::Int64, m::Int64, C::PsdCone{T}) where {T} = mat_to_vec_ind(i, j, m)
vectorized_ind(i::Int64, j::Int64, m::Int64, C::PsdConeTriangle{T}) where {T} = mat_to_svec_ind(i, j)

function svec_to_mat_ind(k::Int64)
  j = isqrt(2 * k)
  i = k - div((j - 1) * j, 2)
  return i, j
end


text:
@testset "Genera" begin
  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = PolynomialRing(K, "t")
  L, b = number_field(t^2 - gen(K) * t + 1)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  G = @inferred local_genera_hermitian(L, p, 4, 2, 4)
  @test length(G) == 15
  for i in 1:length(G)
    @test rank(G[i]) == 4
    @test (@inferred representative(G[i])) in G[i]
  end

  GG = G[1]
  u = @inferred uniformizer(GG)
  @assert parent(u) == K

  p = prime_decomposition(maximal_order(K), 17)[1][1]
  G = @inferred local_genera_hermitian(L, p, 5, 5, 5)
  @test length(G) == 7
  for i in 1:length(G)
    @test rank(G[i]) == 5
    @test (@inferred representative(G[i])) in G[i]
  end

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  L, b = number_field(t^2 - gen(K) * t + 1)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  l =  Vector{Tuple{Int, Int, Int, Int}}[[(0, 3, 1, 0), (4, 1, 1, 2)],
       [(0, 3, -1, 0), (4, 1, 1, 2)],
       [(0, 3, 1, 0), (4, 1, -1, 2)],
       [(0, 3, -1, 0), (4, 1, -1, 2)],
       [(0, 2, 1, 0), (2, 2, 1, 1)],
       [(0, 2, -1, 0), (2, 2, 1, 1)],
       [(0, 2, 1, 1), (2, 2, 1, 1)],
       [(0, 2, 1, 0), (2, 2, 1, 2)],
       [(0, 2, 1, 1), (2, 2, -1, 1)],
       [(0, 2, -1, 0), (2, 2, 1, 2)],
       [(0, 2, 1, 1), (2, 2, 1, 2)],
       [(0, 1, 1, 0), (1, 2, 1, 1), (2, 1, 1, 1)],
       [(0, 1, -1, 0), (1, 2, 1, 1), (2, 1, 1, 1)],
       [(1, 4, 1, 1)],
       [(1, 4, -1, 1)]]
  # The compiler does not like the following:
  # Gs = map(x -> genus(HermLat, L, p, x), l)
  Gs = Hecke.LocalGenusHerm{typeof(L), typeof(p)}[ genus(HermLat, L, p, x) for x in l ]
  myG = @inferred local_genera_hermitian(L, p, 4, 2, 4)
  @test length(Gs) == length(myG)
  @test all(x -> x in Gs, myG)
  @test all(x -> x in myG, Gs)

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  L, b = number_field(t^2 - gen(K) * t + 1)
  rlp = real_places(K)
  G = @inferred genera_hermitian(L, 3, Dict(rlp[1] => 1, rlp[2] => 1), 100 * maximal_order(L))
  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    M, = @inferred orthogonal_sum(representative(g1), representative(g2))
    @test M in (g1 + g2)
  end

  # Representatives

  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  E, b = NumberField(t^2 + 1, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  G = genus(HermLat, E, p, [(0, 3, -1, 0)])
  L = @inferred representative(G)
  @test length(@inferred Hecke.genus_representatives(L)) == 1

  K, a = NumberField(x^2 - 15)
  gens = [[a-1, 2*a-12, 0], [-471//70*a+881//14, 12*a-471, 1//70*a+39//14], [-7367*a+33891, 38904*a-212340, -194*a+1164], [2858191//5*a-1701731, -3700688*a+8438412, 103014//5*a-40352]]
  D = matrix(K, 3, 3, [38*a+150, 0, 0, 0, 2*a+8, 0, 0, 0, 302*a+1170])
  L = @inferred quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(@inferred genus_representatives(L)) == 1

  K, a = NumberField(x - 1)
  gens = [[-768, -1024, 0], [-768, -928, 160], [256, 320, -64]]
  D = matrix(K, 3, 3, [2, 0, 0, 0, 2, 0, 0, 0, 2])
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 2

  K, a = NumberField(x - 1)
  gens = [[-128, -1024, 0], [0, -2016, -32], [0, -512, 0]]
  D = matrix(K, 3, 3, [2, 0, 0, 0, 2, 0, 0, 0, 2])
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 2

  K, a = NumberField(x^2 - 2)
  gens = [[32*a-48, 32*a-64, 0], [0, 208*a+272, -24], [0, 512*a+768, -40*a-8]]
  D = matrix(K, 3, 3, [2, 0, 0, 0, 2, 0, 0, 0, 2])
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 2

  K, a = NumberField(x^3-x^2-2*x+1)
  gens = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 1//2, 0, 1//2]]
  D = matrix(K, 4, 4, [1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2])
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 2

  K, a = NumberField(x - 1)
  OK = maximal_order(K)
  p = prime_decomposition(OK, 2)[1][1]
  @test length(Hecke.local_genera_quadratic(K, p, rank = 2, det_val = 1)) == 8

  # Addition of genera

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = PolynomialRing(K, "t")
  L, b = number_field(t^2 - gen(K) * t + 1)

  p = prime_decomposition(maximal_order(K), 2)[1][1]
  G = local_genera_hermitian(L, p, 4, 2, 4)
  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    g3 = @inferred orthogonal_sum(g1, g2)
    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
  end

  p = prime_decomposition(maximal_order(K), 3)[1][1]
  G = local_genera_hermitian(L, p, 4, 2, 4)
  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    g3 = @inferred orthogonal_sum(g1, g2)
    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
  end

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [2, 0, 0, 0, -1, 0, 0, 0, -36]);
  gens = [[1, 0, 0], [1, 0, 0], [0, 1, 0], [0, 1, 0], [1//2, 1//2, 1//4], [1//2, 1//2, 1//4]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  p = prime_decomposition(maximal_order(K), 3)[1][1]
  fl, LL = @inferred Hecke.ismaximal_integral(L, p)
  @test !fl
  fl, _ = Hecke.ismaximal_integral(LL, p)
  @test fl

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [2, 0, 0, 0, -1, 0, 0, 0, -36]);
  V = quadratic_space(K, D)

  W = quadratic_space(K, diagonal_matrix([K(2), K(-36)]))

  @test @inferred Hecke.isrepresented_by(W, V)

  W = quadratic_space(K, diagonal_matrix([K(2), K(36)]))

  @test @inferred !Hecke.isrepresented_by(W, V)

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2-3
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(1)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(a), 0, 0, 0, (2 * a), 0, 0, 0, (a + 5)])
  V = hermitian_space(E, D)
  DD = matrix(E, 2, 2, [a + 3, 0, 0, (a + 2)])
  W = hermitian_space(E, DD)
  @test @inferred Hecke.isrepresented_by(W, V)
  @test @inferred Hecke.isrepresented_by(V, V)
  DD = E(a) * identity_matrix(E, 5)
  W = hermitian_space(E, DD)
  @test @inferred Hecke.isrepresented_by(V, W)

  # they are not isometric at some dyadic prime ideal:
  DD = identity_matrix(E, 4)
  W = hermitian_space(E, DD)
  @test !(@inferred Hecke.isrepresented_by(W, V))
end


text:
"""
Modified Morris-Lecar model from [Dhooge, Govaerts, Kuznetsov (2003)].

* [Dhooge, Govaerts, Kuznetsov (2003)].
  Numerical Continuation of Fold Bifurcations of Limit Cycles in MATCONT

[Dhooge, Govaerts, Kuznetsov (2003)]: https://doi.org/10.1007/3-540-44860-8_72
"""
module MorrisLecar

using DiffEqBase: ODEProblem
using Parameters: @with_kw, @unpack
using StaticArrays: SVector
using Setfield: @lens

using ...Bifurcations: BifurcationProblem
using ...Codim2: DiffEqCodim2Problem

@with_kw struct MorrisLecarParam{yType, zType}
    y::yType = 0.1104711840729923
    z::zType = 0.1
end

function f(u, p::MorrisLecarParam, t)
    @unpack y, z = p
    v, w = u
    m∞ = (1 + tanh((v + 0.01) / 0.15)) / 2
    w∞ = (1 + tanh((v - z) / 0.145)) / 2
    τ = cosh((v - 0.1) / 0.29)
    return SVector(
        y - 0.5(v + 0.5) - 2w * (v + 0.7) - m∞ * (v - 1),
        1.15(w∞ - w) * τ,
    )
end

function f(du, u, p::MorrisLecarParam, t)
    du .= f(SVector{2}(u), p, t)
    nothing
end

function make_prob(
        p = MorrisLecarParam();
        u0 = SVector(0.047222144637300754, 0.32564026295750237),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis = (@lens _.y),
        t_domain = (-0.05, 0.2),
        kwargs...)
    return BifurcationProblem(
        ode, param_axis, t_domain;
        kwargs...)
end

prob = make_prob()
ode = prob.p.de_prob
p = ode.p
u0 = ode.u0
param_axis = prob.p.param_axis

end  # module


text:
using Test
using NQCDynamics
using StochasticDiffEq
using LinearAlgebra
using DiffEqNoiseProcess
using DiffEqDevTools
using Random
using Logging
Random.seed!(100)

u0 = zeros(2,2)
v0 = ones(2,2)
γ = 1
f1_harmonic(v,u,p,t) = -u
f2_harmonic(v,u,p,t) = v
g(u,p,t) = diagm(ones(length(u)))

p = Simulation(Atoms([:H,:H]), NQCModels.Free(2); temperature=100)

ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,0.5),p)

sol1 = solve(prob1,DynamicsMethods.IntegrationAlgorithms.MDEF_BAOAB();dt=1/10,save_noise=true)

f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
g_iip(du,u,p,t) = du .= g(u,p,t)

prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,0.5),p;
                           noise=NoiseWrapper(sol1.W))

sol2 = solve(prob2,DynamicsMethods.IntegrationAlgorithms.MDEF_BAOAB();dt=1/10)
@test sol1[:] ≈ sol2[:]

disable_logging(Logging.Info)

dts = (1/2) .^ (8:-1:3)
res = analyticless_test_convergence(dts, prob1, DynamicsMethods.IntegrationAlgorithms.MDEF_BAOAB(), (1/2)^12;
    trajectories=Int(1e2), use_noise_grid=false)
@test res.𝒪est[:weak_final] ≈ 1 atol=0.5


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.BalanceLaws: vars_state, Prognostic, Auxiliary
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates: flattenednames
using ClimateMachine.Ocean.SplitExplicit01
using ClimateMachine.Ocean.OceanProblems
using ClimateMachine.GenericCallbacks
using ClimateMachine.VTK

using Test
using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

const ArrayType = ClimateMachine.array_type()

function main(BC)
    mpicomm = MPI.COMM_WORLD

    brickrange_2D = (xrange, yrange)
    topl_2D =
        BrickTopology(mpicomm, brickrange_2D, periodicity = (false, false))
    grid_2D = DiscontinuousSpectralElementGrid(
        topl_2D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    brickrange_3D = (xrange, yrange, zrange)
    topl_3D = StackedBrickTopology(
        mpicomm,
        brickrange_3D;
        periodicity = (false, false, false),
        boundary = ((1, 1), (1, 1), (2, 3)),
    )
    grid_3D = DiscontinuousSpectralElementGrid(
        topl_3D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    prob = OceanGyre{FT}(Lˣ, Lʸ, H; τₒ = τₒ, λʳ = λʳ, θᴱ = θᴱ, BC = BC)
    gravity::FT = grav(param_set)

    #- set model time-step:
    dt_fast = 240
    dt_slow = 5400

    nout = ceil(Int64, tout / dt_slow)
    dt_slow = tout / nout
    numImplSteps > 0 ? ivdc_dt = dt_slow / FT(numImplSteps) : ivdc_dt = dt_slow

    model = OceanModel{FT}(
        param_set,
        prob,
        cʰ = cʰ,
        add_fast_substeps = add_fast_substeps,
        numImplSteps = numImplSteps,
        ivdc_dt = ivdc_dt,
        κᶜ = FT(0.1),
    )
    # model = OceanModel{FT}(prob, cʰ = cʰ, fₒ = FT(0), β = FT(0) )
    # model = OceanModel{FT}(prob, cʰ = cʰ, νʰ = FT(1e3), νᶻ = FT(1e-3) )
    # model = OceanModel{FT}(prob, cʰ = cʰ, νʰ = FT(0), fₒ = FT(0), β = FT(0) )

    barotropicmodel = BarotropicModel(model)

    minΔx = min_node_distance(grid_3D, HorizontalDirection())
    minΔz = min_node_distance(grid_3D, VerticalDirection())
    #- 2 horiz directions
    gravity_max_dT = 1 / (2 * sqrt(gravity * H) / minΔx)
    # dt_fast = minimum([gravity_max_dT])

    #- 2 horiz directions + harmonic visc or diffusion: 2^2 factor in CFL:
    viscous_max_dT = 1 / (2 * model.νʰ / minΔx^2 + model.νᶻ / minΔz^2) / 4
    diffusive_max_dT = 1 / (2 * model.κʰ / minΔx^2 + model.κᶻ / minΔz^2) / 4
    # dt_slow = minimum([diffusive_max_dT, viscous_max_dT])

    @info @sprintf(
        """Update
           Gravity Max-dT = %.1f
           Timestep       = %.1f""",
        gravity_max_dT,
        dt_fast
    )

    @info @sprintf(
        """Update
       Viscous   Max-dT = %.1f
       Diffusive Max-dT = %.1f
       Timestep      = %.1f""",
        viscous_max_dT,
        diffusive_max_dT,
        dt_slow
    )

    barotropic_dg = DGModel(
        barotropicmodel,
        grid_2D,
        # CentralNumericalFluxFirstOrder(),
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q_2D = init_ode_state(barotropic_dg, FT(0); init_on_cpu = true)

    dg = OceanDGModel(
        model,
        grid_3D,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        modeldata = (barotropic_dg, Q_2D),
    )

    Q_3D = init_ode_state(dg, FT(0); init_on_cpu = true)

    lsrk_ocean = LSRK54CarpenterKennedy(dg, Q_3D, dt = dt_slow, t0 = 0)
    lsrk_barotropic =
        LSRK54CarpenterKennedy(barotropic_dg, Q_2D, dt = dt_fast, t0 = 0)

    odesolver = SplitExplicitLSRK2nSolver(lsrk_ocean, lsrk_barotropic)

    #-- Set up State Check call back for config state arrays, called every ntFreq time steps
    ntFreq = 1
    cbcs_dg = ClimateMachine.StateCheck.sccreate(
        [
            (Q_3D, "oce Q_3D"),
            (dg.state_auxiliary, "oce aux"),
            (Q_2D, "baro Q_2D"),
            (barotropic_dg.state_auxiliary, "baro aux"),
            # (dg.diffstate,"oce diff",),
            # (lsrk_ocean.dQ,"oce_dQ",),
            # (dg.modeldata.tendency_dg.state_auxiliary,"tend Int aux",),
            # (dg.modeldata.conti3d_Q,"conti3d_Q",),
            # (barotropic_dg.diffstate,"baro diff",),
            # (lsrk_barotropic.dQ,"baro_dQ",)
        ],
        ntFreq;
        prec = 12,
    )

    step = [0, 0]
    cbvector = make_callbacks(
        vtkpath,
        step,
        nout,
        mpicomm,
        odesolver,
        dg,
        model,
        Q_3D,
        barotropic_dg,
        barotropicmodel,
        Q_2D,
    )

    eng0 = norm(Q_3D)
    @info @sprintf """Starting
    norm(Q₀) = %.16e
    ArrayType = %s""" eng0 ArrayType

    cbv = (cbvector..., cbcs_dg)
    solve!(Q_3D, odesolver; timeend = timeend, callbacks = cbv)

    ## Enable the code block below to print table for use in reference value code
    ## reference value code sits in a file named $(@__FILE__)_refvals.jl. It is hand
    ## edited using code generated by block below when reference values are updated.
    regenRefVals = false
    if regenRefVals
        ## Print state statistics in format for use as reference values
        println(
            "# SC ========== Test number ",
            1,
            " reference values and precision match template. =======",
        )
        println("# SC ========== $(@__FILE__) test reference values ======================================")
        ClimateMachine.StateCheck.scprintref(cbcs_dg)
        println("# SC ====================================================================================")
    end

    ## Check results against reference if present
    checkRefVals = true
    if checkRefVals
        include("../refvals/simple_box_ivd_refvals.jl")
        refDat = (refVals[1], refPrecs[1])
        checkPass = ClimateMachine.StateCheck.scdocheck(cbcs_dg, refDat)
        checkPass ? checkRep = "Pass" : checkRep = "Fail"
        @test checkPass
        @info @sprintf("""Compare vs RefVals: %s""", checkRep)
    end

    return nothing
end

function make_callbacks(
    vtkpath,
    step,
    nout,
    mpicomm,
    odesolver,
    dg_slow,
    model_slow,
    Q_slow,
    dg_fast,
    model_fast,
    Q_fast,
)

    if isdir(vtkpath)
        rm(vtkpath, recursive = true)
    end
    mkpath(vtkpath)
    mkpath(vtkpath * "/slow")
    mkpath(vtkpath * "/fast")

    function do_output(span, step, model, dg, Q)
        outprefix = @sprintf(
            "%s/%s/mpirank%04d_step%04d",
            vtkpath,
            span,
            MPI.Comm_rank(mpicomm),
            step
        )
        @info "doing VTK output" outprefix
        statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
        auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
        writevtk(outprefix, Q, dg, statenames, dg.state_auxiliary, auxnames)

        mycomm = Q.mpicomm
        ## Generate the pvtu file for these vtk files
        if MPI.Comm_rank(mpicomm) == 0 && MPI.Comm_size(mpicomm) > 1
            ## name of the pvtu file
            pvtuprefix = @sprintf("%s/%s/step%04d", vtkpath, span, step)
            ## name of each of the ranks vtk files
            prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
                @sprintf("mpirank%04d_step%04d", i - 1, step)
            end
            writepvtu(
                pvtuprefix,
                prefixes,
                (statenames..., auxnames...),
                eltype(Q),
            )
            @info "Done writing VTK: $pvtuprefix"
        end

    end

    do_output("slow", step[1], model_slow, dg_slow, Q_slow)
    cbvtk_slow = GenericCallbacks.EveryXSimulationSteps(nout) do (init = false)
        do_output("slow", step[1], model_slow, dg_slow, Q_slow)
        step[1] += 1
        nothing
    end

    do_output("fast", step[2], model_fast, dg_fast, Q_fast)
    cbvtk_fast = GenericCallbacks.EveryXSimulationSteps(nout) do (init = false)
        do_output("fast", step[2], model_fast, dg_fast, Q_fast)
        step[2] += 1
        nothing
    end

    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q_slow)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end

    # return (cbvtk_slow, cbvtk_fast, cbinfo)
    return (cbinfo,)
end

#################
# RUN THE TESTS #
#################
FT = Float64
vtkpath = abspath(joinpath(ClimateMachine.Settings.output_dir, "vtk_split"))

const timeend = 5 * 24 * 3600 # s
const tout = 24 * 3600 # s

const N = 4
const Nˣ = 20
const Nʸ = 20
const Nᶻ = 20
const Lˣ = 4e6  # m
const Lʸ = 4e6  # m
const H = 1000  # m

xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

#const cʰ = sqrt(gravity * H)
const cʰ = 1  # typical of ocean internal-wave speed
const cᶻ = 0

#- inverse ratio of additional fast time steps (for weighted average)
#  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
# e.g., = 1 --> 100% more ; = 2 --> 50% more ; = 3 --> 33% more ...
add_fast_substeps = 2

#- number of Implicit vertical-diffusion sub-time-steps within one model full time-step
# default = 0 : disable implicit vertical diffusion
numImplSteps = 5

#const τₒ = 2e-1  # (Pa = N/m^2)
# since we are using old BC (with factor of 2), take only half:
const τₒ = 1e-1
const λʳ = 10 // 86400 # m/s
const θᴱ = 10    # deg.C

BC = (
    ClimateMachine.Ocean.SplitExplicit01.CoastlineNoSlip(),
    ClimateMachine.Ocean.SplitExplicit01.OceanFloorNoSlip(),
    ClimateMachine.Ocean.SplitExplicit01.OceanSurfaceStressForcing(),
)

@testset "$(@__FILE__)" begin
    main(BC)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA