text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    UniversalFallback

The `UniversalFallback` can be applied on a [`MathOptInterface.ModelLike`](@ref)
`model` to create the model `UniversalFallback(model)` supporting *any*
constraint and attribute. This allows to have a specialized implementation in
`model` for performance critical constraints and attributes while still
supporting other attributes with a small performance penalty. Note that `model`
is unaware of constraints and attributes stored by `UniversalFallback` so this
is not appropriate if `model` is an optimizer (for this reason,
[`MathOptInterface.optimize!`](@ref) has not been implemented). In that case,
optimizer bridges should be used instead.
"""
mutable struct UniversalFallback{MT} <: MOI.ModelLike
    model::MT
    objective::Union{MOI.AbstractScalarFunction,Nothing}
    # See https://github.com/jump-dev/JuMP.jl/issues/1152 and
    # https://github.com/jump-dev/JuMP.jl/issues/2238 for why we use an
    # `OrderedDict`
    single_variable_constraints::OrderedDict{Type,OrderedDict}
    constraints::OrderedDict{Tuple{Type,Type},VectorOfConstraints}
    con_to_name::Dict{MOI.ConstraintIndex,String}
    name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
    optattr::Dict{MOI.AbstractOptimizerAttribute,Any}
    modattr::Dict{MOI.AbstractModelAttribute,Any}
    varattr::Dict{MOI.AbstractVariableAttribute,Dict{VI,Any}}
    conattr::Dict{MOI.AbstractConstraintAttribute,Dict{MOI.ConstraintIndex,Any}}
    function UniversalFallback{MT}(model::MOI.ModelLike) where {MT}
        return new{typeof(model)}(
            model,
            nothing,
            OrderedDict{Type,OrderedDict}(),
            OrderedDict{Tuple{Type,Type},VectorOfConstraints}(),
            Dict{MOI.ConstraintIndex,String}(),
            nothing,
            Dict{MOI.AbstractOptimizerAttribute,Any}(),
            Dict{MOI.AbstractModelAttribute,Any}(),
            Dict{MOI.AbstractVariableAttribute,Dict{VI,Any}}(),
            Dict{
                MOI.AbstractConstraintAttribute,
                Dict{MOI.ConstraintIndex,Any},
            }(),
        )
    end
end

function UniversalFallback(model::MOI.ModelLike)
    return UniversalFallback{typeof(model)}(model)
end

function Base.show(io::IO, U::UniversalFallback)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    MOIU.print_with_acronym(io, summary(U))
    !(U.objective === nothing) && print(io, "\n$(indent)with objective")
    for (attr, name) in (
        (U.single_variable_constraints, "`VariableIndex` constraint"),
        (U.constraints, "constraint"),
        (U.optattr, "optimizer attribute"),
        (U.modattr, "model attribute"),
        (U.varattr, "variable attribute"),
        (U.conattr, "constraint attribute"),
    )
        n = length(attr)
        if n > 0
            print(io, "\n$(indent)with $n $name$(s(n))")
        end
    end
    print(io, "\n$(indent)fallback for ")
    return show(IOContext(io, :indent => get(io, :indent, 0) + 2), U.model)
end

function MOI.is_empty(uf::UniversalFallback)
    return MOI.is_empty(uf.model) &&
           uf.objective === nothing &&
           isempty(uf.single_variable_constraints) &&
           isempty(uf.constraints) &&
           isempty(uf.modattr) &&
           isempty(uf.varattr) &&
           isempty(uf.conattr)
end

function MOI.empty!(uf::UniversalFallback)
    MOI.empty!(uf.model)
    uf.objective = nothing
    empty!(uf.single_variable_constraints)
    empty!(uf.constraints)
    empty!(uf.con_to_name)
    uf.name_to_con = nothing
    empty!(uf.modattr)
    empty!(uf.varattr)
    empty!(uf.conattr)
    return
end

function pass_nonvariable_constraints(
    dest::UniversalFallback,
    src::MOI.ModelLike,
    idxmap::IndexMap,
    constraint_types,
)
    supported_types = eltype(constraint_types)[]
    unsupported_types = eltype(constraint_types)[]
    for (F, S) in constraint_types
        if MOI.supports_constraint(dest.model, F, S)
            push!(supported_types, (F, S))
        else
            push!(unsupported_types, (F, S))
        end
    end
    pass_nonvariable_constraints(dest.model, src, idxmap, supported_types)
    return pass_nonvariable_constraints_fallback(
        dest,
        src,
        idxmap,
        unsupported_types,
    )
end

function MOI.copy_to(dest::UniversalFallback, src::MOI.ModelLike; kwargs...)
    return MOIU.default_copy_to(dest, src; kwargs...)
end

function MOI.supports_incremental_interface(uf::UniversalFallback)
    return MOI.supports_incremental_interface(uf.model)
end

function final_touch(uf::UniversalFallback, index_map)
    return final_touch(uf.model, index_map)
end

"""
    throw_unsupported(uf::UniversalFallback; excluded_attributes = MOI.AnyAttribute[])

Throws [`MathOptInterface.UnsupportedAttribute`](@ref) if there are any model,
variable or constraint attribute such that 1) `is_copyable(attr)` returns
`true`, 2) the attribute was set to `uf` but not to `uf.model` and 3) the
attribute is not in `excluded_attributes`.

Suppose `Optimizer` supports only the constraints and attributes
that `OptimizerCache` supports in addition to
[`MathOptInterface.VariablePrimalStart`](@ref).
Then, this function can be used in the implementation of the following method:
```julia
function MOI.copy_to(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache},
)
    attr = MOI.VariablePrimalStart()
    MOI.Utilities.throw_unsupported(
        src,
        excluded_attributes = MOI.AnyAttribute[attr],
    )
    index_map = MOI.copy_to(dest, src.model)
    if attr in MOI.get(src, MOI.ListOfModelAttributesSet())
        for vi_src in MOI.get(src, MOI.ListOfVariableIndices())
            vi_dest = index_map[vi_src]
            value = MOI.get(src, attr, vi_src)
            MOI.set(dest, attr, vi_dest, value)
        end
    end
    return index_map
end
```
"""
function throw_unsupported(
    uf::UniversalFallback;
    excluded_attributes = MOI.AnyAttribute[],
)
    for attr in keys(uf.modattr)
        if !(attr in excluded_attributes)
            MOI.throw(MOI.UnsupportedAttribute(attr))
        end
    end
    for (attr, val) in uf.varattr
        if !isempty(val) && !(attr in excluded_attributes)
            MOI.throw(MOI.UnsupportedAttribute(attr))
        end
    end
    for (attr, val) in uf.conattr
        if !isempty(val) && !(attr in excluded_attributes)
            MOI.throw(MOI.UnsupportedAttribute(attr))
        end
    end
    for (S, val) in uf.single_variable_constraints
        if !isempty(val)
            throw(MOI.UnsupportedConstraint{MOI.VariableIndex,S}())
        end
    end
    for ((F, S), val) in uf.constraints
        if !MOI.is_empty(val)
            throw(MOI.UnsupportedConstraint{F,S}())
        end
    end
    return
end

# References
function MOI.is_valid(uf::UniversalFallback, vi::MOI.VariableIndex)
    return MOI.is_valid(uf.model, vi)
end

function MOI.is_valid(
    uf::UniversalFallback,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        return MOI.is_valid(uf.model, ci)
    end
    return haskey(uf.single_variable_constraints, S) &&
           haskey(uf.single_variable_constraints[S], ci)
end

function MOI.is_valid(
    uf::UniversalFallback,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if !MOI.supports_constraint(uf.model, F, S) &&
       !haskey(uf.constraints, (F, S))
        return false
    end
    return MOI.is_valid(constraints(uf, ci), ci)
end

function _delete(
    uf::UniversalFallback,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        MOI.delete(uf.model, ci)
    else
        MOI.is_valid(uf, ci) || throw(MOI.InvalidIndex(ci))
        delete!(uf.single_variable_constraints[S], ci)
    end
    return
end

function _delete(uf::UniversalFallback, ci::MOI.ConstraintIndex)
    MOI.delete(constraints(uf, ci), ci)
    return
end

function MOI.delete(
    uf::UniversalFallback,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    _delete(uf, ci)
    if !MOI.supports_constraint(uf.model, F, S)
        delete!(uf.con_to_name, ci)
        uf.name_to_con = nothing
    end
    for d in values(uf.conattr)
        delete!(d, ci)
    end
    return
end

function _remove_variable(
    uf::UniversalFallback,
    constraints::OrderedDict{<:MOI.ConstraintIndex{MOI.VariableIndex}},
    vi::MOI.VariableIndex,
)
    MOI.delete(uf, [ci for ci in keys(constraints) if ci.value == vi.value])
    return
end

function MOI.delete(uf::UniversalFallback, vi::MOI.VariableIndex)
    vis = [vi]
    for constraints in values(uf.constraints)
        _throw_if_cannot_delete(constraints, vis, vi)
    end
    MOI.delete(uf.model, vi)
    for d in values(uf.varattr)
        delete!(d, vi)
    end
    if uf.objective !== nothing
        uf.objective = remove_variable(uf.objective, vi)
    end
    for constraints in values(uf.single_variable_constraints)
        _remove_variable(uf, constraints, vi)
    end
    for constraints in values(uf.constraints)
        _deleted_constraints(constraints, vi) do ci
            delete!(uf.con_to_name, ci)
            uf.name_to_con = nothing
            for d in values(uf.conattr)
                delete!(d, ci)
            end
        end
    end
    return
end

function MOI.delete(uf::UniversalFallback, vis::Vector{MOI.VariableIndex})
    fast_in_vis = Set(vis)
    for constraints in values(uf.constraints)
        _throw_if_cannot_delete(constraints, vis, fast_in_vis)
    end
    MOI.delete(uf.model, vis)
    for d in values(uf.varattr)
        for vi in vis
            delete!(d, vi)
        end
    end
    if uf.objective !== nothing
        uf.objective = remove_variable(uf.objective, vis)
    end
    for constraints in values(uf.single_variable_constraints)
        for vi in vis
            _remove_variable(uf, constraints, vi)
        end
    end
    for constraints in values(uf.constraints)
        _deleted_constraints(constraints, vis) do ci
            delete!(uf.con_to_name, ci)
            uf.name_to_con = nothing
            for d in values(uf.conattr)
                delete!(d, ci)
            end
        end
    end
    return
end

# Attributes
function _get(uf::UniversalFallback, attr::MOI.AbstractOptimizerAttribute)
    return get(uf.optattr, attr, nothing)
end

function _get(uf::UniversalFallback, attr::MOI.AbstractModelAttribute)
    return get(uf.modattr, attr, nothing)
end

function _get(
    uf::UniversalFallback,
    attr::MOI.AbstractVariableAttribute,
    vi::MOI.VariableIndex,
)
    attribute_dict = get(uf.varattr, attr, nothing)
    if attribute_dict === nothing
        # It means the attribute is not set to any variable so in particular, it
        # is not set for `vi`
        return
    end
    return get(attribute_dict, vi, nothing)
end

function _get(
    uf::UniversalFallback,
    attr::MOI.AbstractConstraintAttribute,
    ci::MOI.ConstraintIndex,
)
    attribute_dict = get(uf.conattr, attr, nothing)
    if attribute_dict === nothing
        # It means the attribute is not set to any constraint so in particular,
        # it is not set for `ci`
        return
    end
    return get(attribute_dict, ci, nothing)
end

function _get(
    uf::UniversalFallback,
    attr::MOI.CanonicalConstraintFunction,
    ci::MOI.ConstraintIndex,
)
    return MOI.get_fallback(uf, attr, ci)
end

function MOI.get(
    uf::UniversalFallback,
    attr::Union{MOI.AbstractOptimizerAttribute,MOI.AbstractModelAttribute},
)
    if !MOI.is_copyable(attr) || MOI.supports(uf.model, attr)
        return MOI.get(uf.model, attr)
    end
    return _get(uf, attr)
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.AbstractConstraintAttribute,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if MOI.supports_constraint(uf.model, F, S) &&
       (!MOI.is_copyable(attr) || MOI.supports(uf.model, attr, typeof(ci)))
        return MOI.get(uf.model, attr, ci)
    end
    return _get(uf, attr, ci)
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.AbstractVariableAttribute,
    vi::MOI.VariableIndex,
)
    if !MOI.is_copyable(attr) || MOI.supports(uf.model, attr, MOI.VariableIndex)
        return MOI.get(uf.model, attr, vi)
    end
    return _get(uf, attr, vi)
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.NumberOfConstraints{MOI.VariableIndex,S},
)::Int64 where {S}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        return MOI.get(uf.model, attr)
    elseif !haskey(uf.single_variable_constraints, S)
        return 0
    else
        return length(uf.single_variable_constraints[S])
    end
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.NumberOfConstraints{F,S},
)::Int64 where {F,S}
    return MOI.get(constraints(uf, F, S), attr)
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
) where {S}
    F = MOI.VariableIndex
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, listattr)
    end
    indices = get(
        uf.single_variable_constraints,
        S,
        OrderedDict{MOI.ConstraintIndex{F,S},S}(),
    )
    return collect(keys(indices))
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfConstraintIndices{F,S},
) where {F,S}
    return MOI.get(constraints(uf, F, S), listattr)
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfConstraintTypesPresent,
)
    list = MOI.get(uf.model, listattr)
    for (S, constraints) in uf.single_variable_constraints
        if !isempty(constraints)
            push!(list, (MOI.VariableIndex, S))
        end
    end
    for (FS, constraints) in uf.constraints
        if !MOI.is_empty(constraints)
            push!(list, FS)
        end
    end
    return list
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfOptimizerAttributesSet,
)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.optattr)
        push!(list, attr)
    end
    return list
end

function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfModelAttributesSet)
    list = MOI.get(uf.model, listattr)
    if uf.objective !== nothing
        push!(list, MOI.ObjectiveFunction{typeof(uf.objective)}())
    end
    for attr in keys(uf.modattr)
        push!(list, attr)
    end
    return list
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfVariableAttributesSet,
)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.varattr)
        push!(list, attr)
    end
    return list
end

function MOI.get(
    uf::UniversalFallback,
    listattr::MOI.ListOfConstraintAttributesSet{F,S},
) where {F,S}
    list = MOI.get(uf.model, listattr)
    for (attr, dict) in uf.conattr
        if any(k -> k isa MOI.ConstraintIndex{F,S}, keys(dict))
            push!(list, attr)
        end
    end
    # ConstraintName isn't stored in conattr, but in the .con_to_name dictionary
    # instead. Only add it to `list` if it isn't already, and if a constraint
    # of the right type has a name set. This avoids, for example, returning
    # ConstraintName for VariableIndex constraints.
    if !(MOI.ConstraintName() in list)
        if any(k -> k isa MOI.ConstraintIndex{F,S}, keys(uf.con_to_name))
            push!(list, MOI.ConstraintName())
        end
    end
    return list
end

# Objective
function MOI.set(
    uf::UniversalFallback,
    attr::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    if sense == MOI.FEASIBILITY_SENSE
        uf.objective = nothing
    end
    MOI.set(uf.model, attr, sense)
    return
end

function MOI.get(uf::UniversalFallback, attr::MOI.ObjectiveFunctionType)
    if uf.objective === nothing
        return MOI.get(uf.model, attr)
    end
    return typeof(uf.objective)
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.ObjectiveFunction{F},
)::F where {F}
    if uf.objective === nothing
        return MOI.get(uf.model, attr)
    end
    return uf.objective
end

function MOI.set(
    uf::UniversalFallback,
    attr::MOI.ObjectiveFunction,
    func::MOI.AbstractScalarFunction,
)
    if MOI.supports(uf.model, attr)
        MOI.set(uf.model, attr, func)
        # Clear any fallback objective
        uf.objective = nothing
    else
        uf.objective = copy(func)
        # Clear any `model` objective
        sense = MOI.get(uf.model, MOI.ObjectiveSense())
        MOI.set(uf.model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
        MOI.set(uf.model, MOI.ObjectiveSense(), sense)
    end
    return
end

function MOI.modify(
    uf::UniversalFallback,
    obj::MOI.ObjectiveFunction,
    change::MOI.AbstractFunctionModification,
)
    if uf.objective === nothing
        MOI.modify(uf.model, obj, change)
    else
        uf.objective = modify_function(uf.objective, change)
    end
    return
end

# Name
# The names of constraints not supported by `uf.model` need to be handled
function MOI.set(
    uf::UniversalFallback,
    attr::MOI.ConstraintName,
    ci::MOI.ConstraintIndex{F,S},
    name::String,
) where {F,S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, attr, ci, name)
    else
        uf.con_to_name[ci] = name
        uf.name_to_con = nothing # Invalidate the name map.
    end
    return
end

function MOI.supports(
    ::UniversalFallback,
    ::MOI.ConstraintName,
    ::Type{<:MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.AbstractScalarSet}},
)
    return throw(MOI.VariableIndexConstraintNameError())
end

function MOI.set(
    ::UniversalFallback,
    ::MOI.ConstraintName,
    ::MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.AbstractScalarSet},
    ::String,
)
    return throw(MOI.VariableIndexConstraintNameError())
end

function MOI.get(
    uf::UniversalFallback,
    attr::MOI.ConstraintName,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, attr, ci)
    end
    return get(uf.con_to_name, ci, EMPTYSTRING)
end

function MOI.get(uf::UniversalFallback, ::Type{VI}, name::String)
    return MOI.get(uf.model, VI, name)
end

check_type_and_multiple_names(::Type, ::Nothing, ::Nothing, name) = nothing

function check_type_and_multiple_names(
    ::Type{T},
    value::T,
    ::Nothing,
    name,
) where {T}
    return value
end

check_type_and_multiple_names(::Type, ::Any, ::Nothing, name) = nothing

function check_type_and_multiple_names(
    ::Type{T},
    ::Nothing,
    value::T,
    name,
) where {T}
    return value
end

check_type_and_multiple_names(::Type, ::Nothing, ::Any, name) = nothing

function check_type_and_multiple_names(T::Type, ::Any, ::Any, name)
    return throw_multiple_name_error(T, name)
end

function MOI.get(
    uf::UniversalFallback,
    ::Type{MOI.ConstraintIndex{F,S}},
    name::String,
) where {F,S}
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end
    ci = if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, MOI.ConstraintIndex{F,S}, name)
    else
        # There is no `F`-in-`S` constraint in `b.model`, `ci` is only queried
        # to check for duplicate names.
        MOI.get(uf.model, MOI.ConstraintIndex, name)
    end
    ci_uf = get(uf.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci_uf, name)
    c = check_type_and_multiple_names(MOI.ConstraintIndex{F,S}, ci_uf, ci, name)
    return c
end

function MOI.get(
    uf::UniversalFallback,
    ::Type{MOI.ConstraintIndex},
    name::String,
)
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end
    ci_uf = get(uf.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci_uf, name)
    return check_type_and_multiple_names(
        MOI.ConstraintIndex,
        ci_uf,
        MOI.get(uf.model, MOI.ConstraintIndex, name),
        name,
    )
end

function _set(
    uf::UniversalFallback,
    attr::MOI.AbstractOptimizerAttribute,
    value,
)
    uf.optattr[attr] = value
    return
end

function _set(uf::UniversalFallback, attr::MOI.AbstractModelAttribute, value)
    uf.modattr[attr] = value
    return
end

function _set(
    uf::UniversalFallback,
    attr::MOI.AbstractVariableAttribute,
    vi::MOI.VariableIndex,
    value,
)
    if !haskey(uf.varattr, attr)
        uf.varattr[attr] = Dict{MOI.VariableIndex,Any}()
    end
    uf.varattr[attr][vi] = value
    return
end

function _set(
    uf::UniversalFallback,
    attr::MOI.AbstractConstraintAttribute,
    ci::MOI.ConstraintIndex,
    value,
)
    if !haskey(uf.conattr, attr)
        uf.conattr[attr] = Dict{MOI.ConstraintIndex,Any}()
    end
    uf.conattr[attr][ci] = value
    return
end

function MOI.supports(
    ::UniversalFallback,
    ::Union{MOI.AbstractModelAttribute,MOI.AbstractOptimizerAttribute},
)
    return true
end

function MOI.set(
    uf::UniversalFallback,
    attr::Union{MOI.AbstractOptimizerAttribute,MOI.AbstractModelAttribute},
    value,
)
    if MOI.supports(uf.model, attr)
        MOI.set(uf.model, attr, value)
    else
        _set(uf, attr, value)
    end
    return
end

function MOI.supports(
    ::UniversalFallback,
    ::Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute},
    ::Type{<:MOI.Index},
)
    return true
end

function MOI.set(
    uf::UniversalFallback,
    attr::MOI.AbstractVariableAttribute,
    vi::MOI.VariableIndex,
    value,
)
    if MOI.supports(uf.model, attr, MOI.VariableIndex)
        MOI.set(uf.model, attr, vi, value)
    else
        _set(uf, attr, vi, value)
    end
    return
end

function MOI.set(
    uf::UniversalFallback,
    attr::MOI.AbstractConstraintAttribute,
    ci::MOI.ConstraintIndex{F,S},
    value,
) where {F,S}
    if MOI.supports_constraint(uf.model, F, S) &&
       MOI.supports(uf.model, attr, MOI.ConstraintIndex{F,S})
        MOI.set(uf.model, attr, ci, value)
    else
        _set(uf, attr, ci, value)
    end
    return
end

# Constraints
function MOI.supports_constraint(
    ::UniversalFallback,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return true
end

function constraints(
    uf::UniversalFallback,
    ::Type{F},
    ::Type{S},
    getter::Function = get,
) where {F,S}
    if MOI.supports_constraint(uf.model, F, S)
        return uf.model
    end
    return getter(uf.constraints, (F, S)) do
        return VectorOfConstraints{F,S}()
    end::VectorOfConstraints{F,S}
end

function constraints(
    uf::UniversalFallback,
    ::MOI.ConstraintIndex{F,S},
) where {F,S}
    return constraints(uf, F, S)
end

function MOI.add_constraint(
    uf::UniversalFallback,
    func::MOI.VariableIndex,
    set::S,
) where {S<:MOI.AbstractScalarSet}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        return MOI.add_constraint(uf.model, func, set)
    end
    constraints = get!(
        uf.single_variable_constraints,
        S,
    ) do
        return OrderedDict{
            MOI.ConstraintIndex{MOI.VariableIndex,S},
            S,
        }()
    end::OrderedDict{MOI.ConstraintIndex{MOI.VariableIndex,S},S}
    ci = MOI.ConstraintIndex{MOI.VariableIndex,S}(func.value)
    constraints[ci] = set
    return ci
end

function MOI.add_constraint(
    uf::UniversalFallback,
    func::MOI.AbstractFunction,
    set::MOI.AbstractSet,
)
    return MOI.add_constraint(
        constraints(uf, typeof(func), typeof(set), get!),
        func,
        set,
    )
end

function MOI.modify(
    uf::UniversalFallback,
    ci::MOI.ConstraintIndex,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(constraints(uf, ci), ci, change)
    return
end

function MOI.get(
    uf::UniversalFallback,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
)
    return MOI.get(constraints(uf, ci), attr, ci)
end

function MOI.set(
    uf::UniversalFallback,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
    func_or_set,
)
    return MOI.set(constraints(uf, ci), attr, ci, func_or_set)
end

function MOI.get(
    uf::UniversalFallback,
    ::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{MOI.VariableIndex},
)
    MOI.throw_if_not_valid(uf, ci)
    return MOI.VariableIndex(ci.value)
end

function MOI.get(
    uf::UniversalFallback,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        return MOI.get(uf.model, MOI.ConstraintSet(), ci)
    end
    MOI.throw_if_not_valid(uf, ci)
    return uf.single_variable_constraints[S][ci]
end

function MOI.set(
    ::UniversalFallback,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex},
    ::MOI.VariableIndex,
)
    return throw(MOI.SettingVariableIndexNotAllowed())
end

function MOI.set(
    uf::UniversalFallback,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
    set::S,
) where {S}
    if MOI.supports_constraint(uf.model, MOI.VariableIndex, S)
        MOI.set(uf.model, MOI.ConstraintSet(), ci, set)
    else
        MOI.throw_if_not_valid(uf, ci)
        uf.single_variable_constraints[S][ci] = set
    end
    return
end

# Variables

MOI.add_variable(uf::UniversalFallback) = MOI.add_variable(uf.model)

MOI.add_variables(uf::UniversalFallback, n) = MOI.add_variables(uf.model, n)

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Calc
using Test

@testset "Calc.jl" begin
    @test add(1.0, 2.0) == 3.0
    @test subtract(1.0, 2.0) == -1.0
    @test multiply(1.0, 2.0) == 2.0
    @test divide(1.0, 2.0) == 0.5
    @test_throws DivideError() divide(1.0, 0.0)
end

ext: jl


text:
macro print_and_throw(args...)
    quote
        @cuprintln "ERROR: " $(args...) "."
        throw(nothing)
    end
end

# math.jl
@device_override @noinline Base.Math.throw_complex_domainerror(f::Symbol, x) =
    @print_and_throw "This operation requires a complex input to return a complex result"
@device_override @noinline Base.Math.throw_exp_domainerror(f::Symbol, x) =
    @print_and_throw "Exponentiation yielding a complex result requires a complex argument"

# intfuncs.jl
@device_override @noinline Base.throw_domerr_powbysq(::Any, p) =
    @print_and_throw "Cannot raise an integer to a negative power"
@device_override @noinline Base.throw_domerr_powbysq(::Integer, p) =
    @print_and_throw "Cannot raise an integer to a negative power"
@device_override @noinline Base.throw_domerr_powbysq(::AbstractMatrix, p) =
    @print_and_throw "Cannot raise an integer to a negative power"

# checked.jl
@device_override @noinline Base.Checked.throw_overflowerr_binaryop(op, x, y) =
    @print_and_throw "Binary operation overflowed"
@device_override @noinline Base.Checked.throw_overflowerr_negation(op, x, y) =
    @print_and_throw "Negation overflowed"

# boot.jl
@device_override @noinline Core.throw_inexacterror(f::Symbol, ::Type{T}, val) where {T} =
    @print_and_throw "Inexact conversion"

# abstractarray.jl
@device_override @noinline Base.throw_boundserror(A, I) =
    @print_and_throw "Out-of-bounds array access"

# trig.jl
@device_override @noinline Base.Math.sincos_domain_error(x) =
    @print_and_throw "sincos(x) is only defined for finite x."


text:
# ------------------------------------------------------------------------------------------
# # 多重派发
# 在这个 notebook 中我们将探索 Julia 的重要特性之一：**多重派发**。
#
# 多重派发让软件更**通用（generic）、更快**！
#
# ## 目录
#
# - [从熟悉的开始讲](#从熟悉的开始讲)
# - [指定传入变量的类型](#指定传入变量的类型)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 从熟悉的开始讲
#
# 从我们已经接触过的知识来理解多重派发。
#
# 我们可以声明一个 Julia 函数而不用给出传入变量的任何类型信息：
# ------------------------------------------------------------------------------------------

f(x) = x^2

# ------------------------------------------------------------------------------------------
# Julia 自会判断哪些传入变量是合理的，哪些又是不合理的：
# ------------------------------------------------------------------------------------------

f(10)

f([1, 2, 3])

# ------------------------------------------------------------------------------------------
# ## 指定传入变量的类型
#
# 但是，我们*也可以*显式地指定传入变量的类型。
#
# 比如，我们写一个只接受字符串类型传入变量的函数 `foo`。
# ------------------------------------------------------------------------------------------

foo(x::String, y::String) = println("My inputs x and y are both strings!")

# ------------------------------------------------------------------------------------------
# 可以看到，要限制 `x` 和 `y` 的类型为 `String`，只要在形参名后面加上两个冒号和关键字 `String` 就行了。
#
# 现在 `foo` 只有在传入变量为 `String` 类型才正常工作。
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# 要让 `foo` 函数接受整型（`Int`）的传入变量，只要在声明 `foo` 函数时在形参后面加上 `::Int`。
# ------------------------------------------------------------------------------------------

foo(x::Int, y::Int) = println("My inputs x and y are both integers!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# 现在函数 `foo` 能接受整型变量了！但是注意，函数 `foo` 仍然能接受字符串变量！
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

# ------------------------------------------------------------------------------------------
# 这就涉及到多重派发的核心了。当我们声明
#
# ```julia
# foo(x::Int, y::Int) = println("My inputs x and y are both integers!")
# ```
# 我们并没有重写或覆盖
# ```julia
# foo(y::String, y::String)
# ```
# 而是给叫做 `foo` 的**通用函数（generic function）**添加了一个**新方法（method）**。
#
# 一个**通用函数** 是一个表示特定操作的抽象概念。
# 举个例子，通用函数 `+` 代表了加法这个概念
#
# 一个**方法** 是一个**通用函数** 接受*特定参数类型* 的一个具体实现。
# 举个例子，`+` 有不同的方法以接受浮点数、整型、矩阵等等。
#
# 我们可以通过函数 `methods` 来查看函数 `foo` 有多少种方法。
# ------------------------------------------------------------------------------------------

methods(foo)

methods(+)

# ------------------------------------------------------------------------------------------
# 所以，我们现在可以在调用 `foo` 的时候传入整型或字符串了。
# 当你调用 `foo` 时传入特定组合的参数时，Julia 将根据传入参数的类型派发对应的方法。*这* 就是多重派发。
#
# 多重派发让我们的代码更加通用更加快。
# 代码更加通用更灵活是因为我们是在写代码描述抽象操作例如加法和乘法，而不是描述具体的实现。
# 同时，因为 Julia 可以为相关类型调用高效的方法，代码可以运行得很快。
#
# 想查看我们调用一个通用函数时派发的是哪个方法，我们可以用宏 `@which`：
# ------------------------------------------------------------------------------------------

@which foo(3, 4)

@which 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# 考虑到一个特别写给浮点数的方法派发给了 `3.0 + 3.0`，生成的 LLVM 代码极其简练：
# ------------------------------------------------------------------------------------------

@code_llvm 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# 要注意的是 Julia 在我们写一般的函数定义时也很快，因为具体的、定制的方法最终是在后台调用的。
#
# 举个例子，我们声明加法函数 `myadd` 时没有类型注释——
# ------------------------------------------------------------------------------------------

myadd(x, y) = x + y

# ------------------------------------------------------------------------------------------
# 虽然我们没有限制 `x` 和 `y` 的类型，但是 `myadd(3.0, 3.0)` 生成的 LLVM 代码和 `3.0 + 3.0` 差不多
# ------------------------------------------------------------------------------------------

@code_llvm myadd(3.0, 3.0)

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
using ReinforcementLearningEnvironmentClassicControl, ReinforcementLearning

env = CartPole()
rlsetup = RLSetup(ActorCriticPolicyGradient(ns = 4, na = 2, α = .02,
                                            nsteps = 25),
                  env, ConstantNumberSteps(400),
                  callbacks = [Visualize(wait = .02)])
@info("Before learning.")
run!(rlsetup)
rlsetup.callbacks = [EvaluationPerEpisode(TotalReward())]
rlsetup.stoppingcriterion = ConstantNumberSteps(10^5)
@time learn!(rlsetup)
getvalue(rlsetup.callbacks[1])
@info("After learning.")
rlsetup.callbacks = [Visualize(wait = .02)]
rlsetup.stoppingcriterion = ConstantNumberSteps(400)
run!(rlsetup)


text:
# Autogenerated wrapper script for Glibc_jll for armv6l-linux-gnueabihf
JLLWrappers.@generate_wrapper_header("Glibc")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using DFTK

include("testcases.jl")

# TODO Once we have converged SCF densities in a file it would be better to instead / also
#      test the energies of these densities and compare them directly to the reference
#      energies obtained in the data files

@testset "Evaluate energies of guess density" begin
    Ecut = 15
    n_bands = 8
    fft_size = [27, 27, 27]
    kgrid  = (1, 2, 3)
    kshift = (0, 1/2, 0)

    model = model_DFT(silicon.lattice, silicon.atoms, silicon.positions,
                      [:lda_x, :lda_c_vwn]; symmetries=false)
    basis = PlaneWaveBasis(model; Ecut, kgrid, fft_size, kshift)

    ρ0 = guess_density(basis)
    E, H = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)

    @test E["Hartree"] ≈  0.3527293727197568  atol=5e-8
    @test E["Xc"]      ≈ -2.3033165870558165  atol=5e-8

    # Run one diagonalization and compute energies
    res = diagonalize_all_kblocks(lobpcg_hyper, H, n_bands, tol=1e-9)
    occupation = [[2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]
                  for i in 1:length(basis.kpoints)]
    ρ = compute_density(H.basis, res.X, occupation)
    E, H = energy_hamiltonian(basis, res.X, occupation; ρ)

    @test E["Kinetic"]        ≈  3.3824289861522194  atol=5e-8
    @test E["AtomicLocal"]    ≈ -2.4178712046759157  atol=5e-8
    @test E["AtomicNonlocal"] ≈  1.664289455206788   atol=5e-8
    @test E["Hartree"]        ≈  0.6712993199211524  atol=5e-8
    @test E["Xc"]             ≈ -2.4489960475309056  atol=5e-8
    @test E["Ewald"]          ≈ -8.397893578467201   atol=5e-8
    @test E["PspCorrection"]  ≈ -0.294622067031369   atol=5e-8

    # Now we have a reasonable set of ψ, we make up a crazy model, and check the energies
    V(x, p) = 4*p.ε * ((p.σ/x)^12 - (p.σ/x)^6)
    params = Dict( (:Si, :Si) => (; ε=1e5, σ=0.5) )
    model = model_PBE(silicon.lattice, silicon.atoms, silicon.positions;
                      extra_terms=[ExternalFromReal(X -> cos(1.2 * (X[1] + X[3]))),
                                   ExternalFromFourier(X -> cos(1.3 * (X[1] + X[3]))),
                                   LocalNonlinearity(ρ -> 1.2 * ρ^2.4),
                                   Magnetic(X -> [1, cos(1.4 * X[2]), exp(X[3])]),
                                   PairwisePotential(V, params)],
                      )
    basis = PlaneWaveBasis(model; Ecut, kgrid, fft_size, kshift)
    E, H = energy_hamiltonian(basis, res.X, occupation; ρ)

    @test E["Kinetic"]             ≈  3.3824289861522194  atol=5e-8
    @test E["AtomicLocal"]         ≈ -2.4178712046759157  atol=5e-8
    @test E["AtomicNonlocal"]      ≈  1.664289455206788   atol=5e-8
    @test E["Hartree"]             ≈  0.6712993199211524  atol=5e-8
    @test E["Xc"]                  ≈ -2.469375219486637   atol=5e-8
    @test E["Ewald"]               ≈ -8.397893578467201   atol=5e-8
    @test E["PspCorrection"]       ≈ -0.294622067031369   atol=5e-8
    @test E["ExternalFromReal"]    ≈ -0.01756831422361496 atol=5e-8
    @test E["ExternalFromFourier"] ≈  0.06493077052321815 atol=5e-8
    @test E["LocalNonlinearity"]   ≈  0.14685350034704006 atol=5e-8
    @test E["PairwisePotential"]   ≈ -4.151269801749716   atol=5e-8

    # TODO This is not really a test ... and it does not really work stably.
    # @test E["Magnetic"]            ≈  1.99901120545585e-7 atol=5e-8
end


text:
function gradient_bundle_1st(mechanism::Mechanism, z0, u0, idx;
    N::Int=100,
    Σ=1e-6I,
    opts=SolverOptions(rtol=1e-10, btol=1e-10, undercut=1.0, no_progress_undercut=1.0))

    nx = 6
    nu = 3
    ∇ = zeros(nx,nu)

    for i = 1:N
        # u = normal_sample(u0, Σ)
        u = zeros(3)
        u[idx:idx] = normal_sample(u0[idx:idx], Σ)
        ∇x, ∇u = get_minimal_gradients!(mechanism, z0, u,
            opts=opts)
        # ∇u = (u[idx] / mechanism.timestep >= 10.0) * ones(nx,nu) * mechanism.timestep
        ∇ += ∇u
    end
    return ∇ ./ N
end

function gradient_bundle_0th(mechanism::Mechanism, z0, u0, idx;
    N::Int=100,
    Σ=1e-6*I,
    opts=SolverOptions(rtol=1e-10, btol=1e-10, undercut=1.0, no_progress_undercut=1.0))
    nx = 6
    nu = 1
    δx = [zeros(nx) for i=1:N]
    δu = [zeros(nu) for i=1:N]

    step!(mechanism, z0, u0, opts=opts)
    z10 = get_next_state(mechanism)
    x10 = maximal_to_minimal(mechanism, z10)

    for i = 1:N
        u = zeros(3)
        u[idx:idx] = normal_sample(u0[idx:idx], Σ)
        step!(mechanism, z0, u, opts=opts)
        z1 = get_next_state(mechanism)
        x1 = maximal_to_minimal(mechanism, z1)
        # x11 = copy(x10)
        # x11[idx] = max(0, u[idx]  - 10.0) / 10
        # x11[3+idx] = max(0, u[idx - 10.0)
        δu[i] = (u - u0)[idx:idx]
        δx[i] = x1 - x10
    end
    ∇ = zeros(nx,3)
    ∇[:,idx:idx] = leastsquares(δx, δu)
    return ∇
end


text:
module STC
using ProgressMeter
using Printf
import Base: setindex!,getindex,copy!
export additivestc, extract

function readlsbq end
function readlsb end
function trypmone end
function updatedistortion! end
abstract type AbstractImageDistortion end

ind2sub(image::AbstractImageDistortion, i::Int) = Tuple(CartesianIndices((image.height,image.width))[i])
Base.getindex(image::AbstractImageDistortion,i,j) = image.image[i,j]
Base.getindex(image::AbstractImageDistortion,i) = getindex(image,ind2sub(image, i)...)

include("lsbop.jl")
include("additive.jl")
include("common.jl")
include("randomizedtrellis.jl")
include("dynamic.jl")
include("minchanges.jl")

include("SUniward/SUniward.jl")

end


text:
using Random

struct LazyDataLoader
	preprocess::Function
	datasource::Tuple
	batchsize::Int
	shuffle::Bool
	nobs::Int
	LazyDataLoader(preprocess, datasource; batchsize=1, shuffle=true, partial=false) = begin

		batchsize > 0 || throw(ArgumentError("Batchsize must be a positive number"))

		# Wrap the datasource in a tuple if isn't already
		datasource isa Tuple || (datasource = (datasource,))

		nobs = min(map(x->size(x)[end], datasource)...)
		if nobs < batchsize
			@warn "Number of data points less than batchsize, a single batch of $nobs points will be generated"
			batchsize = nobs
		end
		partial && (nobs -= nobs % batchsize)
		new(preprocess, datasource, batchsize, shuffle, nobs)

	end
end

LazyDataLoader(datasource; kwargs...) = LazyDataLoader(identity, datasource; kwargs...)
Base.length(loader::LazyDataLoader) = loader.nobs

function Base.iterate(loader::LazyDataLoader, (i, iterorder)=(1, nothing))

	i > loader.nobs && return nothing

	# During the first iteration define the access order of observations
	if (i === 1)
		iterorder = [1:loader.nobs;]
		loader.shuffle && shuffle!(iterorder)
	end

	# Get slices from data tuple following iterorder
	currentdata = map(loader.datasource) do data
		return view(
			data,
			(Base.Colon() for dim in 1:ndims(data)-1)...,
			iterorder[i:min(i + loader.batchsize - 1, loader.nobs)]
		)
	end |> loader.preprocess

	return (currentdata, ((i + loader.batchsize), iterorder))

end

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
"""
Created on Sun May 03 2020

@author: Yoann Pradat

Structs used throughout NMF code.
"""

#### # PARAMETERS
#### # #########################################################################################################

abstract type AbstractNMFParams end

"""
    NMFParams([,rank, init, seed, max_iter, stopping_crit, stopping_tol, verbose])

Set of basic parameters passed to the various NMF algorithms.

# Arguments
- `rank::Integer`         : rank of the factorisation
- `init::Symbol`          : initialization method. See [`_initialize_nmf`](@ref)
- `dist::Distribution`    : distribution for random initialization. See [`_initialize_nmf`](@ref)
- `seed::Integer`         : seed for random initializations
- `max_iter::Integer`     : maximum number of iterations (in main procedure)
- `stopping_crit::Symbol` : criterion for stopping iterations. One of :none, :conn, :rel_cost or :abs_cost
- `stopping_iter::Symbol` : number of successive iterations with criterion satisfied for calling convergence
- `stopping_tol::Real`    : tolerance for the stopping crit if crit is :cost
- `verbose::Integer`      : 0 for silent, 1 for warning and 2 and more for to show intermediate information
"""
mutable struct NMFParams <: AbstractNMFParams
    rank::Integer
    init::Symbol
    dist::Distribution
    seed::Integer
    max_iter::Integer
    stopping_crit::Symbol
    stopping_iter::Integer
    stopping_tol::Real
    verbose::Integer

    function NMFParams(;
        rank::Integer          = 5,
        init::Symbol           = :random,
        dist::Distribution     = truncated(Normal(0,1), 0, Inf),
        seed::Integer          = 0,
        max_iter::Integer      = 1000,
        stopping_crit::Symbol  = :none,
        stopping_iter::Integer = 20,
        stopping_tol::Real     = 1e-5,
        verbose::Integer       = 1)

        valid_init = [:default, :random, :nndsvd, :nndsvda, :nndsvdar]
        valid_stopping_crit = [:none, :conn, :rel_cost, :abs_cost]

        rank > 1  || throw(ArgumentError("rank must be at least 2"))
        max_iter > 1  || throw(ArgumentError("max_iter must be at least 2"))
        stopping_tol > 0 || throw(ArgumentError("stopping_tol must be > 0"))
        stopping_iter > 0 || throw(ArgumentError("stopping_iter must be > 0"))
        init in valid_init  || throw(ArgumentError(string("init should be one of ", valid_init)))
        stopping_crit in valid_stopping_crit || throw(ArgumentError(string("stopping_crit should be one of ", valid_stopping_crit)))

        new(rank, init, dist, seed, max_iter, stopping_crit, stopping_iter, stopping_tol, verbose)
    end
end

#### # METRICS
#### # #########################################################################################################

"""
    NMFMetrics(V, W_init, H_init, func_cost, params, local_params)

Struct for tracking cost-related metrics.
"""

mutable struct NMFMetrics
    cur_iter::Integer
    cur_cost::Real
    cur_divg::Real
    cur_conn::Matrix{Int64}
    unchange_conn::Array{Bool,1}
    absolute_cost::Array{Real,1}
    absolute_divg::Array{Real,1}
    relative_cost::Array{Real,1}
    func_cost::Function
    func_divg::Function

    function NMFMetrics(
        V::Matrix{T},
        W_init::Matrix{T},
        H_init::Matrix{T},
        func_cost::Function,
        func_divg::Function) where T <: Real

        unchange_conn = Array{Bool,1}()
        absolute_cost = Array{Real,1}()
        absolute_divg = Array{Real,1}()
        relative_cost = Array{Real,1}()

        cur_iter = 1
        cur_cost = func_cost(V, W_init, H_init)
        cur_divg = func_divg(V, W_init, H_init)
        cur_conn = ones(Int64, size(V, 2), size(V, 2))
        push!(absolute_cost, cur_cost)
        push!(absolute_divg, cur_divg)

        new(cur_iter, cur_cost, cur_divg, cur_conn, unchange_conn, absolute_cost, absolute_divg, relative_cost, func_cost, func_divg)
    end
end

#### # RESULTS
#### # #########################################################################################################

abstract type AbstractNMFResults end

"""
    NMFResults(W, H, metrics, global_params, local_params)

NMFResults returned by NMF algorithms.

# Arguments
- `W             :  : Matrix{T}`         : W matrix after last iteration
- `H             :  : Matrix{T}`         : H matrix after last iteration
- `metrics       :  : NMFMetrics`           : see [`NMFMetrics`](@ref)
- `global_params :  : AbstractNMFParams` : see [`NMFParams`](@ref)
- `local_params  :  : AbstractNMFParams` : see [`FIParams`](@ref), [`CParams`](@ref), [`MUParams`](@ref)
"""
mutable struct NMFResults{T} <: AbstractNMFResults
    W::Matrix{T}
    H::Matrix{T}
    metrics::NMFMetrics
    global_params::AbstractNMFParams
    local_params::AbstractNMFParams
    converged::Bool

    function NMFResults{T}(
        W::Matrix{T},
        H::Matrix{T},
        metrics::NMFMetrics,
        global_params::AbstractNMFParams,
        local_params::AbstractNMFParams,
        converged::Bool) where T <: Real

        new{T}(W, H, metrics, global_params, local_params, converged)
    end
end

#### # STRUCT FOR STORING NMF
#### # #########################################################################################################

mutable struct NMF
    solver::Function
    global_params::AbstractNMFParams
    local_params::AbstractNMFParams

    function NMF(;solver::Function, global_params::AbstractNMFParams, local_params::AbstractNMFParams)
        new(solver, global_params, local_params)
    end
end

#### overload copy function for NMF struct
Base.copy(nmf::NMF) = NMF(;Dict(n => deepcopy(getfield(nmf, n)) for n ∈ fieldnames(NMF))...)


text:
using LightGraphs, SparseArrays, LinearAlgebra

import Base: copy
import LightGraphs: nv, has_edge, add_edge!, rem_edge!, rem_vertex!,
    rem_vertices!, add_vertex!, add_vertices!, outneighbors, inneighbors, neighbors,
    vertices, adjacency_matrix, ne, is_directed, degree, indegree, outdegree, edges,
    has_vertex, all_neighbors

export DiMultigraph

mutable struct DiMultigraph{T<:Integer} <: AbstractMultigraph{T}
    adjlist::Dict{T, Vector{T}}
    _idmax::T
    function DiMultigraph{T}(d::Dict{T, Vector{T}}, _idmax::T) where {T<:Integer}
        adjlist = deepcopy(d)
        vs = keys(adjlist)
        for (v, l) in adjlist
            if l ⊆ vs
                sort!(l)
            else
                error("Some vertices connected to $v is not in the multigraph!")
            end
        end
        _idmax >= 0 || (_idmax = isempty(adjlist) ? 0 : maximum(vs))
        new{T}(adjlist, _idmax)
    end
end

DiMultigraph(adjlist::Dict{T, Vector{T}}) where {T<:Integer} = DiMultigraph{T}(adjlist, isempty(adjlist) ? 0 : maximum(keys(adjlist)))
function DiMultigraph(adjmx::AbstractMatrix{U}) where {U<:Integer}
    m, n = size(adjmx)
    if m != n
        error("Adjacency matrices should be square!")
    end
    if sum(adjmx .!= 0) != sum(adjmx .> 0)
        error("All elements in adjacency matrices should be non-negative!")
    end
    adjlist = Dict(zip((1:m), [Int[] for _ = 1:m]))
    for v1 = 1:m
        for v2 = 1:m
            for i = 1:adjmx[v1, v2]
                push!(adjlist[v1], v2)
            end
        end
    end
    DiMultigraph{Int}(adjlist, m)
end
function DiMultigraph(n::T) where {T<:Integer}
    n >= 0 || error("Number of vertices should be non-negative")
    adjlist = Dict{T, Vector{T}}()
    for i = 1:n
        adjlist[i] = T[]
    end
    return DiMultigraph(adjlist)
end
DiMultigraph(g::SimpleDiGraph{T}) where {T<:Integer} = DiMultigraph(Dict(zip(T(1):nv(g), LightGraphs.SimpleGraphs.fadj(g))))

copy(mg::DiMultigraph{T}) where {T} = DiMultigraph{T}(deepcopy(mg.adjlist), mg._idmax)

nv(mg::DiMultigraph{T}) where {T<:Integer} = T(length(mg.adjlist))
vertices(mg::DiMultigraph) = collect(keys(mg.adjlist))
has_vertex(mg::DiMultigraph, v::Integer) = haskey(mg.adjlist, v)

function adjacency_matrix(mg::DiMultigraph)
    adjmx = spzeros(Int, nv(mg), nv(mg))

    ids = sort!(vertices(mg))
    for id1 in ids
        v1 = searchsortedfirst(ids, id1)
        for id2 in mg.adjlist[id1]
            v2 = searchsortedfirst(ids, id2)
            @inbounds adjmx[v1, v2] += 1
        end
    end
    return adjmx
end

function add_edge!(mg::DiMultigraph, me::AbstractMultipleEdge)
    s = src(me)
    d = dst(me)
    m = mul(me)
    if has_vertex(mg, s) && has_vertex(mg, d)
        for i = 1:m
            insert!(mg.adjlist[s], searchsortedfirst(mg.adjlist[s], d), d)
        end
        return true
    end
    return false
end

function rem_edge!(mg::DiMultigraph, me::AbstractMultipleEdge)
    if has_edge(mg, me)
        s = src(me)
        d = dst(me)
        m = mul(me)
        for i = 1:m
            deleteat!(mg.adjlist[s], searchsortedfirst(mg.adjlist[s], d))
        end
        return true
    else
        return false
    end
end

function rem_vertices!(mg::DiMultigraph{T}, vs::Vector{T}) where {T<:Integer}
    if all(has_vertex(mg, v) for v in vs)
        for v in vs
            for u in neighbors(mg, v)
                l = mg.adjlist[u]
                deleteat!(l, searchsorted(l, v))
            end
            delete!(mg.adjlist, v)
        end
        if mg._idmax in vs
            mg._idmax = maximum(keys(mg.adjlist))
        end
        return true
    end
    return false
end

function add_vertices!(mg::DiMultigraph{T}, n::Integer) where {T<:Integer}
    idmax = mg._idmax
    mg._idmax += n
    new_ids = collect((idmax+1):(idmax+n))
    for i in new_ids
        mg.adjlist[i] = T[]
    end
    return new_ids
end

function outneighbors(mg::DiMultigraph, v::Integer; count_mul::Bool = false)
    has_vertex(mg, v) || error("Vertex not found!")
    if count_mul
        return copy(mg.adjlist[v])
    else
        return sort!(collect(Set(mg.adjlist[v])))
    end
end
function inneighbors(mg::DiMultigraph{T}, v::Integer; count_mul::Bool = false) where T
    has_vertex(mg, v) || error("Vertex not found!")

    innb = T[]
    for u in vertices(mg)
        mul_u_v = length(searchsorted(outneighbors(mg, u), v))
        if mul_u_v > 0
            if count_mul
                for i = 1:mul_u_v
                    push!(innb, u)
                end
            else
                push!(innb, u)
            end
        end
    end
    sort!(innb)
    return innb
end
neighbors(mg::DiMultigraph, v::Integer; count_mul::Bool = false) = outneighbors(mg, v, count_mul = count_mul)
all_neighbors(mg::DiMultigraph, v::Integer) = sort!(union(outneighbors(mg, v), inneighbors(mg, v)))

function mul(mg::DiMultigraph, s::Integer, d::Integer)
    (has_vertex(mg, s) && has_vertex(mg, d)) || error("Vertices not found!")
    return length(searchsorted(mg.adjlist[s], d))
end

is_directed(mg::DiMultigraph) = true
function ne(mg::DiMultigraph; count_mul::Bool = false)
    if count_mul
        return sum([length(mg.adjlist[v]) for v in vertices(mg)])
    else
        return sum([length(Set(mg.adjlist[v])) for v in vertices(mg)])
    end
end

function outdegree(mg::DiMultigraph{T}) where T
    degs = Dict{T, Int}()
    for v in vertices(mg)
        degs[v] = length(mg.adjlist[v])
    end
    return degs
end
function indegree(mg::DiMultigraph{T}) where T
    degs = Dict{T, Int}()
    for v in vertices(mg)
        degs[v] = 0
    end
    for v in vertices(mg)
        for u in mg.adjlist[v]
            degs[u] += 1
        end
    end
    return degs
end
degree(mg::DiMultigraph) = outdegree(mg)
degree(mg::DiMultigraph, v::Integer) = outdegree(mg, v)
indegree(mg::DiMultigraph, v::Integer) = length(inneighbors(mg, v; count_mul = true))
outdegree(mg::DiMultigraph, v::Integer) = length(outneighbors(mg, v; count_mul = true))


text:
# Benchmark the parsing of a MMTF file given as an argument

using BioStructures

mmtf_filepath = ARGS[1]

# Run to JIT compile
read(mmtf_filepath, MMTF)

elapsed = @elapsed struc = read(mmtf_filepath, MMTF)

println(elapsed)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks