text:
################################################################################
#                               VARIABLE DEFINITION
################################################################################
# Define types for different variables
abstract type InfOptVariableType end
struct Infinite <: InfOptVariableType end
struct Point <: InfOptVariableType end
struct Hold <: InfOptVariableType end
struct Deriv <: InfOptVariableType end

# Fallback _make_variable (the methods are defined in the variable files)
function _make_variable(_error::Function, info::JuMP.VariableInfo, type;
                        extra_kw_args...)
    _error("Unrecognized variable type `$(typeof(type).parameters[1])`, " *
           "should be Infinite, Point, or Hold.")
end

"""
    JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                        var_type::Type{<:InfOptVariableType};
                        [parameter_refs::Union{GeneralVariableRef,
                                              AbstractArray{<:GeneralVariableRef},
                                              Tuple, Nothing} = nothing,
                        infinite_variable_ref::Union{GeneralVariableRef,
                                                     Nothing} = nothing,
                        parameter_values::Union{Number, AbstractArray{<:Real},
                                                Tuple, Nothing} = nothing,
                        parameter_bounds::Union{ParameterBounds{GeneralVariableRef},
                                                Nothing} = nothing]
                        )::InfOptVariable

Extend the `JuMP.build_variable` function to accomodate `InfiniteOpt`
variable types. Returns the appropriate variable Datatype (i.e.,
[`InfiniteVariable`](@ref), [`PointVariable`](@ref), and
[`HoldVariable`](@ref)). Primarily, this method is to be used internally by the
appropriate constructor macros [`@infinite_variable`](@ref),
[`@point_variable`](@ref), and [`@hold_variable`](@ref). However, it can be
called manually to build `InfiniteOpt` variables. Errors if an unneeded keyword
argument is given or if the keywoard arguments are formatted incorrectly (e.g.,
`parameter_refs` contains repeated parameter references when an infinite variable
is defined). Also errors if needed keyword arguments are negated.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel())
julia> @independent_parameter(m, t in [0, 1])
t

julia> info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite, parameter_refs = t)
InfiniteVariable{GeneralVariableRef}(VariableInfo{Float64,Float64,Float64,Function}(false, 0.0, false, 0.0, false, 0.0, false, start_func, false, false), (t,), Int64[], Int64[], true)

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)

julia> pt_var = build_variable(error, info, Point, infinite_variable_ref = ivref,
                               parameter_values = 0.5)
PointVariable{GeneralVariableRef}(VariableInfo{Float64,Float64,Float64,Float64}(false, 0.0, false, 0.0, false, 0.0, true, 0.0, false, false), var_name(t), [0.5])

julia> hd_var = build_variable(error, info, Hold)
HoldVariable{GeneralVariableRef}(VariableInfo{Float64,Float64,Float64,Float64}(false, 0.0, false, 0.0, false, 0.0, false, 0.0, false, false), Subdomain bounds (0): )
```
"""
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                             var_type::Type{<:InfOptVariableType};
                             macro_error::Union{Function, Nothing} = nothing,
                             kw_args...)::InfOptVariable
    if macro_error !== nothing
        _error = macro_error # replace with macro error function
    end
    # make the variable and conduct necessary checks
    return _make_variable(_error, info, var_type; kw_args...)
end

# Fallback
function _check_and_make_variable_ref(model::InfiniteModel, v::T, name) where {T}
    throw(ArgumentError("Invalid variable object type `$T`."))
end

# helper function for setting the info constraints (#TODO add test)
function _set_info_constraints(info::JuMP.VariableInfo, gvref::GeneralVariableRef,
                               dvref::DispatchVariableRef)::Nothing
    model = JuMP.owner_model(gvref)
    if info.has_lb
        newset = MOI.GreaterThan(info.lower_bound)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_lower_bound_index(dvref, JuMP.index(cref))
    end
    if info.has_ub
        newset = MOI.LessThan(info.upper_bound)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_upper_bound_index(dvref, JuMP.index(cref))
    end
    if info.has_fix
        newset = MOI.EqualTo(info.fixed_value)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_fix_index(dvref, JuMP.index(cref))
    end
    if info.binary
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, MOI.ZeroOne()),
                                   is_info_constr = true)
        _set_binary_index(dvref, JuMP.index(cref))
    elseif info.integer
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, MOI.Integer()),
                                   is_info_constr = true)
        _set_integer_index(dvref, JuMP.index(cref))
    end
    return
end

"""
    JuMP.add_variable(model::InfiniteModel, var::InfOptVariable,
                      [name::String = ""])::GeneralVariableRef

Extend the [`JuMP.add_variable`](@ref JuMP.add_variable(::JuMP.Model, ::JuMP.ScalarVariable, ::String))
function to accomodate `InfiniteOpt` variable types. Adds a variable to an
infinite model `model` and returns a [`GeneralVariableRef`](@ref).
Primarily intended to be an internal function of the
constructor macros [`@infinite_variable`](@ref), [`@point_variable`](@ref), and
[`@hold_variable`](@ref). However, it can be used in combination with
[`JuMP.build_variable`](@ref) to add variables to an infinite model object.
Errors if invalid parameters reference(s) or an invalid infinite variable
reference is included in `var`.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel())
julia> @infinite_parameter(m, t in [0, 10]);

julia> info = VariableInfo(false, 0, false, 0, false, 0, true, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite, parameter_refs = t);

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)

julia> pt_var = build_variable(error, info, Point, infinite_variable_ref = ivref,
                               parameter_values = 0.5);

julia> pvref = add_variable(m, pt_var, "var_alias")
var_alias

julia> hd_var = build_variable(error, info, Hold);

julia> hvref = add_variable(m, hd_var, "var_name")
var_name
```
"""
function JuMP.add_variable(model::InfiniteModel, var::InfOptVariable,
                           name::String = ""; kwargs...)::GeneralVariableRef
    dvref = _check_and_make_variable_ref(model, var, name; kwargs...)
    vindex = JuMP.index(dvref)
    gvref = _make_variable_ref(model, vindex)
    _set_info_constraints(var.info, gvref, dvref)
    return gvref
end

################################################################################
#                            VARIABLE DEPENDENCIES
################################################################################
# Extend _measure_dependencies
function _measure_dependencies(vref::DecisionVariableRef)::Vector{MeasureIndex}
    return _data_object(vref).measure_indices
end

# Extend _constraint_dependencies
function _constraint_dependencies(vref::DecisionVariableRef)::Vector{ConstraintIndex}
    return _data_object(vref).constraint_indices
end

"""
    used_by_measure(vref::DecisionVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a measure.

**Example**
```julia-repl
julia> used_by_measure(vref)
true
```
"""
function used_by_measure(vref::DecisionVariableRef)::Bool
    return !isempty(_measure_dependencies(vref))
end

"""
    used_by_constraint(vref::DecisionVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a constraint.

**Example**
```julia-repl
julia> used_by_constraint(vref)
false
```
"""
function used_by_constraint(vref::DecisionVariableRef)::Bool
    return !isempty(_constraint_dependencies(vref))
end

"""
    used_by_objective(vref::DecisionVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by the objective.

**Example**
```julia-repl
julia> used_by_objective(vref)
true
```
"""
function used_by_objective(vref::DecisionVariableRef)::Bool
    return _data_object(vref).in_objective
end

"""
    is_used(vref::DecisionVariableRef)::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```julia-repl
julia> is_used(vref)
true
```
"""
function is_used(vref::DecisionVariableRef)::Bool
    return used_by_measure(vref) || used_by_constraint(vref) || used_by_objective(vref)
end

################################################################################
#                                VARIABLE NAMING
################################################################################
"""
    JuMP.name(vref::DecisionVariableRef)::String

Extend [`JuMP.name`](@ref JuMP.name(::JuMP.VariableRef)) to return the names of
`InfiniteOpt` variables.

**Example**
```julia-repl
julia> name(vref)
"var_name"
```
"""
function JuMP.name(vref::DecisionVariableRef)::String
    object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(vref::DecisionVariableRef, name::String)::Nothing

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
names of decision variables.

**Example**
```julia-repl
julia> set_name(vref, "var_name")

julia> name(vref)
"var_name"
```
"""
function JuMP.set_name(vref::DecisionVariableRef, name::String)::Nothing
    _data_object(vref).name = name
    JuMP.owner_model(vref).name_to_var = nothing
    return
end

# Make a variable reference
function _make_variable_ref(model::InfiniteModel, index::ObjectIndex)::GeneralVariableRef
    return GeneralVariableRef(model, index.value, typeof(index))
end

# Get the name_to_var Dictionary
function _var_name_dict(model::InfiniteModel)::Union{Nothing, Dict{String, ObjectIndex}}
    return model.name_to_var
end

# Update name_to_var
function _update_var_name_dict(model::InfiniteModel, var_dict::MOIUC.CleverDict)::Nothing
    name_dict = _var_name_dict(model)
    for (index, data_object) in var_dict
        var_name = data_object.name
        if haskey(name_dict, var_name)
            name_dict[var_name] = HoldVariableIndex(-1) # dumby value
        else
            name_dict[var_name] = index
        end
    end
    model.name_to_var = name_dict
    return
end

"""
    JuMP.variable_by_name(model::InfiniteModel,
                          name::String)::Union{GeneralVariableRef, Nothing}

Extend [`JuMP.variable_by_name`](@ref JuMP.variable_by_name(::JuMP.Model, ::String))
for `InfiniteModel` objects. Return the variable reference assoociated with a
variable name. Errors if multiple variables have the same name. Returns nothing
if no such name exists.

**Examples**
```julia-repl
julia> variable_by_name(m, "var_name")
var_name

julia> variable_by_name(m, "fake_name")

```
"""
function JuMP.variable_by_name(model::InfiniteModel,
                               name::String)::Union{GeneralVariableRef, Nothing}
    if _var_name_dict(model) === nothing
        model.name_to_var = Dict{String, ObjectIndex}()
        _update_var_name_dict(model, model.infinite_vars)
        _update_var_name_dict(model, model.reduced_vars)
        _update_var_name_dict(model, model.point_vars)
        _update_var_name_dict(model, model.hold_vars)
    end
    index = get(_var_name_dict(model), name, nothing)
    if index isa Nothing
        return nothing
    elseif index == HoldVariableIndex(-1)
        error("Multiple variables have the name $name.")
    else
        return _make_variable_ref(model, index)
    end
end

################################################################################
#                              PARAMETER REFERENCES
################################################################################
# Extend parameter_refs for variables (this serves as a fallback for finite types)
function parameter_refs(vref::FiniteVariableRef)::Tuple
    return ()
end

################################################################################
#                         VARIABLE OBJECT MODIFICATION
################################################################################
# Extend _set_core_variable_object
function _set_core_variable_object(vref::DecisionVariableRef,
                                   object::InfOptVariable)::Nothing
    _data_object(vref).variable = object
    return
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Get info
function _variable_info(vref::UserDecisionVariableRef)::JuMP.VariableInfo
    return _core_variable_object(vref).info
end

# Make a temporary constraint reference to use the data accessors
function _temp_constraint_ref(model::InfiniteModel,
                              cindex::ConstraintIndex)::InfOptConstraintRef
    return InfOptConstraintRef(model, cindex)
end

"""
    JuMP.has_lower_bound(vref::UserDecisionVariableRef)::Bool

Extend [`JuMP.has_lower_bound`](@ref JuMP.has_lower_bound(::JuMP.VariableRef)) to
return a `Bool` whether an `InfiniteOpt` variable has a lower bound.

**Example**
```julia-repl
julia> has_lower_bound(vref)
true
```
"""
JuMP.has_lower_bound(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).has_lb

"""
    JuMP.lower_bound(vref::UserDecisionVariableRef)::Float64

Extend [`JuMP.lower_bound`](@ref JuMP.lower_bound(::JuMP.VariableRef)) to return the
lower bound of an `InfiniteOpt` variable. Errors if `vref` doesn't have a lower
bound.

**Example**
```julia-repl
julia> lower_bound(vref)
0.0
```
"""
function JuMP.lower_bound(vref::UserDecisionVariableRef)::Float64
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return _variable_info(vref).lower_bound
end

# Extend to return the index of the lower bound constraint associated with `vref`.
function _lower_bound_index(vref::UserDecisionVariableRef)::ConstraintIndex
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return _data_object(vref).lower_bound_index
end

# Extend to specify the index `cindex` of a lower bound constraint for `vref`.
function _set_lower_bound_index(vref::UserDecisionVariableRef,
                                cindex::Union{ConstraintIndex, Nothing})::Nothing
    _data_object(vref).lower_bound_index = cindex
    return
end

"""
    JuMP.set_lower_bound(vref::UserDecisionVariableRef, lower::Real)::Nothing

Extend [`JuMP.set_lower_bound`](@ref JuMP.set_lower_bound(::JuMP.VariableRef, ::Number))
to specify the lower bound of an `InfiniteOpt` variable `vref`. Errors if `vref`
is fixed.

**Example**
```julia-repl
julia> set_lower_bound(vref, -1)

julia> lower_bound(vref)
-1.0
```
"""
function JuMP.set_lower_bound(vref::UserDecisionVariableRef,
                              lower::Real)::Nothing
    newset = MOI.GreaterThan(convert(Float64, lower))
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, newset)
    model = JuMP.owner_model(vref)
    if JuMP.has_lower_bound(vref)
        cindex = _lower_bound_index(vref)
        cref = _temp_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else
        @assert !JuMP.is_fixed(vref) "$vref is fixed, cannot set lower bound."
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_lower_bound_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(true, convert(Float64, lower),
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.LowerBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef

Extend [`JuMP.LowerBoundRef`](@ref JuMP.LowerBoundRef(::JuMP.VariableRef)) to extract
a constraint reference for the lower bound of `vref`.

**Example**
```julia-repl
var ≥ 0.0
```
"""
function JuMP.LowerBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _lower_bound_index(vref)
    model = JuMP.owner_model(vref)
    return _temp_constraint_ref(model, cindex)
end

"""
    JuMP.delete_lower_bound(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.delete_lower_bound`](@ref JuMP.delete_lower_bound(::JuMP.VariableRef))
to delete lower bound of `vref`. Errors if it doesn't have a lower bound.

**Example**
```julia-repl
julia> delete_lower_bound(vref)

julia> has_lower_bound(vref)
false
```
"""
function JuMP.delete_lower_bound(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.LowerBoundRef(vref))
    _set_lower_bound_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.has_upper_bound(vref::UserDecisionVariableRef)::Bool

Extend [`JuMP.has_upper_bound`](@ref JuMP.has_upper_bound(::JuMP.VariableRef)) to
return a `Bool` whether an `InfiniteOpt` variable has an upper bound.

**Example**
```julia-repl
julia> has_upper_bound(vref)
true
```
"""
JuMP.has_upper_bound(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).has_ub

"""
    JuMP.upper_bound(vref::UserDecisionVariableRef)::Float64

Extend [`JuMP.upper_bound`](@ref JuMP.upper_bound(::JuMP.VariableRef)) to return the
upper bound of an `InfiniteOpt` variable. Errors if `vref` doesn't have a upper bound.

**Example**
```julia-repl
julia> upper_bound(vref)
0.0
```
"""
function JuMP.upper_bound(vref::UserDecisionVariableRef)::Float64
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return _variable_info(vref).upper_bound
end

# Extend to return the index of the upper bound constraint associated with `vref`.
function _upper_bound_index(vref::UserDecisionVariableRef)::ConstraintIndex
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return _data_object(vref).upper_bound_index
end

# Extend to specify the index `cindex` of a upper bound constraint for `vref`
function _set_upper_bound_index(vref::UserDecisionVariableRef,
                                cindex::Union{ConstraintIndex, Nothing})::Nothing
    _data_object(vref).upper_bound_index = cindex
    return
end

"""
    JuMP.set_upper_bound(vref::UserDecisionVariableRef, upper::Real)::Nothing

Extend [`JuMP.set_upper_bound`](@ref JuMP.set_upper_bound(::JuMP.VariableRef, ::Number))
to specify the upper bound of an `InfiniteOpt` variable `vref`. Errors if `vref`
is fixed.

**Example**
```julia-repl
julia> set_upper_bound(vref, 1)

julia> upper_bound(vref)
1.0
```
"""
function JuMP.set_upper_bound(vref::UserDecisionVariableRef,
                              upper::Real)::Nothing
    newset = MOI.LessThan(convert(Float64, upper))
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, newset)
    model = JuMP.owner_model(vref)
    if JuMP.has_upper_bound(vref)
        cindex = _upper_bound_index(vref)
        cref = _temp_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else
        @assert !JuMP.is_fixed(vref) "$vref is fixed, cannot set upper bound."
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_upper_bound_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  true, convert(Float64, upper),
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.UpperBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef

Extend [`JuMP.UpperBoundRef`](@ref JuMP.UpperBoundRef(::JuMP.VariableRef)) to extract
a constraint reference for the upper bound of `vref`.

**Example**
```julia-repl
julia> cref = UpperBoundRef(vref)
var ≤ 1.0
```
"""
function JuMP.UpperBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _upper_bound_index(vref)
    model = JuMP.owner_model(vref)
    return _temp_constraint_ref(model, cindex)
end

"""
    JuMP.delete_upper_bound(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.delete_upper_bound`](@ref JuMP.delete_upper_bound(::JuMP.VariableRef))
to delete the upper bound of `vref`. Errors if it doesn't have an upper bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 >= vref))
julia> delete_upper_bound(vref)

julia> has_upper_bound(vref)
false
```
"""
function JuMP.delete_upper_bound(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.UpperBoundRef(vref))
    _set_upper_bound_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  false, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.is_fixed(vref::UserDecisionVariableRef)::Bool

Extend [`JuMP.is_fixed`](@ref JuMP.is_fixed(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is fixed.

**Example**
```julia-repl
julia> is_fixed(vref)
true
```
"""
JuMP.is_fixed(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).has_fix

"""
    JuMP.fix_value(vref::UserDecisionVariableRef)::Float64

Extend [`JuMP.fix_value`](@ref JuMP.fix_value(::JuMP.VariableRef)) to return the fix
value of an `InfiniteOpt` variable. Errors if variable is not fixed.

**Example**
```julia-repl
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::UserDecisionVariableRef)::Float64
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return _variable_info(vref).fixed_value
end

# Extend to return the index of the fix constraint associated with `vref`.
function _fix_index(vref::UserDecisionVariableRef)::ConstraintIndex
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return _data_object(vref).fix_index
end

# Extend to set the index of the fix constraintassociated with `vref`.
function _set_fix_index(vref::UserDecisionVariableRef,
                        cindex::Union{ConstraintIndex, Nothing})::Nothing
    _data_object(vref).fix_index = cindex
    return
end

"""
    JuMP.fix(vref::UserDecisionVariableRef, value::Real;
             force::Bool = false)::Nothing

Extend [`JuMP.fix`](@ref JuMP.fix(::JuMP.VariableRef, ::Number)) to fix the value of
an `InfiniteOpt` variable. Errors if variable has a lower and/or an upper
bound(s) unless `force = true`.

**Examples**
```julia-repl
julia> fix(vref, 3)

julia> fix_value(vref)
3.0

julia> fix(vref2, 2, force = true)

julia> fix_value(vref2)
2.0
```
"""
function JuMP.fix(vref::UserDecisionVariableRef, value::Real;
                  force::Bool = false)::Nothing
    new_set = MOI.EqualTo(convert(Float64, value))
    model = JuMP.owner_model(vref)
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, new_set)
    if JuMP.is_fixed(vref)  # Update existing fixing constraint.
        cindex = _fix_index(vref)
        cref = _temp_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else  # Add a new fixing constraint.
        if  JuMP.has_upper_bound(vref) ||  JuMP.has_lower_bound(vref)
            if !force
                error("Unable to fix $(vref) to $(value) because it has " *
                      "existing variable bounds. Consider calling " *
                      "`JuMP.fix(variable, value; force=true)` which will " *
                      "delete existing bounds before fixing the variable.")
            end
            if  JuMP.has_upper_bound(vref)
                 JuMP.delete_upper_bound(vref)
            end
            if  JuMP.has_lower_bound(vref)
                 JuMP.delete_lower_bound(vref)
            end
        end
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_fix_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  false, info.upper_bound,
                                                  true, convert(Float64, value),
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.FixRef(vref::UserDecisionVariableRef)::InfOptConstraintRef

Extend [`JuMP.FixRef`](@ref JuMP.FixRef(::JuMP.VariableRef)) to return the constraint
reference of the fix constraint associated with `vref`. Errors `vref` is not
fixed.

**Examples**
```julia-repl
julia> cref = FixRef(vref)
var = 1.0
```
"""
function JuMP.FixRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _fix_index(vref)
    model = JuMP.owner_model(vref)
    return _temp_constraint_ref(model, cindex)
end

"""
    JuMP.unfix(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.unfix`](@ref JuMP.unfix(::JuMP.VariableRef)) to unfix `vref`. Errors
if it is not fixed.

**Example**
```julia-repl
julia> unfix(vref)

julia> is_fixed(vref)
false
```
"""
function JuMP.unfix(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.FixRef(vref))
    _set_fix_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  false, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.start_value(vref::UserDecisionVariableRef)::Union{Nothing, Float64}

Extend [`JuMP.start_value`](@ref JuMP.start_value(::JuMP.VariableRef)) to return
starting value of `InfiniteOpt` variable if it has one. Returns `nothing` otherwise.

**Example**
```julia-repl
julia> start_value(vref)
0.0
```
"""
function JuMP.start_value(vref::UserDecisionVariableRef)::Union{Nothing, Float64}
    if _variable_info(vref).has_start
        return _variable_info(vref).start
    else
        return
    end
end

"""
    JuMP.set_start_value(vref::UserDecisionVariableRef, value::Real)::Nothing

Extend [`JuMP.set_start_value`](@ref JuMP.set_start_value(::JuMP.VariableRef, ::Number))
to specify the start value of `InfiniteOpt` variables.

**Example**
```julia-repl
julia> set_start_value(vref, 1)

julia> start_value(vref)
1.0
```
"""
function JuMP.set_start_value(vref::UserDecisionVariableRef,
                              value::Real)::Nothing
    info = _variable_info(vref)
    set_optimizer_model_ready(JuMP.owner_model(vref), false)
    _update_variable_info(vref,
                          JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                            info.has_ub, info.upper_bound,
                                            info.has_fix, info.fixed_value,
                                            true, Float64(value),
                                            info.binary, info.integer))
    return
end

"""
    JuMP.is_binary(vref::UserDecisionVariableRef)::Bool

Extend [`JuMP.is_binary`](@ref JuMP.is_binary(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is binary.

**Example**
```julia-repl
true
```
"""
JuMP.is_binary(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).binary

# Extend to return the index of the binary constraint associated with `vref`.
function _binary_index(vref::UserDecisionVariableRef)::ConstraintIndex
    if !JuMP.is_binary(vref)
        error("Variable $(vref) is not binary.")
    end
    return _data_object(vref).zero_one_index
end

# Extend to specify the index of the binary constraint associated with `vref`.
function _set_binary_index(vref::UserDecisionVariableRef,
                           cindex::Union{ConstraintIndex, Nothing})::Nothing
    _data_object(vref).zero_one_index = cindex
    return
end

"""
    JuMP.set_binary(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.set_binary`](@ref JuMP.set_binary(::JuMP.VariableRef)) to specify an
`InfiniteOpt` variable as a binary variable. Errors if `vref` is an integer
variable.

**Example**
```julia-repl
julia> set_binary(vref)

julia> is_binary(vref)
true
```
"""
function JuMP.set_binary(vref::UserDecisionVariableRef)::Nothing
    if JuMP.is_binary(vref)
        return
    elseif JuMP.is_integer(vref)
        error("Cannot set the variable_ref $(vref) to binary as it " *
              "is already integer.")
    end
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    cref = JuMP.add_constraint(JuMP.owner_model(vref),
                               JuMP.ScalarConstraint(gvref, MOI.ZeroOne()),
                               is_info_constr = true)
    _set_binary_index(vref, JuMP.index(cref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  true, info.integer))
    return
end

"""
    JuMP.BinaryRef(vref::UserDecisionVariableRef)::InfOptConstraintRef

Extend [`JuMP.BinaryRef`](@ref JuMP.BinaryRef(::JuMP.VariableRef)) to return a
constraint reference to the constraint constrainting `vref` to be binary. Errors
if one does not exist.

**Example**
```julia-repl
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _binary_index(vref)
    model = JuMP.owner_model(vref)
    return _temp_constraint_ref(model, cindex)
end

"""
    JuMP.unset_binary(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.unset_binary`](@ref JuMP.unset_binary(::JuMP.VariableRef)) to unset
`vref` as a binary variable. Errors if it is not binary.

```julia-repl
julia> unset_binary(vref)

julia> is_binary(vref)
false
```
"""
function JuMP.unset_binary(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.BinaryRef(vref))
    _set_binary_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  false, info.integer))
    return
end

"""
    JuMP.is_integer(vref::UserDecisionVariableRef)::Bool

Extend [`JuMP.is_integer`](@ref JuMP.is_integer(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is integer.

**Example**
```julia-repl
julia> is_integer(vref)
true
```
"""
JuMP.is_integer(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).integer

# Extend to return the index of the integer constraintassociated with `vref`.
function _integer_index(vref::UserDecisionVariableRef)::ConstraintIndex
    if !JuMP.is_integer(vref)
        error("Variable $(vref) is not an integer.")
    end
    return _data_object(vref).integrality_index
end

# Extend to specify the index of the integer constraint associated with `vref`.
function _set_integer_index(vref::UserDecisionVariableRef,
                            cindex::Union{ConstraintIndex, Nothing})::Nothing
    _data_object(vref).integrality_index = cindex
    return
end

"""
    JuMP.set_integer(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.set_integer`](@ref JuMP.set_integer(::JuMP.VariableRef)) to specify an
`InfiniteOpt` variable as a integer variable. Errors if `vref` is an binary
variable.

**Example**
```julia-repl
julia> set_integer(vref)

julia> is_integer(vref)
true
```
"""
function JuMP.set_integer(vref::UserDecisionVariableRef)::Nothing
    if JuMP.is_integer(vref)
        return
    elseif JuMP.is_binary(vref)
        error("Cannot set the variable_ref $(vref) to integer as it " *
              "is already binary.")
    end
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    cref = JuMP.add_constraint(JuMP.owner_model(vref),
                               JuMP.ScalarConstraint(gvref, MOI.Integer()),
                               is_info_constr = true)
    _set_integer_index(vref, JuMP.index(cref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, true))
    return
end

"""
    JuMP.IntegerRef(vref::UserDecisionVariableRef)::InfOptConstraintRef

Extend [`JuMP.IntegerRef`](@ref JuMP.IntegerRef(::JuMP.VariableRef)) to return a
constraint reference to the constraint constrainting `vref` to be integer.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _integer_index(vref)
    model = JuMP.owner_model(vref)
    return _temp_constraint_ref(model, cindex)
end

"""
    JuMP.unset_integer(vref::UserDecisionVariableRef)::Nothing

Extend [`JuMP.unset_integer`](@ref JuMP.unset_integer(::JuMP.VariableRef)) to unset
`vref` as an integer variable. Errors if it is not an integer variable.

```julia-repl
julia> unset_integer(vref)

julia> is_integer(vref)
false
```
"""
function JuMP.unset_integer(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.IntegerRef(vref))
    _set_integer_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, false))
    return
end

################################################################################
#                          MODEL VARIABLE QUERIES
################################################################################
"""
    JuMP.num_variables(model::InfiniteModel,
                       [type::Type{InfOptVariable} = InfOptVariable])::Int

Extend [`JuMP.num_variables`](@ref JuMP.num_variables(::JuMP.Model)) to return the
number of `InfiniteOpt` variables assigned to `model`. By default, the total
number of infinite, reduced, point, and hold variables is returned. The amount
of a particular type is obtained by specifying the concrete variable type
of [`InfOptVariable`](@ref) via `type`. Type options include:
 - `InfOptVariable`: all variables
 - `InfiniteVariable`: all infinite variables
 - `ReducedVariable`: all reduced variables
 - `PointVariable`: all point variables
 - `HoldVariable`: all hold variables

**Example**
```julia-repl
julia> num_variables(model)
3

julia> num_variables(model, InfiniteVariable)
2
```
"""
function JuMP.num_variables(model::InfiniteModel,
                            type::Type{InfOptVariable} = InfOptVariable
                            )::Int
    num_vars = JuMP.num_variables(model, InfiniteVariable)
    num_vars += JuMP.num_variables(model, ReducedVariable)
    num_vars += JuMP.num_variables(model, PointVariable)
    num_vars += JuMP.num_variables(model, HoldVariable)
    return num_vars
end

# Particular variable types
function JuMP.num_variables(model::InfiniteModel,
                            type::Type{C})::Int where {C <: InfOptVariable}
    return length(_data_dictionary(model, type))
end

"""
    JuMP.all_variables(model::InfiniteModel,
                       type::Type{InfOptVariable} = InfOptVariable
                       )::Vector{GeneralVariableRef}

Extend [`JuMP.all_variables`](@ref JuMP.all_variables(::JuMP.Model)) to return a
list of all the variable references associated with `model`. By default, all
of the infinite, reduced, point, and hold variables is returned. Those
of a particular type is obtained by specifying the concrete variable type
of [`InfOptVariable`](@ref) via `type`. Type options include:
 - `InfOptVariable`: all variables
 - `InfiniteVariable`: all infinite variables
 - `ReducedVariable`: all reduced variables
 - `PointVariable`: all point variables
 - `HoldVariable`: all hold variables

**Examples**
```julia-repl
julia> all_variables(model)
4-element Array{GeneralVariableRef,1}:
 y(t)
 w(t, x)
 y(0)
 z

julia> all_variables(model, PointVariable)
1-element Array{GeneralVariableRef,1}:
 y(0)
```
"""
function JuMP.all_variables(model::InfiniteModel,
                            type::Type{InfOptVariable} = InfOptVariable
                            )::Vector{GeneralVariableRef}
    vrefs_list = JuMP.all_variables(model, InfiniteVariable)
    append!(vrefs_list, JuMP.all_variables(model, ReducedVariable))
    append!(vrefs_list, JuMP.all_variables(model, PointVariable))
    append!(vrefs_list, JuMP.all_variables(model, HoldVariable))
    return vrefs_list
end

# Particular variable types
function JuMP.all_variables(model::InfiniteModel,
                            type::Type{C}
                            )::Vector{GeneralVariableRef} where {C <: InfOptVariable}
    vrefs_list = Vector{GeneralVariableRef}(undef, JuMP.num_variables(model, type))
    for (i, (index, _)) in enumerate(_data_dictionary(model, type))
        vrefs_list[i] = _make_variable_ref(model, index)
    end
    return vrefs_list
end

################################################################################
#                                 DELETION
################################################################################
# Helper function to delete the info constraints
function _delete_info_constraints(vref::UserDecisionVariableRef)::Nothing
    # remove variable info constraints associated with vref
    if JuMP.has_lower_bound(vref)
        JuMP.delete_lower_bound(vref)
    end
    if JuMP.has_upper_bound(vref)
        JuMP.delete_upper_bound(vref)
    end
    if JuMP.is_fixed(vref)
        JuMP.unfix(vref)
    end
    if JuMP.is_binary(vref)
        JuMP.unset_binary(vref)
    elseif JuMP.is_integer(vref)
        JuMP.unset_integer(vref)
    end
    return
end

"""
    JuMP.delete(model::InfiniteModel, vref::DecisionVariableRef)::Nothing

Extend [`JuMP.delete`](@ref JuMP.delete(::JuMP.Model, ::JuMP.VariableRef)) to delete
`InfiniteOpt` variables and their dependencies. Errors if variable is invalid,
meaning it has already been deleted or it belongs to another model.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]
 g(0.5) = 0

julia> delete(model, g)

julia> print(model)
Min measure(t) + z
Subject to
 z ≥ 0.0
 z ≥ 42.0
```
"""
function JuMP.delete(model::InfiniteModel, vref::DecisionVariableRef)::Nothing
    @assert JuMP.is_valid(model, vref) "Variable is invalid."
    # update the optimizer model status
    if is_used(vref)
        set_optimizer_model_ready(model, false)
    end
    # delete attributes specific to the variable type
    _delete_variable_dependencies(vref)
    gvref = _make_variable_ref(model, JuMP.index(vref))
    # remove from measures if used
    for mindex in _measure_dependencies(vref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        data = measure_data(mref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(mref, new_meas)
    end
    # remove from constraints if used
    for cindex in _constraint_dependencies(vref)
        cref = _temp_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
        else
            _remove_variable(func, gvref)
            # update the object numbers if vref is infinite
            if vref isa Union{InfiniteVariableRef, ReducedVariableRef}
                _data_object(cref).object_nums = sort(_object_numbers(func))
            end
        end
    end
    # remove from objective if vref is in it
    if used_by_objective(vref)
        if JuMP.objective_function(model) isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            JuMP.set_objective_function(model, new_func)
            JuMP.set_objective_sense(model, MOI.FEASIBILITY_SENSE)
        else
            _remove_variable(JuMP.objective_function(model), gvref)
        end
    end
    # delete the variable information
    _delete_data_object(vref)
    return
end

################################################################################
#                           INTEGRALITY RELAXATION
################################################################################
"""
    JuMP.relax_integrality(model::InfiniteModel)::Function

Modifies `model` to "relax" all binary and integrality constraints on
variables. Specifically,

- Binary constraints are deleted, and variable bounds are tightened if
  necessary to ensure the variable is constrained to the interval ``[0, 1]``.
- Integrality constraints are deleted without modifying variable bounds.
- All other constraints are ignored (left in place). This includes discrete
  constraints like SOS and indicator constraints.

Returns a function that can be called without any arguments to restore the
original model. The behavior of this function is undefined if additional
changes are made to the affected variables in the meantime.

**Example**
```julia-repl
julia> undo_relax = relax_integrality(model);

julia> print(model)
Min x + ∫{t ∈ [0, 10]}(y(t))
Subject to
 x ≥ 0.0
 y(t) ≥ 1.0
 x ≤ 1.0
 y(t) ≤ 10.0

julia> undo_relax()

julia> print(model)
Min x + ∫{t ∈ [0, 10]}(y(t))
Subject to
 y(t) ≥ 1.0
 y(t) ≤ 10.0
 y(t) integer
 x binary
```
"""
function JuMP.relax_integrality(model::InfiniteModel)
    # TODO ensure variables are not semi-continous/integer
    # get the variable info
    info_pre_relaxation = map(v -> (dispatch_variable_ref(v),
                              _variable_info(dispatch_variable_ref(v))),
                              JuMP.all_variables(model))
    # relax the appropriate variables
    for (vref, info) in info_pre_relaxation
        if info.integer
            JuMP.unset_integer(vref)
        elseif info.binary
            JuMP.unset_binary(vref)
            if isnan(info.lower_bound)
                JuMP.set_lower_bound(vref, 0.0)
            else
                JuMP.set_lower_bound(vref, max(0.0, info.lower_bound))
            end
            if isnan(info.upper_bound)
                JuMP.set_upper_bound(vref, 1.0)
            else
                JuMP.set_upper_bound(vref, min(1.0, info.upper_bound))
            end
        end
    end
    function unrelax()::Nothing
        for (vref, info) in info_pre_relaxation
            if info.integer
                JuMP.set_integer(vref)
            elseif info.binary
                JuMP.set_binary(vref)
                if info.has_lb
                    JuMP.set_lower_bound(vref, info.lower_bound)
                else
                    JuMP.delete_lower_bound(vref)
                end
                if info.has_ub
                    JuMP.set_upper_bound(vref, info.upper_bound)
                else
                    JuMP.delete_upper_bound(vref)
                end
            end
        end
        return
    end
    return unrelax
end


text:
using Bibliography
using ReferenceTests

@testset "cff" begin
    bib = Bibliography.import_cff("../examples/CITATION.cff")
    rm("CITATION.cff", force=true)
    Bibliography.export_cff(bib)
    # TODO - generate a test file for 32 bits architecture
    if Sys.WORD_SIZE == 64
        @test_reference "CITATION.cff.txt" read("CITATION.cff", String)
    end
    rm("CITATION.cff")
end

ext: jl


text:
@doc raw"""
    cyclic_proximal_point(M, F, proxes, x)

perform a cyclic proximal point algorithm.

# Input
* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `proxes` – an Array of proximal maps (`Function`s) `(λ,x) -> y` for the summands of $F$
* `x` – an initial value $x ∈ \mathcal M$

# Optional
the default values are given in brackets
* `evaluationOrder` – ( [`LinearEvalOrder`](@ref) ) – whether
  to use a randomly permuted sequence ([`FixedRandomEvalOrder`](@ref)), a per
  cycle permuted sequence ([`RandomEvalOrder`](@ref)) or the default linear one.
* `λ` – ( `iter -> 1/iter` ) a function returning the (square summable but not
  summable) sequence of λi
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(5000),`[`StopWhenChangeLess`](@ref)`(10.0^-8))`) a [`StoppingCriterion`](@ref).
* `return_options` – (`false`) – if actiavated, the extended result, i.e. the
  complete [`Options`](@ref) are returned. This can be used to access recorded values.
  If set to false (default) just the optimal value `xOpt` if returned
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function cyclic_proximal_point(
    M::MT,
    F::TF,
    proxes::Union{Tuple,AbstractVector},
    x0;
    evaluationOrder::EvalOrder=LinearEvalOrder(),
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(5000), StopWhenChangeLess(10.0^-12)
    ),
    λ=i -> 1 / i,
    return_options=false,
    kwargs..., #decorator options
) where {MT<:Manifold,TF}
    p = ProximalProblem(M, F, proxes)
    o = CyclicProximalPointOptions(x0, stopping_criterion, λ, evaluationOrder)

    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(p::ProximalProblem, o::CyclicProximalPointOptions)
    c = length(p.proxes)
    o.order = collect(1:c)
    return update_cpp_order!(o.order, c, 0, o.orderType)
end
function step_solver!(p::ProximalProblem, o::CyclicProximalPointOptions, iter)
    c = length(p.proxes)
    λi = o.λ(iter)
    for k in o.order
        o.x = get_proximal_map(p, λi, o.x, k)
    end
    return update_cpp_order!(o.order, c, iter, o.orderType)
end
get_solver_result(o::CyclicProximalPointOptions) = o.x
update_cpp_order!(o, n, i, ::LinearEvalOrder) = o
update_cpp_order!(o, n, i, ::RandomEvalOrder) = shuffle!(o)
update_cpp_order!(o, n, i, ::FixedRandomEvalOrder) = (i == 0) ? shuffle!(o) : o


text:
##############################
# Observation-based learning #
##############################
"""
Observation-based network learning model type.
"""
mutable struct NetworkLearnerObs{T,U,S,V,
				    R<:Vector{<:AbstractRelationalLearner},
				    C<:AbstractCollectiveInferer,
				    A<:Vector{<:AbstractAdjacency},
				    L<:Union{Nothing, MLLabelUtils.LabelEncoding},
				    OD<:LearnBase.ObsDimension} <: AbstractNetworkLearner
	Ml::T											# local model
	fl_exec::U										# local model execution function
	Mr::S											# relational model
	fr_exec::V										# relational model execution function
	RL::R											# relational learner
	Ci::C											# collective inferer
	Adj::A											# adjacency information
	use_local_data::Bool							# whether to use local data
	target_enc::L									# target encoding
	size_in::Int									# expected input dimensionality
	size_out::Int									# expected output dimensionality
	obsdim::OD										# observation dimension
end

# Printers
Base.show(io::IO, m::NetworkLearnerObs) = begin
	println(io,"NetworkLearner, I/O size $(m.size_in)×$(m.size_out), observation based")
	print(io,"`- local model: "); println(io, m.Ml)
	print(io,"`- relational model: "); println(io, m.Mr)
	print(io,"`- relational learners: "); println(io, m.RL)
	print(io,"`- collective inferer: "); println(io, m.Ci)
	print(io,"`- adjacency: "); println(io, m.Adj)
	println(io,"`- use local data: $(m.use_local_data)");
	println(io,"`- observations are: $(m.obsdim == ObsDim.Constant{1} ? "rows" : "columns")");
	println(io,"`- targets: $(m.target_enc isa Nothing ? "not encoded" : "encoded")");
end

####################
# Training methods #
####################
"""
	fit(::Type{NetworkLearnerObs}, X, y, Adj, fl_train, fl_exec, fr_train, fr_exec [;kwargs])

Training method for the observation-based network learning framework.

# Arguments
  * `X::AbstractMatrix` training data (used by `fl_train`, `fl_exec`; if `use_local_data==true`, it is also used by `fr_train`)
  * `Y::AbstractAray` data targets (used by `fl_train`, `fr_train`)
  * `Adj::Vector{AbstractAdjacency}` a vector containing the observation relational structures (adjacency objects)
  * `fl_train` local model training 'function'; can be anything that supports the call `fl_train((X,y))`
  * `fl_exec` local model prediction 'function'; can be anything that supports the call `fl_exec(Ml,X)` where `Ml = fl_train((X,y))`
  * `fr_train` relational model training `function`; can be anything that suports the call `fr_train((Xr,y))`
  * `fr_exec` relational model prediction `function`; can be anything that suports the call `fr_exec(Mr,Xr)` where `Mr = fr_train((Xr,y))`
  and `Xr` is a dataset of relational variables generated by the relational learner using the results of the local model prediction
  function and adjacency structures.

# Keyword arguments
  * `priors::Vector{Float64}` class priors (if applicable)
  * `learner::Symbol` relational learner (i.e. variable generator); available options `:rn`, `:wrn`, `:bayesrn` and `:cdrn` (default `:wrn`)
  * `inference::Symbol` collective inference method; available options `:rl`, `:ic` and `:gs` (default `:rl`)
  * `normalize::Bool` whether to normalize the relational variables per-observation to the L1 norm (default `true`)
  * `use_local_data::Bool` whether the relational model should use the local data provided (i.e. in `X`) (default `true`)
  * `f_targets::Function` function that extracts targets from estimates generated by the local/relational models
  (default `f_targets = x->MLDataPattern.targets(indmax,x)`)
  * `obsdim::Int` observation dimension (default `2`)
  * `tol::Float64` maximum admissible mean estimate error for collective inference convergence (default `1e-6`)
  * `κ::Float64` relaxation labeling starting constant, used if `learner == :rl` (default `1.0`)
  * `α::Float64` relaxation labeling decay constant, used if `learner == :rl` (default `0.99`)
  * `maxiter::Int` maximum number of iterations for collective inference (default `100`)
  * `bratio::Float64` percentage of iterations i.e. `maxiter` used for Gibbs sampling burn-in (default `0.1`)
"""
function fit(::Type{NetworkLearnerObs}, X::AbstractMatrix, y::AbstractArray, Adj::A where A<:Vector{<:AbstractAdjacency},
		fl_train, fl_exec, fr_train, fr_exec;
		priors::Vector{Float64}=getpriors(y), learner::Symbol=:wvrn, inference::Symbol=:rl,
		normalize::Bool=true, use_local_data::Bool=true, f_targets::Function=x->targets(indmax,x),
		obsdim::Int = 2, tol::Float64=1e-6, κ::Float64=1.0, α::Float64=0.99, maxiter::Int=100, bratio::Float64=0.1)

	# Parse, transform input arguments
	κ = clamp(κ, 1e-6, 1.0)
	α = clamp(α, 1e-6, 1.0-1e-6)
	tol = clamp(tol, 0.0, Inf)
	maxiter = ifelse(maxiter<=0, 1, maxiter)
	bratio = clamp(bratio, 1e-6, 1.0-1e-6)

	@assert all((priors.>=0.0) .& (priors .<=1.0)) "All priors have to be between 0.0 and 1.0."
	@assert obsdim in [1,2] "Observation dimension can have only two values 1 (row-major) or 2 (column-major)."

	# Parse relational learner argument and generate relational learner type
	if learner == :rn
		Rl = SimpleRN
	elseif learner == :wrn
		Rl = WeightedRN
	elseif learner == :cdrn
		Rl = ClassDistributionRN
	elseif learner == :bayesrn
		Rl = BayesRN
	else
		@print_verbose 1 "Unknown relational learner. Defaulting to :wrn."
		Rl = WeightedRN
	end

	# Parse collective inference argument and generate collective inference objects
	if inference == :rl
		Ci = RelaxationLabelingInferer(maxiter, tol, f_targets, κ, α)
	elseif inference == :ic
		Ci = IterativeClassificationInferer(maxiter, tol, f_targets)
	elseif inference == :gs
		Ci = GibbsSamplingInferer(maxiter, tol, f_targets, ceil(Int, maxiter*bratio))
	else
		@print_verbose 1 "Unknown collective inferer. Defaulting to :rl."
		Ci = RelaxationLabelingInferer(maxiter, tol, f_targets, κ, α)
	end

	fit(NetworkLearnerObs, X, y, Adj, Rl, Ci, fl_train, fl_exec, fr_train, fr_exec;
		priors=priors, normalize=normalize, use_local_data=use_local_data,
		obsdim=ObsDim.Constant{obsdim}())
end

"""
Training method for the network learning framework. This method should not be called directly.
"""
function fit(::Type{NetworkLearnerObs}, X::T, y::S, Adj::A, Rl::R, Ci::C, fl_train::U, fl_exec::U2, fr_train::U3, fr_exec::U4;
		priors::Vector{Float64}=getpriors(y), normalize::Bool=true, use_local_data::Bool=true,
		obsdim::LearnBase.ObsDimension=ObsDim.Constant{2}()) where {
			T<:AbstractMatrix,
			S<:AbstractArray,
			A<:Vector{<:AbstractAdjacency},
			R<:Type{<:AbstractRelationalLearner},
			C<:AbstractCollectiveInferer,
			U, U2, U3, U4
		}

	# Step 0: pre-process input arguments and retrieve sizes
	size_in = nvars(X,obsdim)								# number of local variables
	size_out = get_size_out(y)								# number of relational variables / adjacency
	n = nobs(X,obsdim)									# number of observations
	m = length(Adj) * size_out								# number of relational variables

	@assert size_out == length(priors) "Found $c classes, the priors indicate $(length(priors))."

	# Pre-allocate relational variables array
	if use_local_data									# Local observation variable data is used
		Xr = matrix_prealloc(n, size_in+m, obsdim, 0.0)
		_Xr = datasubset(Xr, 1:size_in, oppdim(obsdim))
		_Xr[:] = X
		offset = size_in
	else											# Only relational variables are used
		Xr = matrix_prealloc(n, m, obsdim, 0.0)
		offset = 0
	end

	# Encode targets
	(t_enc, yₑ) = encode_targets(y)

	# Step 1: train and execute local model
	@print_verbose 2 "Training and applying local model ..."
	Dl = (X,yₑ)
	Ml = fl_train(Dl);
	Xl = fl_exec(Ml,X);

	# Step 2: Get relational variables by training and executing the relational learner
	@print_verbose 2 "Calculating relational variables ..."
	RL = [fit(Rl, Aᵢ, Xl, yₑ; obsdim=obsdim, priors=priors, normalize=normalize)
       		for Aᵢ in Adj]									# Train relational learners

	Xrᵢ = matrix_prealloc(n, size_out, obsdim, 0.0)						# Initialize temporary storage
	for (i,(RLᵢ,Aᵢ)) in enumerate(zip(RL,Adj))

		# Apply relational learner
		transform!(Xrᵢ, RLᵢ, Aᵢ, Xl, yₑ)

		# Update relational data output
		_Xr = datasubset(Xr, offset+(i-1)*size_out+1:offset+i*size_out, oppdim(obsdim))
		_Xr[:] = Xrᵢ
	end

	# Step 3 : train relational model
	@print_verbose 2 "Training relational model ..."
	Mr = fr_train((Xr,yₑ))

	# Step 4: remove adjacency data
	Adj_s = AbstractAdjacency[];
	for i in 1:length(Adj)
		push!(Adj_s, strip_adjacency(Adj[i]))
	end

	# Step 5: return network learner
	@print_verbose 2 "Done."
	return NetworkLearnerObs(Ml, fl_exec, Mr, fr_exec, RL, Ci, Adj_s, use_local_data, t_enc, size_in, size_out, obsdim)
end

#####################
# Execution methods #
#####################
"""
Prediction method for the network learning framework.
"""
function predict(model::M, X::T, update::BitVector=trues(nobs(X,model.obsdim))) where {M<:NetworkLearnerObs, T<:AbstractMatrix}
	Xo = matrix_prealloc(nobs(X,model.obsdim), model.size_out, model.obsdim, 0.0)
	predict!(Xo, model, X, update)
	return Xo
end

"""
In-place prediction method for the network learning framework.
"""
function predict!(Xo::S, model::M, X::T, update::BitVector=trues(nobs(X,model.obsdim))) where {M<:NetworkLearnerObs, T<:AbstractMatrix, S<:AbstractMatrix}

	# Step 0: Make initializations and pre-allocations
	obsdim = model.obsdim									# observation dimension
	m = nvars(X,obsdim)									# number of input variables
	n = nobs(X,obsdim)									# number of observations
	l = length(model.Adj)									# number of adjacencies in the model

	@assert model.size_in == m "Expected input dimensionality $(model.size_in), got $m."
	@assert (nvars(Xo,obsdim),nobs(Xo,obsdim)) == (model.size_out, n) "Output dataset size does not match expected size ($(model.size_out) variables×$n observations)."

	if model.use_local_data									# Pre-allocate relational dataset based on local data use
		Xr = matrix_prealloc(n, model.size_out*l+m, obsdim, 0.0)
		_Xr = datasubset(Xr, 1:m, oppdim(obsdim))
		_Xr[:] = X									#  - initialize (local) dimensions
		offset = m
	else											# Skip local dataset dimensions
		Xr = matrix_prealloc(n, model.size_out*l, obsdim, 0.0)
		offset = 0
	end

	# Step 1: Apply local model, initialize output
	@print_verbose 2 "Executing local model ..."
	Xl = model.fl_exec(model.Ml, getobs(datasubset(X, update, obsdim)))
	@assert nvars(Xo,obsdim) == nvars(Xl,obsdim) "Local model outputs $(nvars(Xl,obsdim)) estimates/obs, the output indicates $(nvars(Xo,obsdim))."
	_Xo = datasubset(Xo, update, obsdim);
	_Xo[:] = Xl

	# Step 2: Apply collective inference
	@print_verbose 2 "Collective inference ..."
	transform!(Xo, model.Ci, model.obsdim, model.Mr, model.fr_exec, model.RL, model.Adj, offset, Xr, update)

	# Step 3: Return output estimates
	@print_verbose 2 "Done."
	return Xo
end

# It may be necessary to add adjacency information to the model, regarding the test data
"""
	add_adjacency!(model, Av)

Function that adds or, creates and then adds, adjacency objects contained in a vector `Av` to a network
learning `model`. The method is used in 'ouf-of-graph' learning i.e. the training data is not used in
the predictions for new samples. In such circumstances, the adjacency structures - graphs, matrices etc.
pertinent to the new observations have to be created either using the same functions used in training or,
separately supplied, so that relational variables can be created.
"""
function add_adjacency!(model::M, Av::Vector{T}) where {M<:NetworkLearnerObs, T<:AbstractAdjacency}
	@assert length(Av) == length(model.Adj) "New adjacency vector must have a length of $(length(model.Adj))."
	model.Adj = Av
end

function add_adjacency!(model::M, Av::Vector{T}) where {M<:NetworkLearnerObs, T}
	@assert length(Av) == length(model.Adj) "Adjacency data vector must have a length of $(length(model.Adj))."
	model.Adj = adjacency.(Av)
end


text:
using Boltzmann
using MLDatasets

if :integration in TEST_GROUPS
    @testset "Nets Integration" begin
        X, y = MNIST.traindata()
	X = Float64.(reshape(X, 784, :)[:, 1:1000])    # take only 1000 observations for speed
        X = X / (maximum(X) - (minimum(X)))  # normalize to [0..1]

        layers = [("vis", GRBM(784, 256)),
                  ("hid1", BernoulliRBM(256, 100)),
                  ("hid2", BernoulliRBM(100, 100))]
        dbn = DBN(layers)
        fit(dbn, X)
        @assert size(transform(dbn, X)) == (100, 1000)

        dae = unroll(dbn)
        transform(dae, X)

        save_params("test.hdf5", dbn)
        save_params("test2.hdf5", dae)
        rm("test.hdf5")
        rm("test2.hdf5")
    end
end


text:
fs()

@testset "Scope (#412)" begin
    write(joinpath(td, "config.md"), """
        @def title = "config"
        """)
    write(joinpath(td, "page.md"), """
        @def title = "page"
        """)
    write(joinpath(td, "index.html"), """
        {{insert head.html}}
        """)
    mkpath(joinpath(td, "_layout"))
    mkpath(joinpath(td, "_css"))
    write(joinpath(td, "_layout", "head.html"), "
{{fill title}}
")
    write(joinpath(td, "_layout", "page_foot.html"), "")
    write(joinpath(td, "_layout", "foot.html"), "")
    serve(single=true)

    @test startswith(read(joinpath("__site", "index.html"), String),
                        "
config
")
    @test startswith(read(joinpath("__site", "page", "index.html"), String),
                        "
page
")
end


text:
mutable struct SparseSymmetricMatrix{T<:Number}
    entries::Dict{Tuple{Any, Any}, T}
    indices::Set{Any}
end

SparseSymmetricMatrix(T::Type) = SparseSymmetricMatrix{T}(Dict{Tuple{Any,Any},T}(), Set())

SparseSymmetricMatrix() = SparseSymmetricMatrix(Float64)

entries(m::SparseSymmetricMatrix) = m.entries

indices(m::SparseSymmetricMatrix) = m.indices

function setindex!{T<:Number}(m::SparseSymmetricMatrix{T}, v::T, i, j)
    entries(m)[(i, j)] = v
    push!(indices(m), i)
    push!(indices(m), j)
end

function getindex{T<:Number}(m::SparseSymmetricMatrix{T}, i, j)
    get(m.entries, (i, j), zero(T))
end

start(m::SparseSymmetricMatrix) = start(m.entries)
done(m::SparseSymmetricMatrix, state) = done(m.entries, state)
next(m::SparseSymmetricMatrix, state) = next(m.entries, state)

function size(m::SparseSymmetricMatrix)
    is = Set()
    js = Set()
    for (i, j) in keys(m.entries)
        add!(is, i)
        add!(js, j)
    end
    return max(length(is), length(js))
end

copy(m::SparseSymmetricMatrix) = SparseSymmetricMatrix(copy(entries(m)), copy(indices(m)))

function *{T<:Number}(x::T, m::SparseSymmetricMatrix{T})
    p = copy(m)
    for k in keys(entries(m))
        entries(p)[k] *= x
    end
    p
end

function delete!{T<:Number}(m::SparseSymmetricMatrix{T}, i, j)
    delete!(entries(m), (i, j))
    igone = true
    jgone = true
    for (ti, tj) in keys(entries(m))
        if i == ti
            igone = false
        end
        if j == tj
            jgone = false
        end
    end
    if igone
        delete!(indices(m), i)
    end
    if jgone && contains(indices(m), j)
        delete!(indices(m), j)
    end
end

function ==(A::SparseSymmetricMatrix, B::SparseSymmetricMatrix)
    A.entries == B.entries && A.indices == B.indices
end

function isapprox(A::SparseSymmetricMatrix, B::SparseSymmetricMatrix)
    if indices(A) != indices(B)
        return false
    end
    for ((ri, ci), v) in A
        if !isapprox(v, B[ri, ci])
            return false
        end
    end
    for ((ri, ci), v) in B
        if !isapprox(v, A[ri, ci])
            return false
        end
    end
    true
end


text:
abstract type MayUnroll end

struct Unroll <: MayUnroll end

struct NoUnroll <: MayUnroll end


text:
# Many test files are leading
# to precompile command that
# cause errors, so they
# are commented out

for filename in (
  "tagset.jl",
  "smallstring.jl",
  "index.jl",
  "indexset.jl",
  "not.jl",
  #"itensor_dense.jl",
  #"itensor_diag.jl",
  #"contract.jl",
  #"combiner.jl",
  #"trg.jl",
  #"ctmrg.jl",
  #"iterativesolvers.jl",
  #"dmrg.jl",
  #"tag_types.jl",
  #"phys_site_types.jl",
  #"decomp.jl",
  #"lattices.jl",
  #"mps.jl",
  #"mpo.jl",
  #"autompo.jl",
  #"svd.jl",
  #"qn.jl",
  #"qnindex.jl",
  #"itensor_blocksparse.jl",
  #"itensor_diagblocksparse.jl",
  #"readwrite.jl",
  #"readme.jl",
  #"examples.jl",
)
  include("../../test/$filename")
end


text:
# skip the first prime (2) and try only odd numbers. 10000. item is now 10001. item :)

function buildPrimes(count)
    primes = Array(Int64, count)
    slot = 1

    function check(n)
        i = 1
        while i < slot
            if n%primes[i] == 0
                return false
            end
            i += 1
        end
        return true
    end

    n = 3
    while slot <= count
        if check(n)
            primes[slot] = n
            slot += 1
        end
        n += 2
    end

    return primes
end

println(buildPrimes(10000)[end])

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
function helloWord()
    print("Hello\n")
    print("The\n")
    print("World\n")
end

helloWord()

max_stars_repo_name: PushpneetSingh/Hello-world

max_issues_repo_name: PushpneetSingh/Hello-world

max_forks_repo_name: PushpneetSingh/Hello-world


text:
include("../benchmark.jl")
include("create_report.jl")

overlapping_results = get_CUTEst_results()

its, best, ratios, times = compute_its_etc(overlapping_results,MAX_IT=3000);

using PyPlot
PyPlot.close()
plot_iteration_ratios(its, best, ratios)
savefig("$folder/inf_iter_ratios.pdf")


text:
mutable struct LDA <: TopicModel
	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	alpha::Vector{Float64}
	beta::Matrix{Float64}
	gamma::VectorList{Float64}
	phi::Matrix{Float64}
	Elogtheta::Vector{Float64}
	Elogthetasum::Vector{Float64}
	newbeta::Matrix{Float64}
	elbo::Float64
	newelbo::Float64

	function LDA(corp::Corpus, K::Integer)
		@assert ispositive(K)
		@assert !isempty(corp)
		checkcorp(corp)

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		alpha = ones(K)
		beta = rand(Dirichlet(V, 1.0), K)'
		gamma = [ones(K) for _ in 1:M]
		phi = ones(K, N[1]) / K

		model = new(K, M, V, N, C, copy(corp), topics, alpha, beta, gamma, phi)
		fixmodel!(model, check=false)

		model.newelbo = 0
		for d in 1:M
			model.phi = ones(K, N[d]) / K
			updateNewELBO!(model, d)
		end
		model.phi = ones(K, N[1]) / K
		updateELBO!(model)
		return model
	end
end

function Elogptheta(model::LDA)
	x = lgamma(sum(model.alpha)) - sum(lgamma.(model.alpha)) + dot(model.alpha - 1, model.Elogtheta)
	return x
end

function Elogpz(model::LDA, d::Int)
	counts = model.corp[d].counts
	x = dot(model.phi * counts, model.Elogtheta)
	return x
end

function Elogpw(model::LDA, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi .* log.(@boink model.beta[:,terms]) * counts)
	return x
end

function Elogqtheta(model::LDA, d::Int)
	x = -entropy(Dirichlet(model.gamma[d]))
	return x
end

function Elogqz(model::LDA, d::Int)
	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[:,n])) for (n, c) in enumerate(counts)])
	return x
end

function updateELBO!(model::LDA)
	model.elbo = model.newelbo
	model.newelbo = 0
	return model.elbo
end

function updateNewELBO!(model::LDA, d::Int)
	model.newelbo += (Elogptheta(model)
					+ Elogpz(model, d)
					+ Elogpw(model, d)
					- Elogqtheta(model, d)
					- Elogqz(model, d))
end

function updateAlpha!(model::LDA, niter::Integer, ntol::Real)
	"Interior-point Newton method with log-barrier and back-tracking line search."

	nu = model.K
	for _ in 1:niter
		rho = 1.0
		alphaGrad = [nu / model.alpha[i] + model.M * (digamma(sum(model.alpha)) - digamma(model.alpha[i])) for i in 1:model.K] + model.Elogthetasum
		alphaInvHessDiag = -1 ./ (model.M * trigamma.(model.alpha) + nu ./ model.alpha.^2)
		p = (alphaGrad - dot(alphaGrad, alphaInvHessDiag) / (1 / (model.M * trigamma(sum(model.alpha))) + sum(alphaInvHessDiag))) .* alphaInvHessDiag

		while minimum(model.alpha - rho * p) < 0
			rho *= 0.5
		end
		model.alpha -= rho * p

		if (norm(alphaGrad) < ntol) & (nu / model.K < ntol)
			break
		end
		nu *= 0.5
	end
	@bumper model.alpha
	model.Elogthetasum = zeros(model.K)
end

function updateBeta!(model::LDA)
	model.beta = model.newbeta ./ sum(model.newbeta, 2)
	model.newbeta = zeros(model.K, model.V)
end

function updateNewBeta!(model::LDA, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.newbeta[:,terms] += model.phi .* counts'
end

function updateGamma!(model::LDA, d::Int)
	counts = model.corp[d].counts
	@bumper model.gamma[d] = model.alpha + model.phi * counts
end

function updatePhi!(model::LDA, d::Int)
	terms = model.corp[d].terms
	model.phi = model.beta[:,terms] .* exp.(model.Elogtheta)
	model.phi ./= sum(model.phi, 1)
end

function updateElogtheta!(model::LDA, d::Int)
	model.Elogtheta = digamma.(model.gamma[d]) - digamma(sum(model.gamma[d]))
end

function updateElogthetasum!(model::LDA)
	model.Elogthetasum += model.Elogtheta
end

function train!(model::LDA; iter::Integer=150, tol::Real=1.0, niter::Integer=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, chkelbo::Integer=1)
	@assert all(.!isnegative.([tol, ntol, vtol]))
	@assert all(ispositive.([iter, niter, viter, chkelbo]))

	for k in 1:iter
		chk = (k % chkelbo == 0)
		for d in 1:model.M
			for _ in 1:viter
				oldgamma = model.gamma[d]
				updateElogtheta!(model, d)
				updatePhi!(model, d)
				updateGamma!(model, d)
				if norm(oldgamma - model.gamma[d]) < vtol
					break
				end
			end
			updateElogtheta!(model, d)
			updateElogthetasum!(model)
			updateNewBeta!(model, d)
			chk && updateNewELBO!(model, d)
		end
		updateAlpha!(model, niter, ntol)
		updateBeta!(model)
		if checkELBO!(model, k, chk, tol)
			break
		end
	end
	updatePhi!(model, 1)
	updateElogtheta!(model, 1)
	model.topics = [reverse(sortperm(vec(model.beta[i,:]))) for i in 1:model.K]
	nothing
end


text:
using LinearAlgebra

mutable struct Network

    n_exc::Int64
    n_neurons::Int64

    v::Array{Float64}
    u::Array{Float64}
    d::Array{Float64}
    a::Array{Float64}

    connection::Array{Int64,2} # WARNING read by columns !!!!
    connection_exc::Array{Int64,2} # WARNING read by columns !!!!
    connection_inh::Array{Int64,2} # WARNING read by columns !!!!
    s_exc::Array{Float64,2} # WARNING read by columns !!!!
    s_inh::Array{Float64,2} # WARNING read by columns !!!!
    sd::Array{Float64,2} # WARNING read by columns !!!!
    STDP::Array{Float64}
    ho::Array{Float64}

    exc::Array{Int64}
    inh::Array{Int64}

    I::Array{Float64}
    da::Float64
    arch::Array{Int64}
    #param::Dict

    function connect(n_neurons::Int64, n_exc::Int64, connectivity::Float64,
                     arch::Array{Int64}; inh_per::Float64 = 0.05 ) # parameter
        if arch == [0]
            connection = (rand(n_neurons,n_neurons) .< connectivity)
        else
            architr = [0; arch]
            connection = zeros(n_neurons,n_neurons)
            for i in 1:length(architr)-2
                si1 = sum(architr[1:i])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                ai2 = architr[i+1]
                ai3 = architr[i+2]
                connection[1+si1:si2,1+si2:si3] .= ones(ai2,ai3)
            end
            archinh = [0;arch[2:end-1]]
            for i in 1:length(archinh)-1
                sii1 = sum(archinh[1:i])
                sii2 = sum(archinh[1:i+1])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                a2 = arch[i+1]
                n_con_inh = floor(Int64,a2*inh_per)
                #mat_ie = ones(a2,a2) .- I(a2)
                mat_ie = zeros(a2,a2)
                for lin in 1:a2
                    ind = rand(1:a2,n_con_inh)
                    ind_f = filter(x -> x != i, ind)
                    mat_ie[i,ind_f] .+= 1
                end
                connection[n_exc+sii1+1:n_exc+sii2, si2+1:si3] .= mat_ie
                connection[si2+1:si3, n_exc+sii1+1:n_exc+sii2] .= I(a2)
            end
        end
        exc = collect(1:n_exc)
        inh = collect(n_exc+1:n_neurons)
        return connection, exc, inh
    end

    function Network(arch::Array{Int64})
        n_exc = sum(arch)
        n_neurons = n_exc+sum(arch[2:end-1])
        v = -65.0*ones(n_neurons) #parameter
        u = 0.2 * v # parameter
        d = [(i<=n_exc) ? 8.0 : 2.0 for i in 1:n_neurons] # parameter
        a = [(i<=n_exc) ? 0.02 : 0.1 for i in 1:n_neurons] # parameter
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, 0.1, arch)
        connection = transpose(connection)
        connection_exc = connection[:,exc]
        connection_inh = connection[:,inh]
        # maybe set the Exc-> Inh conn weights to high value
        s_exc = 1.0 .* ones(n_neurons, n_exc) # parameter
        s_inh = -0.5 .* ones(n_neurons, n_inh) # parameter
        s_exc .*= connection_exc
        s_inh .*= connection_inh
        sd = zeros(n_neurons,n_exc)
        STDP = zeros(n_neurons)
        ho = zeros(n_neurons)
        I = 13 .* (rand(n_neurons) .- 0.5) # parameters (noise in input current) random thalamic input
        da = 0.0
        arch = arch
        new(n_exc,n_neurons,v,u,d,a,connection,connection_exc,connection_inh,
            s_exc,s_inh,sd,STDP,ho,exc,inh,I,da,arch)
    end

    function Network(n_neurons::Int64, connectivity::Float64, n_exc::Int64)
        v = -65.0*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? 8.0 : 2.0 for i in 1:n_neurons]
        a = [(i<=n_exc) ? 0.02 : 0.1 for i in 1:n_neurons]
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, connectivity, [0])
        connection = transpose(connection)
        connection_exc = connection[:,exc]
        connection_inh = connection[:,inh]
        # maybe set the Exc-> Inh conn weights to high value
        s_exc = 1.0 .* ones(n_neurons, n_exc) # parameter
        s_inh = -0.5 .* ones(n_neurons, n_inh) # parameter (fixed inhibitory weight)
        s_exc .*= connection_exc
        s_inh .*= connection_inh
        sd = zeros(n_neurons,n_exc)
        STDP = zeros(n_neurons)
        ho = zeros(n_neurons)
        I = Float64[]
        da = 0.0
        arch = [0]
        new(n_exc,n_neurons,v,u,d,a,connection,connection_exc,connection_inh,
            s_exc,s_inh,sd,STDP,ho,exc,inh,I,da,arch)
    end
end

# %% functions

# direct spike of the neuron
function input!(net::Network, input_spikes::Array{Int64})
    net.v[input_spikes] .= 40.0
end

# superthreshold current processed by the neuron
function input!(net::Network, input_spikes::Array{Tuple{Int64,Float64}})
    net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    max_current = 20.0                  # parameter (max suprathreshold current)
    min_current = 0.0                  # parameter (min suprathreshold current)
    indices = [input_spikes[i][1] for i in eachindex(input_spikes)]
    intensity = [input_spikes[i][2] for i in eachindex(input_spikes)]
    net.I[indices] .+= (max_current-min_current) .* intensity .+ min_current
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
end

function spike!(net::Network)
    # spikes
    net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    spiked = findall(x->x>30,net.v .- net.ho) # parameter
    spiked_exc = filter(x->(x in net.exc), spiked)
    spiked_ho = filter(x-> x > net.arch[1], spiked)
    net.ho[spiked_ho] .+= 2.0 # parameter (theta)
    net.v[spiked] .= -65.0 # parameter (c)
    net.u[spiked] .+= net.d[spiked]
    # increment current
    for neuron in spiked
        if neuron in net.exc
            net.I .+= net.s_exc[:,neuron]
        elseif neuron in net.inh
            net.I .+= net.s_inh[:,neuron-net.n_exc]
        end
    end
    #time step (Izhikevic model)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.u .= @. net.u+net.a*(0.2*net.v-net.u)   # parameters (b)
    net.STDP[spiked] .= 0.1  # parameter (STDP_step)
    # LTP & LTD
    for neuron in spiked
        if neuron in net.exc
            net.sd[neuron,:] .+= 1.0 .* net.STDP[1:net.n_exc] #parameter (a_post) post synaptic spike
            net.sd[1:net.n_exc,neuron] .-= 1.5 .* net.STDP[1:net.exc] # parameter (a_pre) pre synaptic spike
        end
    end
    # time decay
    net.ho .*= 0.9999 # parameter (tau_theta)
    net.STDP .*= 0.95 # parameter (tau_STDP)
    net.da *= 0.995 # parameter (tau_da)
    spiked
end

function learn!(net::Network)
    #increase the synaptic weigths
    net.s_exc .+= (0.002+net.da) .* net.sd  # parameters
    net.s_exc .= clamp.(net.s_exc, 0, 4) # parameters (max synaptic weight)
    net.s_exc .*= net.connection_exc
    # time decay
    net.sd .*= 0.99 # parameter
end

# %% main functions

function step!(net::Network, train::Bool=false)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Int64}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

# %% visualization function

function weight_viz(net::Network)
    weights = net.s_exc
    filtered = filter(x-> x!=0.0, weights)
    display(histogram(filtered,
                      legend = false,
                      xlabel = "Synaptic strength",
                      ylabel = "Synapse number"))
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
function DataProvidersManager()
    return DataProvidersManager(())
end

function add_default_providers(obj::DataProvidersManager)
    return jcall(obj, "addDefaultProviders", void, ())
end

function add_filter(obj::DataProvidersManager, arg0::DataFilter)
    return jcall(obj, "addFilter", void, (DataFilter,), arg0)
end

function add_provider(obj::DataProvidersManager, arg0::DataProvider)
    return jcall(obj, "addProvider", void, (DataProvider,), arg0)
end

function apply_all_filters(obj::DataProvidersManager, arg0::NamedData)
    return jcall(obj, "applyAllFilters", NamedData, (NamedData,), arg0)
end

function clear_filters(obj::DataProvidersManager)
    return jcall(obj, "clearFilters", void, ())
end

function clear_loaded_data_names(obj::DataProvidersManager)
    return jcall(obj, "clearLoadedDataNames", void, ())
end

function clear_providers(obj::DataProvidersManager)
    return jcall(obj, "clearProviders", void, ())
end

function feed(obj::DataProvidersManager, arg0::JString, arg1::DataLoader)
    return jcall(obj, "feed", jboolean, (JString, DataLoader), arg0, arg1)
end

function get_loaded_data_names(obj::DataProvidersManager)
    return jcall(obj, "getLoadedDataNames", Set, ())
end

function get_providers(obj::DataProvidersManager)
    return jcall(obj, "getProviders", List, ())
end

function is_supported(obj::DataProvidersManager, arg0::DataProvider)
    return jcall(obj, "isSupported", jboolean, (DataProvider,), arg0)
end

function remove_provider(obj::DataProvidersManager, arg0::DataProvider)
    return jcall(obj, "removeProvider", DataProvider, (DataProvider,), arg0)
end


text:
module testJuliaPackage

greet() = print("Hello World!")

include("extra_file.jl")

# make this available for usage if the package is loaded in by
# using testJuliaPackage
# so only export functions which are of interest to users,
# keep the others internal
export my_f

end # module


text:
using HDF5
using HDF5Plugins
using H5Zlz4
using H5Zzstd
using H5Zbzip2
using HDF5Plugins.TestUtils
using Test

function test_zstd_filter()
    cd_values = Cuint[3] # aggression
    test_filter(H5Z_filter_zstd; cd_values = cd_values)
end

function test_bzip2_filter(data = ones(UInt8, 1024))
    cd_values = Cuint[8]
    test_filter(H5Z_filter_bzip2; cd_values = cd_values, data = data)
end

function test_lz4_filter()
    cd_values = Cuint[1024]
    test_filter(H5Z_filter_lz4; cd_values = cd_values)
end

@testset "HDF5Plugins.jl" begin
    for filter_func in (H5Z_filter_lz4, H5Z_filter_zstd, H5Z_filter_bzip2)
        @info filter_func
        # Decompressed size should be 1024, size of default data
        @test test_filter(filter_func)[2] == 0x400
        # Test random data of different lengths
        @test test_filter(filter_func; data = rand(UInt8, 1024))[2] == 0x400
        @test test_filter(filter_func; data = rand(UInt8, 3583))[2] == 3583
        @test test_filter(filter_func; data = rand(UInt8, 389))[2] == 389
        println()
    end

    # Test specific cd_values for filters
    @test test_lz4_filter()[2] == 0x400
    @test test_zstd_filter()[2] == 0x400
    @test test_bzip2_filter()[2] == 0x400

    # To Do, do actual tests on HDF5 files
end


text:
module BoxSymmetries
export BoxSym
export unit, inverse, ∘

using Random: AbstractRNG, Random
using ArgCheck
################################################################################
#### Perm
################################################################################
"""
Permutation
"""
struct Perm{N}
    perm::NTuple{N,Int}
    function Perm(perm::NTuple{N,Int}) where {N}
        @argcheck ispermutation(perm)
        return new{N}(perm)
    end
end
function Perm(args::Integer...)
    Perm(Int.(args))
end
function Base.show(io::IO, p::Perm)
    print(io, Perm, p.perm)
end
function ispermutation(perm::NTuple{N,Int}) where {N}
    ret = true
    for key in ntuple(identity, Val(N))
        ret = ret & (key in perm)
    end
    ret
end
ispermutation(perm::Perm) = true
function act_tuple(p::Perm{N}, t::NTuple{N}) where {N}
    map(p.perm) do i
        @inbounds t[i]
    end
end
function unit(::Type{Perm{N}}) where {N}
    Perm(ntuple(identity, Val(N)))
end
function Base.:(∘)(p1::Perm{N}, p2::Perm{N})::Perm{N} where {N}
    Perm(act_tuple(p1, p2.perm))
end
function inverse(p::Perm{N})::Perm{N} where {N}
    Perm(map(ntuple(identity, Val(N))) do i
        indexof(p.perm, i)
    end)
end
function indexof(t::NTuple, item)
    for i in eachindex(t)
        if t[i] === item
            return i
        end
    end
    error()
end
const PERM1 = [Perm(1)]
const PERM2 = [Perm(1,2), Perm(2,1)]
const PERM3 = [Perm(1,2,3),
               Perm(1,3,2), Perm(3,2,1), Perm(2,1,3),
               Perm(2,3,1), Perm(3,1,2),
]
Base.instances(::Type{Perm{1}}) = PERM1
Base.instances(::Type{Perm{2}}) = PERM2
Base.instances(::Type{Perm{3}}) = PERM3

################################################################################
#### BoxSym
################################################################################
struct BoxSym{N}
    # BoxSym is the semidirect product Boolⁿ ⋊ Sₙ
    # here axesperm ∈ Sₙ controls the induced action on the set of coordinate axes
    # flipsign controls for each axis, whether the sign is flipped. E.g. whether
    # eᵢ ↦ eⱼ or eᵢ ↦ -eⱼ for the standard basis
    axesperm::Perm
    flipsign::NTuple{N,Bool}
end
Base.ndims(o::BoxSym{N}) where {N} = N
Base.ndims(::Type{BoxSym{N}}) where {N} = N

function outaxes(o::BoxSym, axes::Tuple)
    @argcheck length(axes) == ndims(o)
    act_tuple(o.axesperm, axes)
end

isunit(g) = (g === unit(typeof(g)))

act_array!(out, g::BoxSym, x) = act_array_generic!(out, g, x)
act_array(g::BoxSym, x) = act_array_generic(g, x)

function act_array_generic!(out::AbstractArray{<:Any, N}, o::BoxSym{N}, x) where {N}
    @argcheck axes(out) == outaxes(o, axes(x))
    flipped_axes = map(flipaxis, axes(out), o.flipsign)
    @inbounds for I in CartesianIndices(x)
        inds = Tuple(I)
        inds = act_tuple(o.axesperm, inds)
        inds = map(getindex,flipped_axes, inds)
        Iout = CartesianIndex(inds)
        out[Iout] = x[I]
    end
    out
end
function act_array_generic(o::BoxSym, x)
    out = similar(x, outaxes(o,axes(x)))
    act_array_generic!(out, o, x)
end

function (o::BoxSym)(x)
    act_array(o,x)
end

function flipaxis(ax::AbstractRange, flip::Bool)
    if flip
        last(ax):-step(ax):first(ax)
    else
        first(ax):step(ax):last(ax)
    end
end

"""

    BoxSym(args::Integer...)

Create a box symmetry, according to args.
"""
function BoxSym(args::Integer...)
    perm = Perm(map(Int∘abs, args))
    flip = map(<(0), args)
    ret = BoxSym(perm, flip)
    ret
end
function astuple(o::BoxSym)
    (-1) .^ (o.flipsign) .* o.axesperm.perm
end
function Base.show(io::IO, o::BoxSym)
    print(io, "BoxSym", astuple(o))
end

function unit(::Type{BoxSym{N}})::BoxSym{N} where {N}
    BoxSym(ntuple(identity, Val(N))...)
end
function inverse(g::BoxSym{N})::BoxSym{N} where {N}
    n = g.flipsign
    h = g.axesperm
    h⁻¹ = inverse(h)
    n⁻¹ = n
    BoxSym{N}(
        h⁻¹,
        act_tuple(h⁻¹, n⁻¹),
   )
end
function Base.:(∘)(g1::BoxSym{N}, g2::BoxSym{N})::BoxSym{N} where {N}
    # We have
    # BoxSym = Boolⁿ ⋊ Sₙ
    # and we use the semidirect product formular
    n1 = g1.flipsign
    n2 = g2.flipsign
    h1 = g1.axesperm
    h2 = g2.axesperm
    h = h1∘h2
    n = map(⊻, n1, act_tuple(h1, n2))
    BoxSym{N}(h,n)
end

################################################################################
#### ALIASES
################################################################################
const BOXSYM_FROM_ALIAS2D = Dict(
    :unit      => BoxSym( 1, 2),
    :rot90     => BoxSym(-2, 1),
    :rot180    => BoxSym(-1,-2),
    :rot270    => BoxSym( 2,-1),
    :flipx     => BoxSym(-1, 2),
    :flipy     => BoxSym( 1,-2),
    :flipdiag  => BoxSym( 2, 1),
    :flipadiag => BoxSym(-2,-1),
)

BoxSym{2}(alias::Symbol) = BOXSYM_FROM_ALIAS2D[alias]

# TODO 3d aliases

################################################################################
#### instances
################################################################################

const SYMMETRIES1D = [BoxSym(1), BoxSym(-1)]
const SYMMETRIES2D = [
    BoxSym(1,2), BoxSym(1,-2), BoxSym(-1,2), BoxSym(-1,-2),
    BoxSym(2,1), BoxSym(2,-1), BoxSym(-2,1), BoxSym(-2,-1),
]
const SYMMETRIES3D = [
    BoxSym(1,2,3), BoxSym(1,2,-3), BoxSym(1,-2,3), BoxSym(1,-2,-3), BoxSym(-1,2,3), BoxSym(-1,2,-3), BoxSym(-1,-2,3), BoxSym(-1,-2,-3),
    BoxSym(1,3,2), BoxSym(1,3,-2), BoxSym(1,-3,2), BoxSym(1,-3,-2), BoxSym(-1,3,2), BoxSym(-1,3,-2), BoxSym(-1,-3,2), BoxSym(-1,-3,-2),
    BoxSym(3,2,1), BoxSym(3,2,-1), BoxSym(3,-2,1), BoxSym(3,-2,-1), BoxSym(-3,2,1), BoxSym(-3,2,-1), BoxSym(-3,-2,1), BoxSym(-3,-2,-1),
    BoxSym(2,1,3), BoxSym(2,1,-3), BoxSym(2,-1,3), BoxSym(2,-1,-3), BoxSym(-2,1,3), BoxSym(-2,1,-3), BoxSym(-2,-1,3), BoxSym(-2,-1,-3),
    BoxSym(2,3,1), BoxSym(2,3,-1), BoxSym(2,-3,1), BoxSym(2,-3,-1), BoxSym(-2,3,1), BoxSym(-2,3,-1), BoxSym(-2,-3,1), BoxSym(-2,-3,-1),
    BoxSym(3,1,2), BoxSym(3,1,-2), BoxSym(3,-1,2), BoxSym(3,-1,-2), BoxSym(-3,1,2), BoxSym(-3,1,-2), BoxSym(-3,-1,2), BoxSym(-3,-1,-2),
]

Base.instances(::Type{BoxSym{1}}) = SYMMETRIES1D
Base.instances(::Type{BoxSym{2}}) = SYMMETRIES2D
Base.instances(::Type{BoxSym{3}}) = SYMMETRIES3D

Random.rand(rng::AbstractRNG, G::Type{<:BoxSym}) = rand(rng, instances(G))

################################################################################
#### Precompile
################################################################################
let
    BoxSym(1)(fill(1.0, 1))
    BoxSym(1,-2)(fill(1.0, 1,1))
    BoxSym(1,-2,3)(fill(1.0, 1,1,1))
end

# TODO lazy iterator of all symmetries of given dimension
# TODO isrotation for checking if orientation is preserved
end #module


text:
using HorizonSideRobots

function goToSide!(r::Robot, side::HorizonSide; save::Bool=false)::Union{Tuple{Int, HorizonSide}, Nothing}
    if save
        shags = 0
        while !isborder(r, side)
            move!(r, side)
            shags+=1
        end
        return (shags, reverseSide(side))
    else
        while !isborder(r, side)
            move!(r, side)
        end
    end
end

function goToSideWithShag!(r::Robot, side::HorizonSide, shags::Int; save::Bool=false, checkBorder::Bool = false, mark::Bool = false)::Union{Tuple{Int, HorizonSide}, Nothing}
    if save
        shag = 0
        for _ = 1:shags
            if checkBorder && isborder(r, side) break end
            move!(r, side)
            if mark putmarker!(r) end
            shag+=1
        end
        return (shag, reverseSide(side))
    else
        for _ = 1:shags
            if checkBorder && isborder(r, side) break end
            move!(r, side)
        end
    end
end

function reverseSide(side::HorizonSide)::HorizonSide
    return HorizonSide(
        (Integer(side)+2) % 4
    )
end

function perSide(side::HorizonSide)::HorizonSide
    return HorizonSide(
        (Integer(side)+1) % 4
    )
end

function goToStart!(r::Robot, listOfShags)::Nothing
    reverse!(listOfShags)
    for tup = listOfShags
        for _ = 1:tup[1]
            move!(r, tup[2])
        end
    end
end

function goToWestSud!(r::Robot)::Nothing
    goToSide!(r, Sud)
    goToSide!(r, West)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
pages = [
        "DiffEqOperators.jl: Linear operators for Scientific Machine Learning" => "index.md",
        "Operator Tutorials" => [
            "operator_tutorials/kdv.md"
        ],
        "Operators" => [
            "operators/operator_overview.md",
            "operators/derivative_operators.md",
            "operators/vector_calculus_operators.md",
            "operators/vector_jacobian_product.md",
            "operators/jacobian_vector_product.md",
            "operators/matrix_free_operators.md"
        ],
        "Nonlinear Derivatives" => [
            "nonlinear_derivatives/nonlinear_diffusion.md"
        ]
     ]