text:
################################################################################
#                                  STATUS QUERIES
################################################################################
"""
    JuMP.termination_status(model::InfiniteModel)

Extend [`termination_status`](@ref JuMP.termination_status(::JuMP.Model)) to
return the `MOI.TerminationStatus` in accordance with the optimizer model.

**Example**
```julia-repl
julia> termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4
```
"""
function JuMP.termination_status(model::InfiniteModel)
    return JuMP.termination_status(optimizer_model(model))
end

"""
    JuMP.raw_status(model::InfiniteModel)

Extend [`raw_status`](@ref JuMP.raw_status(::JuMP.Model)) to return the status
reported by the solver in accordance with the optimizer model.

**Example**
```julia-repl
julia> raw_status(model) # Ipopt
"Solve_Succeeded"
```
"""
function JuMP.raw_status(model::InfiniteModel)
    return JuMP.raw_status(optimizer_model(model))
end

"""
     JuMP.primal_status(model::InfiniteModel; [result::Int = 1])

Extend [`primal_status`](@ref JuMP.primal_status(::JuMP.Model)) to return the
`MOI.PrimalStatus` reported in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
```
"""
function  JuMP.primal_status(model::InfiniteModel; result::Int = 1)
    return JuMP.primal_status(optimizer_model(model); result = result)
end

"""
    JuMP.dual_status(model::InfiniteModel; [result::Int = 1])

Extend [`dual_status`](@ref JuMP.dual_status(::JuMP.Model)) to return the
`MOI.DualStatus` reported in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
```
"""
function JuMP.dual_status(model::InfiniteModel; result::Int = 1)
    return JuMP.dual_status(optimizer_model(model); result = result)
end

"""
    JuMP.solve_time(model::InfiniteModel)

Extend [`solve_time`](@ref JuMP.solve_time(::JuMP.Model)) to return the
time used by the solver to terminate reported in accordance with the optimizer
model. This will error if not supported by the solver.

**Example**
```julia-repl
julia> solve_time(model)
0.004999876022338867
```
"""
function JuMP.solve_time(model::InfiniteModel)
    return JuMP.solve_time(optimizer_model(model))
end

"""
    JuMP.has_values(model::InfiniteModel; [result::Int = 1])::Bool

Extend [`has_values`](@ref JuMP.has_values(::JuMP.Model)) to return a `Bool`
whether variable values are available in accordance with the optimizer model and
the result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> has_values(model)
true
```
"""
function JuMP.has_values(model::InfiniteModel; result::Int = 1)::Bool
    return JuMP.primal_status(model; result = result) != MOI.NO_SOLUTION
end

"""
    JuMP.has_duals(model::InfiniteModel; [result::Int = 1])::Bool

Extend [`has_duals`](@ref JuMP.has_duals(::JuMP.Model)) to return a `Bool`
whether constraint duals are available in accordance with the optimizer model and
the result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> has_duals(model)
true
```
"""
function JuMP.has_duals(model::InfiniteModel; result::Int = 1)::Bool
    return JuMP.dual_status(model; result = result) != MOI.NO_SOLUTION
end

"""
    simplex_iterations(model::InfiniteModel)

Gets the cumulative number of simplex iterations during the most-recent optimization.

Solvers must implement `MOI.SimplexIterations()` to use this function.
"""
function JuMP.simplex_iterations(model::InfiniteModel)
    return JuMP.simplex_iterations(optimizer_model(model))
end

"""
    barrier_iterations(model::InfiniteModel)

Gets the cumulative number of barrier iterations during the most-recent optimization.

Solvers must implement `MOI.BarrierIterations()` to use this function.
"""
function JuMP.barrier_iterations(model::InfiniteModel)
    return JuMP.barrier_iterations(optimizer_model(model))
end

"""
    node_count(model::InfiniteModel)

Gets the total number of branch-and-bound nodes explored during the most recent
optimization in a Mixed Integer Program.

Solvers must implement `MOI.NodeCount()` to use this function.
"""
function JuMP.node_count(model::InfiniteModel)
    return JuMP.node_count(optimizer_model(model))
end

################################################################################
#                              OBJECTIVE QUERIES
################################################################################
"""
    JuMP.objective_bound(model::InfiniteModel)::Float64

Extend [`objective_bound`](@ref JuMP.objective_bound(::JuMP.Model)) to return
the objective bound in accordance with the optimizer model.

**Example**
```julia-repl
julia> objective_bound(model)
42.0
```
"""
function JuMP.objective_bound(model::InfiniteModel)::Float64
    return JuMP.objective_bound(optimizer_model(model))
end

"""
    JuMP.objective_value(model::InfiniteModel; [result::Int = 1])::Float64

Extend [`objective_value`](@ref JuMP.objective_value(::JuMP.Model)) to return
the objective value in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> objective_value(model)
42.0
```
"""
function JuMP.objective_value(model::InfiniteModel; result::Int = 1)::Float64
    return JuMP.objective_value(optimizer_model(model); result = result)
end

"""
    JuMP.dual_objective_value(model::InfiniteModel; [result::Int = 1])::Float64

Extend [`dual_objective_value`](@ref JuMP.dual_objective_value(::JuMP.Model)) to
return the dual objective value in accordance with the optimizer model and the
result index `result` of the most recent solution obtained. Errors
if the solver does not support this.

**Example**
```julia-repl
julia> dual_objective_value(model)
42.0
```
"""
function JuMP.dual_objective_value(model::InfiniteModel; result::Int = 1)::Float64
    return JuMP.dual_objective_value(optimizer_model(model); result = result)
end

################################################################################
#                                 VALUE QUERIES
################################################################################
"""
    map_value([ref/expr], key::Val{ext_key_name}, result::Int; kwargs...)

Map the value(s) of `ref` to its counterpart in the optimizer model type that is
distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`, `optimizer_model_expression`, and/or
`optimizer_model_constraint`. Such as is the case with reformuations that do not
have a direct mapping between variables and/or constraints in the original
infinite form. Otherwise, `optimizer_model_variable`,
`optimizer_model_expression`, and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are passed on these functions. Here
`result` is the result index used in `value`.
"""
function map_value end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_value(vref::GeneralVariableRef, key, result::Int; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.value(v; result = result), opt_vref)
    else
        return JuMP.value(opt_vref; result = result)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_value(cref::InfOptConstraintRef, key, result::Int; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.value(c; result = result), opt_cref)
    else
        return JuMP.value(opt_cref; result = result)
    end
end

# Default method that depends on optimizer_model_expression --> making extensions easier
function map_value(expr::JuMP.AbstractJuMPScalar, key, result::Int; kwargs...)
    opt_expr = optimizer_model_expression(expr, key; kwargs...)
    if opt_expr isa AbstractArray
        return map(e -> JuMP.value(e; result = result), opt_expr)
    else
        return JuMP.value(opt_expr; result = result)
    end
end

## Define dispatch methods to collect value of parameters
# InfiniteParameter
function _get_value(pref, ::Type{<:InfiniteParameterIndex}, result; kwargs...)
    label = get(kwargs, :label, PublicLabel)
    return supports(pref, label = label)
end

# FiniteParameter
function _get_value(pref, ::Type{<:FiniteParameterIndex}, result; kwargs...)
    return parameter_value(pref)
end

# Others
function _get_value(vref, index_type, result; kwargs...)
    return map_value(vref, Val(optimizer_model_key(JuMP.owner_model(vref))),
                     result; kwargs...)
end

"""
    JuMP.value(vref::GeneralVariableRef; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Extend [`JuMP.value`](@ref JuMP.value(::JuMP.VariableRef)) to return the value(s)
of `vref` in accordance with its reformulation variable(s) stored in the optimizer
model and the result index `result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check
if a result exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

To provide context for the
results it may be helpful to also query the variable's `parameter_refs` and
`supports` which will have a one-to-one correspondence with the value(s).
It may also be helpful to query via [`optimizer_model_variable`](@ref) to
retrieve the variables(s) that these values are based on. These functions should
all be called with the same keyword arugments for consistency.

For extensions, this only works if
[`optimizer_model_variable`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> value(z)
42.0
```
"""
function JuMP.value(vref::GeneralVariableRef; result::Int = 1, kwargs...)
    return _get_value(vref, _index_type(vref), result; kwargs...)
end

################################################################################
#                                 REDUCED COST
################################################################################

"""
    map_reduced_cost([ref], key::Val{ext_key_name}, result::Int; kwargs...)

Map the reduced cost(s) of `ref` to its counterpart in the optimizer model type that is
distininguished by its extension key `key` as type `Val{ext_key_name}`. Here `ref` needs refer
to methods for variable references.  This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`. Such as is the case with reformuations that do not
have a direct mapping between variables in the original
infinite form. Otherwise, `optimizer_model_variable`, is used to make
these mappings by default where `kwargs` are passed on these functions. Here
`result` is the result index used in `value`.
"""
function map_reduced_cost end

function map_reduced_cost(vref::GeneralVariableRef, key; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.reduced_cost(v), opt_vref)
    else
        return JuMP.reduced_cost(opt_vref)
    end
end

"""
    JuMP.reduced_cost(vref::GeneralVariableRef)
Extending [`JuMP.reduced_cost`](@ref JuMP.reduced_cost(::JuMP.VariableRef)). This returns
the reduced cost of some infinite variable. Returns a scalar value.

**Example**
```
julia> reduced_cost(x)
12.81

```
"""
function JuMP.reduced_cost(vref::GeneralVariableRef; kwargs...)
    return map_reduced_cost(vref, Val(optimizer_model_key(JuMP.owner_model(vref))); kwargs...)
end

"""
    JuMP.value(cref::InfOptConstraintRef; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Extend [`JuMP.value`](@ref JuMP.value(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the value(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model and the result index `result` of the most recent
solution obtained. Use [`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel))
to check if a result exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

To provide context for
the results it may be helpful to also query the constraint's `parameter_refs`
and `supports` which will have a one-to-one correspondence with the value(s).
It may also be helpful to query via [`optimizer_model_constraint`](@ref) to
retrieve the constraint(s) that these values are based on. By default, only the
values corresponding to public supports are returned. These functions should
all be called with the same keyword arugments for consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> value(c1)
4-element Array{Float64,1}:
 -0.0
 20.9
 20.9
 20.9
```
"""
function JuMP.value(cref::InfOptConstraintRef; result::Int = 1,
                    kwargs...)
    return map_value(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                     result; kwargs...)
end

"""
    JuMP.value(expr::JuMP.AbstractJuMPScalar; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Return the value(s) of `expr` in accordance with the optimized variable values
the result index `result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check if a result
exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
expressions are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the expression has multiple
infinite parameter dependencies.

To provide context for the results it may be
helpful to also query the expression's `parameter_refs` and `supports` which
will have a one-to-one correspondence with the value(s). It may also be helpful
to query via [`optimizer_model_expression`](@ref) to retrieve the expression(s)
that these values are based on. These should use the same keyword arguments for
consistency.

For extensions, this only works if [`optimizer_model_expression`](@ref) has been
extended correctly and/or [`map_value`](@ref) has been extended for expressions.

**Example**
```julia-repl
julia> value(my_finite_expr)
23.34

julia> value(my_infinite_expr)
4-element Array{Float64,1}:
 -0.0
 20.9
 20.9
 20.9
```
"""
function JuMP.value(expr::Union{JuMP.GenericAffExpr{C, V}, JuMP.GenericQuadExpr{C, V}};
    result::Int = 1,
    kwargs...
    ) where {C, V <: GeneralVariableRef}
    # get the model
    model = _model_from_expr(expr)
    # if no model then the expression only contains a constant
    if model === nothing
        return JuMP.constant(expr)
    # otherwise let's call map_value
    else
        key = optimizer_model_key(model)
        return map_value(expr, Val(key), result; kwargs...)
    end
end

################################################################################
#                             OPTIMIZER INDEX QUERIES
################################################################################
"""
    map_optimizer_index(ref, key::Val{ext_key_name}; kwargs...)

Map the `MathOptInterface` index(es) of `ref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are passed on as well.
"""
function map_optimizer_index end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_optimizer_index(vref::GeneralVariableRef, key; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.optimizer_index(v), opt_vref)
    else
        return JuMP.optimizer_index(opt_vref)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_optimizer_index(cref::InfOptConstraintRef, key; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.optimizer_index(c), opt_cref)
    else
        return JuMP.optimizer_index(opt_cref)
    end
end

"""
    JuMP.optimizer_index(vref::GeneralVariableRef;
                         [label::Type{<:AbstractSupportLabel} = PublicLabel,
                         ndarray::Bool = false, kwargs...])

Extend [`JuMP.optimizer_index`](@ref JuMP.optimizer_index(::JuMP.VariableRef)) to
return the `MathOptInterface` index(es) of `vref` in accordance with its
reformulation variable(s) stored in the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the optimizer indices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the indices of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_variable`](@ref) to
retrieve the variables(s) that these indices are based on. These should use the
same keyword arguments for consistency.

For extensions, this
only works if [`optimizer_model_variable`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> optimizer_index(x)
4-element Array{MathOptInterface.VariableIndex,1}:
 MathOptInterface.VariableIndex(2)
 MathOptInterface.VariableIndex(3)
 MathOptInterface.VariableIndex(4)
 MathOptInterface.VariableIndex(5)
```
"""
function JuMP.optimizer_index(vref::GeneralVariableRef; kwargs...)
    return map_optimizer_index(vref, Val(optimizer_model_key(JuMP.owner_model(vref)));
                               kwargs...)
end

"""
    JuMP.optimizer_index(cref::InfOptConstraintRef;
                         [label::Type{<:AbstractSupportLabel} = PublicLabel,
                         ndarray::Bool = false, kwargs...])

Extend [`JuMP.optimizer_index`](@ref JuMP.optimizer_index(::JuMP.ConstraintRef{JuMP.Model}))
to return the `MathOptInterface` index(es) of `cref` in accordance with its
reformulation constraints(s) stored in the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the optimizer indices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the indices of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_constraint`](@ref) to
retrieve the constraints(s) that these indices are based on. The same keyword
arguments should be used for consistency.

For extensions, this
only works if [`optimizer_model_constraint`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> optimizer_index(c1)
4-element Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},1}:
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(1)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(2)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(3)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(4)
```
"""
function JuMP.optimizer_index(cref::InfOptConstraintRef; kwargs...)
    return map_optimizer_index(cref, Val(optimizer_model_key(JuMP.owner_model(cref)));
                               kwargs...)
end

################################################################################
#                                  DUAL QUERIES
################################################################################
"""
    map_dual(cref::InfOptConstraintRef, key::Val{ext_key_name}, result::Int;
             kwargs...)

Map the dual(s) of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are also pass on to. Here `result` is
the result index that is used in `dual`.
"""
function map_dual end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_dual(cref::InfOptConstraintRef, key, result::Int; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.dual(c; result = result), opt_cref)
    else
        return JuMP.dual(opt_cref; result = result)
    end
end

"""
    JuMP.dual(cref::InfOptConstraintRef; [result::Int = 1,
              label::Type{<:AbstractSupportLabel} = PublicLabel,
              ndarray::Bool = false, kwargs...])

Extend [`JuMP.dual`](@ref JuMP.dual(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the dual(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model and the result index `result` of the most recent
solution obtained. Use [`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel))
to check if a result exists before asking for duals.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the duals associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the duals of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to
query via [`optimizer_model_constraint`](@ref) to retrieve the constraint(s)
that these duals are based on. Calling `parameter_refs` and `supports` may also
be insightful. Be sure to use the same keyword arguments for consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_dual`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> dual(c1)
4-element Array{Float64,1}:
 -42.0
 -42.0
 32.3
 0.0
```
"""
function JuMP.dual(cref::InfOptConstraintRef; result::Int = 1, kwargs...)
    return map_dual(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                    result; kwargs...)
end

# Error redriect for variable call
function JuMP.dual(vref::GeneralVariableRef; kwargs...)
    return JuMP.dual(JuMP.VariableRef(JuMP.Model(), MOI.VariableIndex(1)))
end

################################################################################
#                              SHADOW PRICE QUERIES
################################################################################
"""
    map_shadow_price(cref::InfOptConstraintRef, key::Val{ext_key_name};
                     kwargs...)

Map the shadow price(s) of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are passed on to.
"""
function map_shadow_price end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_shadow_price(cref::InfOptConstraintRef, key; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.shadow_price(c), opt_cref)
    else
        return JuMP.shadow_price(opt_cref)
    end
end

"""
    JuMP.shadow_price(cref::InfOptConstraintRef;
                      [label::Type{<:AbstractSupportLabel} = PublicLabel,
                      ndarray::Bool = false, kwargs...])

Extend [`JuMP.shadow_price`](@ref JuMP.shadow_price(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the shadow price(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model. Use [`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel))
to check if a result exists before asking for duals.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the shadow prices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the prices of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be
helpful to query via [`optimizer_model_constraint`](@ref) to retrieve the
constraint(s) that these shadow prices are based on. Calling `parameter_refs` and
`supports` may also be insightful. Be sure to use the same keyword arguments for
consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_shadow_price`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> shadow_price(c1)
4-element Array{Float64,1}:
 42.0
 42.0
 -32.3
 -0.0
```
"""
function JuMP.shadow_price(cref::InfOptConstraintRef; kwargs...)
    return map_shadow_price(cref, Val(optimizer_model_key(JuMP.owner_model(cref)));
                            kwargs...)
end

################################################################################
#                               PERTURBATION QUERIES
################################################################################
"""
    map_lp_rhs_perturbation_range(cref::InfOptConstraintRef,
                                  key::Val{ext_key_name}, toler::Float64;
                                  kwargs...)

Map the RHS perturbation range of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `cref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_constraint` is used to make these mappings by default where `kwargs`
are passed on to. Here `toler` corresponds to the `feasibility_tolerance` used by
`lp_rhs_perturbation_range`.
"""
function map_lp_rhs_perturbation_range end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_lp_rhs_perturbation_range(cref::InfOptConstraintRef, key,
                                       toler::Float64; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.lp_rhs_perturbation_range(c; feasibility_tolerance = toler),
                   opt_cref)
    else
        return JuMP.lp_rhs_perturbation_range(opt_cref; feasibility_tolerance = toler)
    end
end

"""
    JuMP.lp_rhs_perturbation_range(cref::InfOptConstraintRef;
                                   [feasibility_tolerance::Float64 = 1e-8,
                                   label::Type{<:AbstractSupportLabel} = PublicLabel,
                                   ndarray::Bool = false, kwargs...])

Extend [`JuMP.lp_rhs_perturbation_range`](@ref JuMP.lp_rhs_perturbation_range(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the range(s) of the RHS of `cref` for which the shadow price(s) are valid
in accordance with its reformulation constraint(s)
stored in the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the ranges associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the ranges of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to
query via [`optimizer_model_constraint`](@ref) to retrieve the constraint(s)
that these ranges are based on. Calling `parameter_refs` and `supports` may also
be insightful. Be sure to use the same keyword arguments for consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_lp_rhs_perturbation_range`](@ref) has been implemented.

**Example**
```julia-repl
julia> lp_rhs_perturbation_range(c1)
4-element Array{Tuple{Float64,Float64},1}:
 (-42.0, Inf)
 (-Inf, 42.0)
 (-Inf, 42.0)
 (-Inf, 42.0)
```
"""
function JuMP.lp_rhs_perturbation_range(cref::InfOptConstraintRef;
    feasibility_tolerance::Float64 = 1e-8,
    kwargs...
    )
    return map_lp_rhs_perturbation_range(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                                         feasibility_tolerance; kwargs...)
end

"""
    map_lp_objective_perturbation_range(vref::DecisionVariableRef,
                                        key::Val{ext_key_name}, toler::Float64;
                                        kwargs...)

Map the reduced cost range(s) of `vref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `vref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` is used to make these mappings by default where `kwargs`
are passed on as well. Here `toler` corresponds to the `optimality_tolerance` used by
`lp_objective_perturbation_range`.
"""
function map_lp_objective_perturbation_range end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_lp_objective_perturbation_range(vref::GeneralVariableRef, key,
                                             toler::Float64; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.lp_objective_perturbation_range(v; optimality_tolerance = toler),
                   opt_vref)
    else
        return JuMP.lp_objective_perturbation_range(opt_vref; optimality_tolerance = toler)
    end
end

"""
    JuMP.lp_objective_perturbation_range(vref::InfOptVariableRef;
                                         [optimality_tolerance::Float64 = 1e-8,
                                         label::Type{<:AbstractSupportLabel} = PublicLabel,
                                         ndarray::Bool = false, kwargs...])

Extend [`JuMP.lp_objective_perturbation_range`](@ref JuMP.lp_objective_perturbation_range(::JuMP.VariableRef))
to return the range(s) that the reduced cost(s) of `vref` remain valid
in accordance with its reformulation variables(s)
stored in the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the ranges associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the ranges of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_variable`](@ref) to retrieve
the variable(s) that these ranges are based on. Calling `parameter_refs` and
`supports` may also be insightful. Be sure to use the same keyword arguments for
conistency.

For extensions, this only
works if [`optimizer_model_variable`](@ref) has been extended correctly and/or
if [`map_lp_objective_perturbation_range`](@ref) has been implemented.

**Example**
```julia-repl
julia> lp_objective_perturbation_range(z)
(-2.0, Inf)
```
"""
function JuMP.lp_objective_perturbation_range(vref::GeneralVariableRef;
    optimality_tolerance::Float64 = 1e-8,
    kwargs...
    )
    return map_lp_objective_perturbation_range(vref, Val(optimizer_model_key(JuMP.owner_model(vref))),
                                               optimality_tolerance; kwargs...)
end


text:
include("euler/euler.jl")
using .Game: dice_combinations
using .Math: fast_max
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(a::Tuple{Int, Int}, b::Tuple{Int, Int})::Float64
    limit = fast_max(a[1] * a[2], b[1] * b[2])
    looses = zeros(Int, limit)
    result = 0
    for i ∈ b[1]:b[1] * b[2] looses[i] = dice_combinations(i, b[1], b[2]) + looses[i - 1] end
    for i ∈ a[1]:a[1] * a[2] result += dice_combinations(i, a[1], a[2]) * looses[i - 1] end
    return round(result / a[2] ^ a[1] / b[2] ^ b[1], digits=7)
end

compute((9, 4), (6, 6))

@benchmark compute((9, 4), (6, 6))

ext: jl

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
function betterbubblesort!(A)
    n = length(A)
    v=0
    j=1
    for i in 1:n-1
        for j in 1:n-i
            if A[j]>A[j+1]
                A[j], A[j+1] = A[j+1], A[j]
                v=j+1
            end
        end
        j=v
    end
    return A
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function vactoair(wave_vac::AbstractFloat)
    if wave_vac >= 2000
        sigma2 = (1e4/wave_vac)^2  # Convert to wavenumber squared
        # Computer conversion factor.
        fact = 1 + 5.792105e-2/(238.0185 - sigma2) +
            1.67917e-3/(57.362 - sigma2)
        return wave_vac/fact # Convert wavelength
    else
        return wave_vac
    end
end

"""
    vactoair(wave_vacuum) -> wave_air

### Purpose ###

Converts vacuum wavelengths to air wavelengths.

### Explanation ###

Corrects for the index of refraction of air under standard conditions.
Wavelength values below \$2000 Å\$ will not be altered.  Uses relation of Ciddor
(1996).

### Arguments ###

* `wave_vacuum`: vacuum wavelength in angstroms.  Wavelengths are corrected for the index of
  refraction of air under standard conditions.  Wavelength values below \$2000 Å\$ will
  *not* be altered, take care within \$[1 Å, 2000 Å]\$.

### Output ###

Air wavelength in angstroms.

### Method ###

Uses relation of Ciddor (1996), Applied Optics 35, 1566
(http://adsabs.harvard.edu/abs/1996ApOpt..35.1566C).

### Example ###

If the vacuum wavelength is `w = 2000`, then `vactoair(w)` yields an air
wavelength of `1999.353`.

```jldoctest
julia> using AstroLib

julia> vactoair(2000)
1999.3526230448367
```

### Notes ###

`airtovac` converts air wavelengths to vacuum wavelengths.

Code of this function is based on IDL Astronomy User's Library.
"""
vactoair(wave_vac::Real) = vactoair(float(wave_vac))


text:
# Separable sum, using arrays of arrays as variables

immutable SeparableSum{S <: AbstractArray} <: ProximableFunction
	fs::S
end

function (f::SeparableSum{S}){S <: AbstractArray}(x::AbstractArray)
	sum = 0.0
  for k in eachindex(f.fs)
		sum += f.fs[k](x[k])
	end
	return sum
end

function prox!{S <: AbstractArray}(f::SeparableSum{S}, x::AbstractArray, y::AbstractArray, gamma::Real=1.0)
  sum = 0.0
  for k in eachindex(f.fs)
	  sum += prox!(f.fs[k], x[k], y[k], gamma)
  end
  return sum
end

function prox(f::SeparableSum, x::AbstractArray, gamma::Real=1.0)
	y = map(similar, x)
	fy = prox!(f, x, y, gamma)
	return y, fy
end

fun_name(f::SeparableSum) = "separable sum"
function fun_dom(f::SeparableSum)
	# s = ""
	# for k in eachindex(f.fs)
	# 	s = string(s, fun_dom(f.fs[k]), " × ")
	# end
	# return s
	return "n/a" # for now
end
fun_expr(f::SeparableSum) = "(x₁,…,xᵣ) ↦ f[1](x₁) + … + f[r](xᵣ)"
function fun_params(f::SeparableSum)
	# s = "r = $(size(f.fs)), "
	# for k in eachindex(f.fs)
	# 	s = string(s, "f[$(k)] = $(typeof(f.fs[k])), ")
	# end
	# return s
	return "n/a" # for now
end


text:
module FiniteDifferences

using ChainRulesCore
using LinearAlgebra
using Printf
using Random
using Richardson
using SparseArrays
using StaticArrays

export to_vec, grad, jacobian, jvp, j′vp

include("deprecated.jl")
include("methods.jl")
include("numerics.jl")
include("to_vec.jl")
include("grad.jl")

end


text:
# Autogenerated wrapper script for VerizonEctoken_jll for x86_64-linux-gnu
export ectoken3

using OpenSSL_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `ectoken3`
const ectoken3_splitpath = ["bin", "ectoken3"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ectoken3_path = ""

# ectoken3-specific global declaration
function ectoken3(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ectoken3_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"VerizonEctoken")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (OpenSSL_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (OpenSSL_jll.LIBPATH_list,))

    global ectoken3_path = normpath(joinpath(artifact_dir, ectoken3_splitpath...))

    push!(PATH_list, dirname(ectoken3_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Pkg.Artifacts
const artifact_toml = find_artifacts_toml(@__DIR__)

# Kernels to be downloaded lazily
const SPICE_url_base = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels"
const all_kernels = Dict(
    "meta"    => ["pck/pck00010.tpc",
                  "pck/gm_de431.tpc",
                  "lsk/latest_leapseconds.tls"],  # Required for all other kernels
    "default" => ["spk/planets/de430.bsp"],
    "mars"    => ["spk/satellites/mar097.bsp"],
    "jupiter" => ["spk/satellites/jup310.bsp"],
    "saturn"  => ["spk/satellites/sat427.bsp"],
    "uranus"  => ["spk/satellites/ura111.bsp"],
    "neptune" => ["spk/satellites/nep095.bsp"],
    "pluto"   => ["spk/satellites/plu055.bsp"]
)
for (system, kernels) in all_kernels
    artifact_name = "$(system)_kernels"
    artifact = artifact_hash(artifact_name, artifact_toml)
    if isnothing(artifact)
        # Only download the kernel if it doesn't exist in Artifacts.toml. Otherwise, we want to
        # download it lazily in SpiceUtils.
        download_info = [(join((SPICE_url_base, kernel), "/"), "") for kernel in kernels]
        artifact = create_artifact() do artifact_dir
            for (download_URL, _) in download_info
                @info "Downloading $(basename(download_URL))..."
                download(download_URL, joinpath(artifact_dir, basename(download_URL)))
            end
        end
        if Sys.islinux()
            bind_artifact!(artifact_toml, artifact_name, artifact; download_info, lazy=true, force=true)
        else
            @warn "Could not bind artifact, need to run on Linux so the hashes are correct!"
        end
    end
end


text:
@testset "The Extended Decorator" begin
    M = Euclidean(2)
    x = RnPoint([1.,0.])
    ξ = RnTVector([0.,1.])
    xE = MPointE(x)
    ξE = TVectorE(ξ,xE)
    # check stripping for double decorate - check that we never decorate double
    @test MPointE(xE).base == x
    @test ξE.base == x
    @test getBasePoint(ξE) == xE #access methods promote
    ξEE = TVectorE(ξE,x)
    @test ξEE.base == x
    @test ξEE.vector == ξ
    ξEE2 = TVectorE(ξE,xE)
    @test ξEE.base == x
    @test ξEE.vector == ξ
    #
    promote_rule(MPointE{RnPoint{Float64}}, RnPoint{Float64}) == MPointE{RnPoint{Float64}}
    # convert MPoints
    @test convert(MPointE{RnPoint{Float64}},x) == xE
    @test_throws ErrorException convert(MPointE{SnPoint{Float64}},x)
    @test convert(MPointE{RnPoint{Float64}},xE) == xE
    @test convert(RnPoint{Float64},xE) == x
    @test_throws ErrorException convert(SnPoint{Float64},xE)
    @test convert(RnPoint{Float64},x) == x
    # convert TVecs
    @test_throws ErrorException convert(TVectorE{SnTVector{Float64},SnPoint{Float64}},ξE)
    @test convert(RnTVector{Float64}, ξE) == ξ
    @test_throws ErrorException convert(SnTVector{Float64},ξE)
    @test convert(TVectorE{RnTVector{Float64}, RnPoint{Float64}},ξE) == ξE
    @test convert(RnTVector{Float64},ξ) == ξ
    @test_throws ErrorException convert(TVectorE{RnTVector{Float64},RnPoint{Float64}},ξ)
    @test_throws ErrorException convert(TVectorE{SnTVector{Float64},SnPoint{Float64}},ξE)
    @test strip(ξE) == ξ
    @test getValue(ξE+ξE) == getValue(ξE)+getValue(ξE)
    @test_throws DomainError checkBasePoint(ξE, TVectorE(ξ,RnPoint([0., 1.]) ) )
    # promote
    @test promote(RnPoint(1.), MPointE(RnPoint(1.))) == (MPointE(RnPoint(1.)), MPointE(RnPoint(1.)))
    # getValue should be transparend
    @test getValue(xE) == getValue(x)
    @test getValue(ξE) == getValue(ξ)
    @test "$(xE)" == "$(x)E"
    @test "$(ξE)" == "$(ξ)E_$x"
    @test strip(x) == x
    @test strip(ξ) == ξ
    @test checkBasePoint(ξE,ξE)
    @test checkBasePoint(ξ,ξE)
    @test checkBasePoint(ξE,ξ)
    @test checkBasePoint(ξE,x)
    @test_throws DomainError checkBasePoint(ξE,RnPoint([0., 1.]))
    @test checkBasePoint(ξ,x)
    @test getValue(2. *ξE) == 2*getValue(ξ)
    @test getValue(ξE*2.) == getValue(ξ)*2
    @test getValue(ξE/2) == getValue(ξ)/2
    @test getValue(ξE/2) == getValue(ξ)/2
    @test getValue(-ξE) == -getValue(ξ)
    @test getValue(ξE-ξE) == getValue(TVectorE(ξ-ξ,x))
    @test getValue(+ξE) == getValue(ξ)
    @test ξE == ξE #check wether == for these types works
    # for the default functinos make sure they check and promote
    @test distance(M,xE,xE) == 0
    @test distance(M,x,xE) == 0
    @test distance(M,xE,x) == 0
    @test distance(M,x,MPointE(x,true)) == 0 # validate works
    @test_throws ErrorException distance( Sphere(2), SnPoint([1., 0., 0.]), MPointE( SnPoint([2., 0., 0.]),true) )
    # dot – all combinations
    @test dot(M,xE,ξE,ξE) == 1
    @test dot(M,xE,ξE,ξ) == 1
    @test dot(M,xE,ξ,ξE) == 1
    @test dot(M,xE,ξ,ξ) == 1
    @test dot(M,x,ξE,ξE) == 1
    @test dot(M,x,ξE,ξ) == 1
    @test dot(M,x,ξ,ξE) == 1
    @test dot(M,x,TVectorE(ξ,x,true),ξE) == 1
    # exp - all combinations
    y = RnPoint([1.,1.])
    yE = MPointE(y)
    @test exp(M,xE,ξE) == yE
    @test exp(M,MPointE(x,true),ξ) == yE
    @test exp(M,x,ξE) == yE
    # log
    @test log(M,x,yE) == ξE
    @test log(M,xE,y) == ξE
    @test log(M,xE,MPointE(y,true)) == ξE
    #
    @test manifoldDimension(xE) == 2
    #
    @test norm(M,x,ξE) == 1
    @test norm(M,xE,ξ) == 1
    @test norm(M,MPointE(x,true),ξE) == 1
    #
    ξEy = TVectorE(ξ,y)
    @test parallelTransport(M,xE,yE,ξE) == ξEy
    @test parallelTransport(M,xE,y,ξE) == ξEy
    @test parallelTransport(M,x,yE,ξE) == ξEy
    @test parallelTransport(M,MPointE(x,true),yE,ξ) == ξEy
    @test parallelTransport(M,xE,y,ξ) == ξEy
    @test parallelTransport(M,x,yE,ξ) == ξEy
    #
    @test validateTVector(M,MPointE(x,true),randomTVector(M,MPointE(x,true)))
    #
    #    @test tangentONB(M,xE,y)
    #
    validateMPoint(M,xE)
    validateTVector(M,xE,ξE)
    validateTVector(M,x,ξE)
    validateTVector(M,xE,ξ)
    #
    A,b = tangentONB(M,x,y)
    AE = TVectorE.(A,Ref(x))
    AE1,b1 = tangentONB(M,xE,yE)
    @test all( [ (AE1 .== AE)... ,b1==b] )
    AE2,b2 = tangentONB(M,x,yE)
    @test all( [ (AE2 .== AE)... ,b2==b] )
    AE3,b3 = tangentONB(M,xE,y)
    @test all( [ (AE3 .== AE)... ,b3==b] )
    AE4, b4 = tangentONB(M,xE,ξE)
    @test all( [ (AE4 .== AE)... ,b4==b] )
    AE5, b5 = tangentONB(M,xE,ξ)
    @test all( [ (AE5 .== AE)... ,b5==b] )
    AE6, b6 = tangentONB(M,x,ξE)
    @test all( [ (AE6 .== AE)... ,b6==b] )
    @test zeroTVector(M,MPointE(x,true)) == TVectorE(RnTVector([0.,0.]),x)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Vulkan_Headers_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Vulkan_Headers")
JLLWrappers.@generate_main_file("Vulkan_Headers", UUID("8d446b21-f3ad-5576-a034-752265b9b6f9"))
end  # module Vulkan_Headers_jll

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
using IterativeSolvers
using LinearAlgebra
using Random

export nlar
export default_eigval_sorter
export residual_eigval_sorter
export threshold_eigval_sorter

"""
    function nlar([eltype],nep::ProjectableNEP,[orthmethod=ModifiedGramSchmidt],[neigs=10],[errmeasure],[tol=eps(real(T))*100],[maxit=100],[λ0=0],[v0=randn(T,size(nep,1))],[logger=0],[linsolvercreator=DefaultLinSolverCreator()],[R=0.01],[eigval_sorter=residual_eigval_sorter],[qrfact_orth=false],[max_subspace=100],[num_restart_ritz_vecs=8],[inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

The function implements the Nonlinear Arnoldi method, which finds `neigs` eigenpairs (or throws a `NoConvergenceException`) by projecting the problem to a subspace that is expanded in the course  of the algorithm.
The basis is orthogonalized either by using the QR method if `qrfact_orth` is `true` or else by an orthogonalization method `orthmethod`).
This entails solving a smaller projected problem using a method specified by `inner_solver_method`.
The logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.
(`λ0`,`v0`) is the initial guess for the eigenpair. `linsolvercreator` specifies how the linear system is created and solved.
`R` is a parameter used by the function specified by `eigval_sorter` to reject those ritz values that are within a distance `R` from any of the converged eigenvalues, so that repeated convergence to the same eigenpair can be avoided.
`max_subspace` is the maximum allowable size of the basis befor the algorithm restarts using a basis made of `num_restart_ritz_vecs` ritz vectors and the eigenvectors that the algorithm has converged to.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("dep0_tridiag");
julia> λ,v=nlar(nep,tol=1e-7,neigs=1,maxit=100,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ[1],v))
8.00192341259751e-7
```

# References
* H. Voss, An Arnoldi method for nonlinear eigenvalue problems. BIT. Numer. Math. 44: 387-401, 2004.

"""
nlar(nep::NEP;params...) = nlar(ComplexF64,nep::NEP;params...)
function nlar(::Type{T},
            nep::ProjectableNEP;
            orthmethod::Type{T_orth} = ModifiedGramSchmidt,
            neigs::Int=10,                                     #Number of eigenvalues required
            errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
            tol = eps(real(T))*100,
            maxit::Int = 100,
            λ::Number = zero(T),
            v::Vector = randn(T,size(nep,1)),
            logger = 0,
            linsolvercreator=DefaultLinSolverCreator(),
            R = 0.01,
            eigval_sorter::Function = residual_eigval_sorter, #Function to sort eigenvalues of the projected NEP
            qrfact_orth::Bool = false,
            max_subspace::Int = 100,                           #Maximum subspace size before we implement restarting
            num_restart_ritz_vecs::Int=8,
            inner_solver_method = DefaultInnerSolver(),
            inner_logger = 0) where {T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod}

        @parse_logger_param!(logger)
        @parse_logger_param!(inner_logger)

        #Check if maxit is larger than problem size
        if (maxit > size(nep,1))
            @warn "Maximum iteration count maxit=$maxit larger than problem size n=$(size(nep,1)). Reducing maxit."
            maxit = size(nep,1);
        end

        #Check if number of ritz vectors for restating is greater than neigs
        if(num_restart_ritz_vecs > neigs)
            @warn "Number of ritz vectors for restarting num_restart_ritz_vecs=$num_restart_ritz_vecs larger than neigs=$neigs. Reducing num_restart_ritz_vecs."
            num_restart_ritz_vecs = neigs;
        end

        #Check if maximum allowable subspace size is lesser than num_restart_ritz_vecs
        if(max_subspace < num_restart_ritz_vecs)
            @warn "Maximum subspace max_subspace=$max_subspace smaller than number of restarting ritz vectors num_restart_ritz_vecs=$num_restart_ritz_vecs. Increasing max_subspace"
            max_subspace = num_restart_ritz_vecs+20; #20 is hardcoded.
        end

        local σ::T = T(λ); #Initial pole
        n = size(nep,1)
        λ::T = T(λ)
        nu::T = λ
        u::Vector{T} = v

        #Initialize the basis V_1
        V::Matrix{T}= zeros(T,n,max_subspace);
        X::Matrix{T} = zeros(T,n,neigs);
        V[:,1] = normalize(v);
        cbs = 1;#Current basis size

        D::Vector{T} = zeros(T,neigs);#To store the converged eigenvalues
        err_hist=eps()*ones(maxit,neigs) ;#Error history

        Z::Matrix{T} = zeros(T,n,neigs+num_restart_ritz_vecs);#The matrix used for constructing the restarted basis
        m = 0; #Number of converged eigenvalues

        k = 1;

        proj_nep = create_proj_NEP(nep,maxit,T);

        local linsolver::LinSolver=create_linsolver(linsolvercreator,nep,λ)

        err = Inf;

        push_info!(logger, "Using inner solver $inner_solver_method");

        while ((m < neigs) && (k < maxit))
            Vk = view(V,:,1:cbs)
            # Construct and solve the small projected PEP projected problem (V^H)T(λ)Vx = 0
            expand_projectmatrices!(proj_nep,Vk,Vk);

            #Use inner_solve() to solve the smaller projected problem
            push_info!(logger, 2, "Solving inner problem",continues=true)
            dd,vv = inner_solve(inner_solver_method,T,proj_nep,neigs=neigs,σ=σ,inner_logger=inner_logger);
            push_info!(logger, 2, ". Done.")
            # Sort the eigenvalues of the projected problem
            nuv,yv = eigval_sorter(nep,dd,vv,σ,D,R,Vk);
            nu = nuv[1]; y=yv[:,1];

            if (isinf(nu))
                error("We did not find any (non-converged) eigenvalues to target")
            end

            #Determine ritz vector and residual
            u[:] = Vk*y; # Note: y and u are vectors (not matrices)

            #Normalize and compute residual
            normalize!(u);
            res::AbstractVector = compute_Mlincomb(nep,nu,u);

            #Check for convergence of one of the eigenvalues
            err = estimate_error(errmeasure,nu,u);

            push_iteration_info!(logger,k,λ=nu,v=u,err=err);
            err_hist[k,m+1]=err;
            if(err < tol)
                mplusone=m+1;
                push_info!(logger,
                           "****** $mplusone converged to eigenvalue: $nu errmeasure:$err");

                #Add to the set of converged eigenvalues and eigenvectors
                D[m+1] = nu;
                X[:,m+1] = u;

                ## Sort and select he eigenvalues of the projected problem as described before
                nuv,yv = eigval_sorter(nep,dd,vv,σ,D,R,Vk);
                nu1=nuv[1];
                y1=yv[:,1];

                #Compute residual again
                u1 = Vk*y1;
                normalize!(u1);
                res = compute_Mlincomb(nep,nu1,u1);

                m = m+1;
            end

            #Check if basis size has exceeded max_subspace. If yes, then restart.
            if(size(Vk,2) >= max_subspace)
                cbs = m+num_restart_ritz_vecs;

                #Construct the new basis
                Zv = view(Z,:,1:cbs);
                Zv[:,1:m] = X[:,1:m]; #Set the first m vectors of the restarted basis to be the converged eigenvectors
                Zv[:,m+1:cbs] = Vk*yv[:,1:num_restart_ritz_vecs]; #Set the rest to be the ritz vectors of the projected problem

                Q,_ = qr(Zv); #Thin QR-factorization
                V[:,1:cbs] = Matrix(Q);

            else
                #Compute new vector Δv to add to the search space V(k+1) = (Vk,Δv)
                Δv=lin_solve(linsolver,res)
                #Orthogonalize and normalize
                if (qrfact_orth)
                    # Orthogonalize the entire basis matrix
                    # together with Δv using QR-method.
                    # Slow but robust.
                    Q,_ = qr(hcat(Vk,Δv))
                    Q = Matrix(Q)

                    cbs = cbs+1;
                    V[:,1:cbs]=Q;
                else
                    h=zeros(T,k);
                    orthogonalize_and_normalize!(Vk,Δv,h,orthmethod);

                    #Expand basis
                    cbs = cbs+1;
                    V[:,cbs] = Δv;
                end
            end
            k = k+1;
        end

        #Throw no convergence exception if enough eigenvalues were not found.
        if(k >= maxit && m < neigs)
            msg="Number of iterations exceeded. maxit=$(maxit) and only $(m) eigenvalues converged out of $(neigs)."
            throw(NoConvergenceException(nu,u,err,msg))
        end

        return D,X,err_hist;
    end

###############################################################################################################
# Ritz value discarder:
# This function is called as a first step by all the sorter functions to discard ritz values that lie within
# a certain distance R of the converged eigenvalues.
function discard_ritz_values!(dd,D,R)
    for i=1:size(dd,1)
        for j=1:size(D,1)
            if (abs(dd[i]-D[j])<R)
                dd[i]=Inf; #Discard all Ritz values within a particular radius R
            end
        end
    end

end

# Default ritzvalue sorter:
# First discard all Ritz values within a distance R of any of the converged eigenvalues(of the original problem).
# Then sort by distance from the shift and select the mm-th furthest value from the pole.

## D = already computed eigenvalues
## dd, vv eigenpairs of projected problem
## σ targets

function  default_eigval_sorter(nep::NEP,dd,vv,σ,D,R,Vk)
    dd2=copy(dd);

    #Discard ritz values within a distance R of the converged eigenvalues
    discard_ritz_values!(dd2,D,R)

    ii = sortperm(abs.(dd2-σ));

    nu = dd2[ii];
    y = vv[:,ii];

    return nu,y;
end

# Residual-based Ritz value sorter:
# First discard all Ritz values within a distance R of any of the converged eigenvalues(of the original problem).
# Then select that Ritz value which gives the mm-th minimum product of (residual and distance from pole).
function residual_eigval_sorter(nep::NEP,dd,vv,σ,D,R,Vk,errmeasure::ErrmeasureType = DefaultErrmeasure(nep))

    eig_res = zeros(size(dd,1));
    dd2=copy(dd);

    #Discard ritz values within a distance R of the converged eigenvalues
    discard_ritz_values!(dd2,D,R)

    #Compute residuals for each Ritz value
    for i=1:size(dd,1)
        eig_res[i] = estimate_error(errmeasure,dd[i],Vk*vv[:,i]);
    end

    #Sort according to methods
    ii = sortperm(eig_res .* abs.(dd2 .- σ))

    nu = dd[ii];
    y = vv[:,ii];

    return nu,y;
end

# Threshold residual based Ritz value sorter:
# Same as residual_eigval_sorter() except that errors above a certain threshold are set to the threshold.
function threshold_eigval_sorter(nep::NEP,dd,vv,σ,D,R,Vk,errmeasure::ErrmeasureType = DefaultErrmeasure(nep),threshold=0.1)

    eig_res = zeros(size(dd,1));
    dd2=copy(dd);

    #Discard ritz values within a distance R of the converged eigenvalues
    discard_ritz_values!(dd2,D,R)

    #Compute residuals for each Ritz value
    temp_res = 0;
    for i=1:size(dd,1)
        temp_res = error_measure(errmeasure,dd[i],Vk*vv[:,i]);
        if(temp_res > threshold)
            eig_res[i] = threshold;
        else
            eig_res[i] = temp_res;
        end
    end

    #Sort according to methods
    ii = sortperm(eig_res.*abs.(dd2-σ));

    nu = dd[ii];
    y = vv[:,ii];

    return nu,y;
end


text:
function LUSolve(A::Array,B::Vector)

      #A = the matrix of coefficients 'A' must be squared
  #B = The Vector 'B' right hand side vector Vertical
    m,n =size(A)
    v = size(B,1)
    if m != n
        error("Matrix must be squared")
    elseif m != v
        error("Vector Length must be equal to Matrix")
    end

L,U,P,p=DolittleLU(A)
    A=P*A
    B=P*B
    X=D=zeros(m)

    #Fordward Substitution
    D[1]=B[1]
    for i=2:m
       D[i]=(B[i] - (L[i,1:i-1]'*D[1:i-1]))

    end

         #Back Substitution
    X[m]=D[end]/U[end,end]

    for i=m-1:-1:1
       X[i]=(D[i] - (U[i,i+1:m]'*X[i+1:m]))/U[i,i]
    end

    return X
end


text:
include("testsetup.jl")

import InteractiveUtils

InteractiveUtils.versioninfo(stdout; verbose = true)

const LOOPVECTORIZATION_TEST = get(ENV, "LOOPVECTORIZATION_TEST", "all")

if LOOPVECTORIZATION_TEST == "all"
  NUMGROUPS = 5
  processes = Vector{Base.Process}(undef, NUMGROUPS)
  paths = Vector{String}(undef, NUMGROUPS)
  ios = Vector{IOStream}(undef, NUMGROUPS)
  tmp = tempdir();
  for i ∈ 1:NUMGROUPS
    path, io = mktemp(tmp)
    paths[i] = path
    ios[i] = io
    env = copy(ENV)
    env["LOOPVECTORIZATION_TEST"] = "part$i"
    env["JULIA_NUM_THREADS"] = string(Threads.nthreads())
    processes[i] = run(pipeline(setenv(`$(Base.julia_cmd()) $(@__FILE__) --project=$(Base.active_project())`, env), stderr = io, stdout = io), wait=false)
  end
  completed = fill(false, NUMGROUPS)
  while true
    for i ∈ eachindex(completed)
      completed[i] && continue
      proc = processes[i]
      process_running(proc) && continue
      completed[i] = true
      close(ios[i])
      println("Test Group: $i")
      run(`cat $(paths[i])`)
      println("\n\n")
    end
    all(completed) && break
    sleep(5)
  end
  @testset verbose=true "All" begin
    for (i,proc) ∈ enumerate(processes)
      @testset "part$i" begin
        @test success(proc)
      end
    end
  end
else
  include("grouptests.jl")
end


text:
function mh_results{T<:Integer}(ncur::T, ncar::T,
                                nruns::T, play_mh::Function)
    stickwins = 0
    switchwins = 0
    for i in 1:nruns
        (isstickwin, isswitchwin) = play_mh(ncur, ncar)
        if isstickwin
            stickwins += 1
        end
        if isswitchwin
            switchwins += 1
        end
    end
    return (stickwins/nruns, switchwins/nruns)
end

function mh_analytic{T<:Integer}(ncur::T, ncar::T)
    stickodds = ncar/ncur
    switchodds = (ncar - stickodds)/(ncur-2)
    return (stickodds, switchodds)
end

function show_odds{T<:Real}(a::T, b::T)
    @sprintf "   %.1f   %.1f     %.2f" 100.0*a 100*b 1.0*b/a
end

function show_simulation{T<:Integer}(ncur::T, ncar::T, nruns::T)
    println()
    print("Simulating a ", ncur, " door, ", ncar, " car ")
    println("Monty Hall problem with ", nruns, " runs.\n")

    println("   Solution   Stick  Switch  Improvement")

    (a, b) = mh_results(ncur, ncar, nruns, play_mh_literal)
    println(@sprintf("%10s: ", "Literal"), show_odds(a, b))

    (a, b) = mh_results(ncur, ncar, nruns, play_mh_clean)
    println(@sprintf("%10s: ", "Clean"), show_odds(a, b))

    (a, b) = mh_analytic(ncur, ncar)
    println(@sprintf("%10s: ", "Analytic"), show_odds(a, b))
    println()
    return nothing
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import LightGraphs
const lg = LightGraphs

const ba8 = lg.barabasi_albert(8, 5, 2)

include("nauty.jl")

g8 = nauty.lg_to_nauty(ba8)


text:
module ExtremeStats

include("anomalies.jl")
include("features.jl")

export get_anomalies, get_anomalies!, Extreme, load_X, label_Extremes, ExtremeList, Features, getFeatures, combineExtremes, sortby, writeExtremes, writeFeatures, writeTimeSeries
import Images.label_components
import NetCDF.ncread

type Extreme{T}
  index::Int64
  locs::Array{Int,2}
  zvalues::Array{T,1}
  tbounds::(Int,Int)
  lonbounds::(Int,Int)
  latbounds::(Int,Int)
end

type ExtremeList{T,U,V}
  extremes::Array{Extreme{T}}
  area::Vector{U}
  lons::Vector{V}
  lats::Vector{V}
end

include("io.jl")

getindex(el::ExtremeList,i...)=getindex(el.extremes,i...)
function sortby!(el::ExtremeList,x::Vector;rev=false)
    o=sortperm(x,rev=rev)
    el.extremes[:]=el.extremes[o]
end

random_x(nlon=200,nlat=200,N_years=7)=rand(Float32,nlon,nlat,N_years*46);

function getSeasStat(series,ilon,ilat,NY,NpY,msc,stdmsc,n)
    fill!(msc,0.0)
    fill!(stdmsc,0.0)
    fill!(n,0.0)
    for iyear=0:(NY-1)
        for iday = 1:NpY
            v=series[ilon,ilat,iyear*NpY+iday]
            if !isnan(v)
                stdmsc[iday]=stdmsc[iday]+v*v
                msc[iday]=msc[iday]+v
                n[iday]=n[iday]+1
            end
        end
    end
    for iday = 1:NpY
        if n[iday]>0
            stdmsc[iday]=sqrt(max(stdmsc[iday]/n[iday]-msc[iday]*msc[iday]/(n[iday]*n[iday]),zero(eltype(series))))
            msc[iday]=msc[iday]/n[iday]
        else
            msc[iday]=nan(eltype(series))
            stdmsc[iday]=nan(eltype(series))
        end
    end
    msc,stdmsc
end

function nanquantile(xtest,q;useabs=false)
    x = Array(eltype(xtest),length(xtest))
    nanquantile!(xtest,q,x,useabs=useabs)
end

function nanquantile!(xtest,q,x;useabs=false)
    useabs ? for i=1:length(xtest) x[i]=abs(xtest[i]) end : copy!(x,xtest)
    nNaN  = 0; for i=1:length(xtest) nNaN = isnan(xtest[i]) ? nNaN+1 : nNaN end
    lv=length(xtest)-nNaN
    lv==0 && return(nan(eltype(x)))
    index = 1 + (lv-1)*q
    lo = ifloor(index)
    hi = iceil(index)
    vals=select!(x,lo:hi)
    h=index - lo
    length(vals)==1 && return(vals[1])
    r = (1.0-h)*vals[1] + h*vals[2]
end

function label_Extremes{T}(x::Array{T,3};quantile::Number=0.0,mask::BitArray{3}=trues(1,1,1),circular::Bool=false,pattern::BitArray=trues(3,3,3),area=ones(Float32,size(x,2)),lons=linspace(0,360,size(x,1)),lats=linspace(90,-90,size(x,2)),normscope="global")

    size(x)==size(mask) || size(mask)==(1,1,1) || error("Sizes of x and mask must be the same")
    nlon,nlat,ntime = size(x)
    offs = circular ? 1 : 0

    if quantile==one(quantile) || quantile==zero(quantile)
        x2=trues(size(x,1)+offs,size(x,2),size(x,3))
    else
        cmpfun = (quantile < 0.5) ? .< : .>
        #Then allocate bitArray that hold the true/falses
        x2=BitArray(size(x,1)+offs,size(x,2),size(x,3))
        if normscope=="global"
          #First calculate threshold
          tres = nanquantile(x,quantile)
          #Fill BitArray
          x2[1:nlon,:,:]=cmpfun(x,tres)
        elseif normscope=="local"
          qar1      = zeros(eltype(x),ntime)
          qar2      = zeros(eltype(x),ntime)
          for ilat=1:nlat, ilon=1:nlon
            for itime = 1:ntime
              qar1[itime] = x[ilon,ilat,itime]
            end
            tres=nanquantile!(qar1,quantile,qar2)
            x2[ilon,ilat,:]=cmpfun(qar1,tres)

          end
        end
    end

    if size(mask) != (1,1,1)
        presum=sum(x2)
        for i=1:size(mask,3), j=1:size(mask,2), k=1:size(mask,1)
            x2[k,j,i]=x2[k,j,i] && mask[k,j,i]
        end
        aftersum=sum(x2)
        println("$(presum-aftersum) of $(presum) elements were removed by additional condition ()$((presum-aftersum)/presum)%)")
    end

    if circular
      #Now attach first slice to the end (circular globe)
      x2[nlon+1,:,:]=x2[1,:,:];
    end
  #Do the labelling
  lx=label_components(x2,pattern);
  x2=0;gc();
  circular && renameLabels(lx,x);
  el=ExtremeList(x,lx,area,lons,lats);
  return(el)
end

function indices2List(larr,xarr,extremeList)
    curind=ones(Int,length(extremeList))
    for i=1:size(xarr,3), j=1:size(xarr,2), k=1:size(xarr,1)
        if larr[k,j,i]>0
            extremeList[larr[k,j,i]].locs[curind[larr[k,j,i]],1]=k
            extremeList[larr[k,j,i]].locs[curind[larr[k,j,i]],2]=j
            extremeList[larr[k,j,i]].locs[curind[larr[k,j,i]],3]=i
            extremeList[larr[k,j,i]].zvalues[curind[larr[k,j,i]]]=xarr[k,j,i]
            curind[larr[k,j,i]]=curind[larr[k,j,i]]+1
    end
  end
    return(curind)
end

function getTbounds(el::ExtremeList)
  for e in el.extremes
    tmin=typemax(Int)
    tmax=0
    lonmin=typemax(Int)
    lonmax=0
    latmin=typemax(Int)
    latmax=0
    for i=1:length(e.zvalues)
      if e.locs[i,3]<tmin tmin=e.locs[i,3] end
      if e.locs[i,3]>tmax tmax=e.locs[i,3] end
      if e.locs[i,1]<lonmin lonmin=e.locs[i,1] end
      if e.locs[i,1]>lonmax lonmax=e.locs[i,1] end
      if e.locs[i,2]<latmin latmin=e.locs[i,2] end
      if e.locs[i,2]>latmax latmax=e.locs[i,2] end
    end
    e.tbounds=(tmin,tmax)
    e.lonbounds=(lonmin,lonmax)
    e.latbounds=(latmin,latmax)
  end
end

typealias FeatureVector{T} Vector{T}

function ExtremeList{T}(x::Array{T,3},lx::Array{Int,3},area=ones(Float32,size(x,2)),lons=linspace(0,360,size(x,1)),lats=linspace(90,-90,size(x,2)))
  nEx=maximum(lx)
  numCells=countNumCell(lx,nEx)
  nempty=sum(numCells.==0);
  extremeList=[Extreme(i,zeros(Int,numCells[i],3),Array(eltype(x),numCells[i]),(0,0),(0,0),(0,0)) for i=1:nEx];
  indices2List(lx,x,extremeList);
  o=sortperm(numCells,rev=true)
  extremeList=extremeList[o];
  deleteat!(extremeList,(nEx-nempty+1):nEx);
  return(ExtremeList(extremeList,area,lons,lats))
end

function combineExtremes(elin::ExtremeList;nEx=length(elin.extremes))
    totlen=0
    for i=1:nEx
        totlen=totlen+length(elin.extremes[i].zvalues)
    end
    locs=Array(Int,totlen,3)
    zvalues=Array(eltype(elin.extremes[1].zvalues),totlen)
    k=1
    for i=1:nEx
        for j=1:length(elin.extremes[i].zvalues)
            locs[k,1]=elin.extremes[i].locs[j,1]
            locs[k,2]=elin.extremes[i].locs[j,2]
            locs[k,3]=elin.extremes[i].locs[j,3]
            zvalues[k]=elin.extremes[i].zvalues[j]
            k=k+1
        end
    end
    mins=minimum(locs,1)
    maxs=maximum(locs,1)
    e=Extreme(1,locs,zvalues,(mins[3],maxs[3]),(mins[1],maxs[1]),(mins[2],maxs[2]))
    elout=ExtremeList([e],elin.area,elin.lons,elin.lats)
end

function countNumCell(labelList,nEx)
  lAr=zeros(Int,nEx)
  for i=1:length(labelList)
    j=labelList[i]
    if j>0
      lAr[j]=lAr[j]+1
    end
  end
  return lAr
end

function renameLabels(lx,x)
  #Check if we have to relabel, if longitudes are padded
  nlon=size(x,1)
  nlat=size(x,2)
  ntime=size(x,3)
  size(lx,1)==size(x,1) && return(lx)
  size(lx,1)==(nlon+1) || error("Something is wrong with the lon dimensions")
  size(lx,2)==nlat || error("Something is wrong with the lat dimensions")
  size(lx,3)==ntime || error("Something is wrong with the time dimensions")
  renames=Dict{Int,Int}()
  for k=1:nlat, t=1:ntime
    if (lx[nlon+1,k,t]>0) && (!haskey(renames,lx[nlon+1,k,t]))
      renames[lx[nlon+1,k,t]]=lx[1,k,t]
    end
  end
  i=nlon+1;sren=1;
  while sren>0 && i>0
    sren=0;
    for k=1:nlat, l=1:(ntime)
      if haskey(renames,lx[i,k,l])
        lx[i,k,l]=renames[lx[i,k,l]]
        sren=sren+1;
      end
    end
    i=i-1;
  end
  lx[nlon+1,:,:]=0
end

function getFeatures(el::ExtremeList,features...)
    el.extremes[1].tbounds==(0,0) && getTbounds(el)
    myf       = ExtremeStats.Features.calcFeatureFunction(features...)
    prearrays = ExtremeStats.Features.getPreArrays(length(el.extremes[1].zvalues),eltype(el.extremes[1].zvalues),features...)
    eval(myf)
    retar     = [Array(Features.rettype(f,el),length(el.extremes)) for f in features]
    for i=1:length(el.extremes)
        ret=getFeatures(el.extremes[i],el.area,el.lons,el.lats,prearrays)
        for j=1:length(ret) retar[j][i]=ret[j] end
    end
    return retar
end

end # module

max_stars_repo_name: CAB-LAB/ExtremeStats

max_issues_repo_name: CAB-LAB/ExtremeStats

max_forks_repo_name: CAB-LAB/ExtremeStats


text:
import MathOptInterface
const MOI = MathOptInterface

MOI.Utilities.@product_of_sets(
	Cones,
	MOI.Zeros,
	MOI.Nonnegatives,
	MOI.Nonpositives,
	MOI.SecondOrderCone)
const OptimizerCache{T} = MOI.Utilities.GenericModel{
	Float64,
	MOI.Utilities.ObjectiveContainer{Float64},
	MOI.Utilities.VariablesContainer{Float64},
	MOI.Utilities.MatrixOfConstraints{
		Float64,
		MOI.Utilities.MutableSparseMatrixCSC{Float64, T, MOI.Utilities.OneBasedIndexing},
		Vector{Float64},
		Cones{Float64}
	}
}

mutable struct Optimizer <: MOI.AbstractOptimizer
	name::String
	cones::Union{Nothing, Cones{Float64}}
	problem::Union{Nothing, Problem}
	objective_constant::Float64
	state::Union{Nothing, State}
	sets::Union{Nothing, PIPG.Cones{Float64}}
	silent::Bool
	elapsed_time::Float64
	max_sense::Bool

	scaling

	niters::Int
	ϵ::Float64
	γ::Float64
	function Optimizer(; niters=1000000, ϵ=1e-9, γ=0.9)
		return new("", nothing, nothing, 0.0, nothing, nothing, false, 0.0, false, nothing, niters, ϵ, γ)
	end
end

# basics
MOI.is_empty(opt::Optimizer) = opt.cones === nothing
function MathOptInterface.empty!(opt::Optimizer)
	opt.cones = nothing
	opt.problem = nothing
	opt.state = nothing
	opt.sets = nothing
end

# silent
MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.set(o::Optimizer, ::MOI.Silent, value::Bool) = (o.silent = value)
MOI.get(o::Optimizer, ::MOI.Silent) = o.silent

MOI.supports(::Optimizer, ::MOI.Name) = true
MOI.set(o::Optimizer, ::MOI.Name, value::String) = (o.name = value)
MOI.get(o::Optimizer, ::MOI.Name) = o.name

MOI.get(o::Optimizer, ::MOI.SolverName) = "PIPG"
MOI.get(o::Optimizer, ::MOI.SolverVersion) = "v0.1.0"
MOI.get(o::Optimizer, ::MOI.RawSolver) = o.problem
MOI.get(o::Optimizer, ::MOI.SolveTimeSec) = o.elapsed_time

MOI.supports(::Optimizer, ::MOI.ConstraintBasisStatus) = false
# supported sets
MOI.supports(::Optimizer, ::Union{
		MOI.ObjectiveSense,
		MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
		MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}}) = true
MOI.supports_constraint(
	::Optimizer,
	::Type{MOI.VectorAffineFunction{Float64}},
	::Type{<: Union{
		MOI.Zeros,
		MOI.Nonnegatives,
		MOI.Nonpositives,
		MOI.SecondOrderCone
	}}) = true

function MOI.modify(
    model::MOI.ModelLike,
    bridge::MOI.Bridges.Constraint.VectorizeBridge,
    change::MOI.ScalarConstantChange)
    MOI.modify(
        model,
        bridge.vector_constraint,
        MOI.VectorConstantChange([-change.new_constant]))
    return
end

convert_cone(z::MOI.Zeros) = Zeros{Float64, z.dimension}()
convert_cone(n::MOI.Nonnegatives) = POCone{Float64, n.dimension}()
convert_cone(n::MOI.Nonpositives) = NOCone{Float64, n.dimension}()
convert_cone(s::MOI.SecondOrderCone) = SOCone{Float64, s.dimension}()

function build_problem(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache{T}}) where {T}
	Ab = src.model.constraints
	A = Ab.coefficients
	dest.cones = Ab.sets

	for (F, S) in keys(src.constraints)
		throw(MOI.UnsupportedConstraint{F,S}())
    end

    model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
    quad_objective_function = MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}()
    aff_objective_function = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
    for attr in model_attributes
        if attr != MOI.ObjectiveSense() && attr != quad_objective_function && attr != aff_objective_function
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    max_sense = false
    if MOI.ObjectiveSense() in model_attributes
        max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    end

    objective_constant = 0.0
    q = zeros(A.n)
    if quad_objective_function in model_attributes
    	obj = MOI.get(src, quad_objective_function)
        objective_constant += MOI.constant(obj)

        I = Vector{Int}()
        J = Vector{Int}()
        V = Vector{Float64}()

        for term in obj.quadratic_terms
        	push!(I, Int(term.variable_1.value))
        	push!(J, Int(term.variable_2.value))
        	push!(V, (max_sense ? -1 : 1) * term.coefficient)
        	if (term.variable_1 != term.variable_2) # off-diagonal
	        	push!(I, Int(term.variable_2.value))
	        	push!(J, Int(term.variable_1.value))
	        	push!(V, (max_sense ? -1 : 1) * term.coefficient)
        	end
        end
        for term in obj.affine_terms
            q[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient
        end
        P = sparse(I, J, V, A.n, A.n)
    else
    	P = spzeros(A.n, A.n)
    end

    if aff_objective_function in model_attributes
    	obj = MOI.get(src, aff_objective_function)
        objective_constant += MOI.constant(obj)
    	for term in obj.terms
            q[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient
    	end
    end

    function _map_sets(f, ::Type{T}, sets, ::Type{S}) where {T,S}
	    F = MOI.VectorAffineFunction{Float64}
	    cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}())
	    return T[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis]
	end
	cones = Cone[]
	for set in MOI.Utilities.set_types(Ab.sets)
		append!(cones, _map_sets(convert_cone, Cone, Ab, set))
	end
	k = PTCone{Float64}((cones..., ))
    d = Reals{Float64, A.n}()
    H = convert(SparseMatrixCSC{Float64, Int64}, A)
    g = -1 .* Ab.constants

    p = Problem(k, d, H, P, MVector{A.n, Float64}(q), MVector{A.m, Float64}(g), objective_constant)
    if dest.scaling != nothing
    	s = State(p; scaling=dest.scaling(p))
	else
		s = State(p)
	end
	dest.max_sense = max_sense
    dest.problem = p
    dest.state = s
    dest.sets = Ab.sets
end

# modification routines
check_constructed(d::Optimizer) = if isnothing(d.problem) || isnothing(d.state) error("Must optimize before modification!") end
function MOI.modify(
	d::Optimizer,
	ci::MOI.ConstraintIndex{T} where T<:Union{MOI.ScalarAffineFunction, MOI.VectorAffineFunction},
	change::Union{MOI.ScalarConstantChange, MOI.VectorConstantChange})
	check_constructed(d)
	rows = MOI.Utilities.rows(d.sets, ci)
	d.problem.g[rows] .= -1.0 .* d.state.row_scale[rows] .* change.new_constant
end
function MOI.modify(
	d::Optimizer,
	ci::MOI.ConstraintIndex{T} where T<:MOI.VectorAffineFunction,
	change::MOI.MultirowChange)
	new_coeiffs = change.new_coefficients
	rows = MOI.Utilities.rows(d.sets, ci)[getindex.(new_coeiffs, 1)]
	scale = d.state.row_scale[rows] .* d.state.col_scale[change.variable.value]
	d.problem.H[change.variable.value, rows] .= scale .* getindex.(new_coeiffs, 2)
end

# optimizer
function MOI.optimize!(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache{T}}) where {T}
	index_map = MOI.Utilities.identity_index_map(src)

	if isnothing(dest.problem) || isnothing(dest.state)
		build_problem(dest, src)
	end

	Ab = src.model.constraints
	A = Ab.coefficients
	PIPG.scale(dest.problem, dest.state)
	α = compute_α(dest.problem, dest.γ)
    res = @timed pipg(dest.problem, dest.state, dest.niters, α, dest.ϵ, SVector{A.n, Float64}(zeros(A.n)), SVector{A.m, Float64}(zeros(A.m)))
    dest.elapsed_time = res[2]
end
function MOI.optimize!(dest::Optimizer)
	if isnothing(dest.problem) || isnothing(dest.state)
		error("Optimizer not initialized!")
	end

	n = length(dest.problem.q)
	m = length(dest.problem.g)
	PIPG.scale(dest.problem, dest.state)
	α = compute_α(dest.problem, dest.γ)
    res = @timed pipg(dest.problem, dest.state, dest.niters, α, dest.ϵ, SVector{n, Float64}(zeros(n)), SVector{m, Float64}(zeros(m)))
    println("niters=$(res[1])")
    dest.elapsed_time = res[2]
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = MOI.Utilities.UniversalFallback(OptimizerCache{Int64}())
    index_map = MOI.copy_to(cache, src)
    MOI.optimize!(dest, cache)
    return index_map, true
end

function MOI.get(o::Optimizer, ::MOI.TerminationStatus)
	if isnothing(o.state) return MOI.OPTIMIZE_NOT_CALLED end
	s = o.state.solver_state
	if s == INDEFINITE
		return MOI.OPTIMIZE_NOT_CALLED
	elseif s == OPTIMAL
		return MOI.OPTIMAL
	elseif s == PRIMAL_INFEASIBLE
		return MOI.INFEASIBLE
	elseif s == DUAL_INFEASIBLE
		return MOI.DUAL_INFEASIBLE
	elseif s == TIMEOUT
		return MOI.INTERRUPTED
	end
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return ""
end

function MOI.get(o::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(o, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif !isnothing(o.state) && o.state.solver_state == OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif !isnothing(o.state) && o.state.solver_state == PRIMAL_INFEASIBLE
        return MOI.INFEASIBLE_POINT
    elseif !isnothing(o.state) && o.state.solver_state == DUAL_INFEASIBLE
        return MOI.INFEASIBILITY_CERTIFICATE
    end
    return MOI.NO_SOLUTION
end
function MOI.get(o::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(o, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif !isnothing(o.state) && o.state.solver_state == OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif !isnothing(o.state) && o.state.solver_state == PRIMAL_INFEASIBLE
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif !isnothing(o.state) && o.state.solver_state == DUAL_INFEASIBLE
        return MOI.INFEASIBLE_POINT
    end
    return MOI.NO_SOLUTION
end

MOI.get(::Optimizer, ::MOI.ResultCount) = 1

function MOI.get(o::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(o, attr)
	raw_value = objective_value(o.problem, o.state)
    return (o.max_sense ? -1.0 : 1.0) * raw_value + o.problem.c
end

function MOI.get(o::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(o, attr)
    return (o.max_sense ? -1.0 : 1.0) * dual_objective_value(o.problem, o.state) + o.problem.c
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.state.primal[vi.value] * optimizer.state.col_scale[vi.value]
end

function MOI.get(
	optimizer::Optimizer,
	attr::MOI.ConstraintDual,
	ci::MOI.ConstraintIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    rows = MOI.Utilities.rows(optimizer.cones, ci)
    return .- optimizer.state.dual[rows] .* optimizer.state.row_scale[rows]
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashlatex

"""
    dashlatex(;kwargs...)
    dashlatex(children::Any;kwargs...)
    dashlatex(children_maker::Function;kwargs...)

A DashLatex component.
DashLatex is a Dash wrapper for ReactLatex (https://github.com/zzish/react-latex) which is itself a wrapper for Katex (https://katex.org/docs/options.html).
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The LaTeX markup to be rendered.
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `displayMode` (Bool; optional): displayMode: boolean (default: false). If true the math will be rendered in display mode. If false the math will be rendered in inline mode. Differences between the two modes include:
Display mode starts in \\displaystyle, so \\int and \\sum are large, for example; while inline mode starts in \\textstyle, where subscripts and superscripts usually don't stack on top of operators like \\sum. You can always manually switch between \\displaystyle and \\textstyle using those commands.
Display mode centers math on its on line and disables automatic line breaking (though you can customize this behavior with custom CSS). In inline mode, KaTeX allows line breaks after outermost relations (like = or <) or binary operators (like + or \\times), the same as TeX.
- `errorColor` (String; optional): errorColor: string. A color string given in the format "#XXX" or "#XXXXXX". This option determines the color that unsupported
commands and invalid LaTeX are rendered in when throwOnError is set to false. (default: #cc0000)
- `fleqn` (Bool; optional): fleqn: boolean. If true, display math renders flush left with a 2em left margin, like \\documentclass[fleqn] in LaTeX with the amsmath package.
- `globalGroup` (Bool; optional): globalGroup: boolean (default: false). Place KaTeX code in the global group. As a consequence,
\\def and \\newcommand persist in macros across render calls. In LaTeX, constructs such as \\begin{equation}
and \$\$ create a local group and prevent definitions from becoming visible outside of those blocks.
- `leqno` (Bool; optional): leqno: boolean. If true, display math has \\tags rendered on the left instead of the right, like \\usepackage[leqno]{amsmath} in LaTeX.
- `macros` (Dict; optional): macros: object. A collection of custom macros. Each macro is a property with a name like \\name (written "\\\\name" in JavaScript)
which maps to a string that describes the expansion of the macro, or a function that accepts an instance of MacroExpander as first
argument and returns the expansion as a string. MacroExpander is an internal API and subject to non-backwards compatible changes.
See src/macros.js for its usage. Single-character keys can also be included in which case the character will be redefined as the
given macro (similar to TeX active characters). This object will be modified if the LaTeX code defines its own macros via \\gdef, which
enables consecutive calls to KaTeX to share state.
- `maxExpand` (Real; optional): maxExpand: number. Limit the number of macro expansions to the specified number, to prevent e.g. infinite macro loops.
If set to Infinity, the macro expander will try to fully expand as in LaTeX. (default: 1000)
- `maxSize` (Real; optional): maxSize: number. All user-specified sizes, e.g. in \\rule{500em}{500em}, will be capped to maxSize ems. If set to Infinity
(the default), users can make elements and spaces arbitrarily large.
- `minRuleThickness` (Real; optional): minRuleThickness: number. Specifies a minimum thickness, in ems, for fraction lines, \\sqrt top lines, {array} vertical lines,
\\hline, \\hdashline, \\underline, \\overline, and the borders of \\fbox, \\boxed, and \\fcolorbox. The usual value for these items
is 0.04, so for minRuleThickness to be effective it should probably take a value slightly above 0.04, say 0.05 or 0.06.
Negative values will be ignored.
- `strict` (Bool | String; optional): strict: boolean or string or function (default: "warn"). If false or "ignore", allow features that make writing LaTeX
convenient but are not actually supported by (Xe)LaTeX (similar to MathJax). If true or "error" (LaTeX faithfulness mode),
throw an error for any such transgressions. If "warn" (the default), warn about such behavior via console.warn.
Provide a custom function handler(errorCode, errorMsg, token) to customize behavior depending on the type of transgression
(summarized by the string code errorCode and detailed in errorMsg); this function can also return "ignore", "error", or "warn"
to use a built-in behavior. A list of such features and their errorCodes:
- `throwOnError` (Bool; optional): throwOnError: boolean. If true (the default), KaTeX will throw a ParseError when it encounters an unsupported command or invalid LaTeX.
If false, KaTeX will render unsupported commands as text, and render invalid LaTeX as its source code with hover text giving the error,
in the color given by errorColor.
"""
function dashlatex(; kwargs...)
        available_props = Symbol[:children, :id, :displayMode, :errorColor, :fleqn, :globalGroup, :leqno, :macros, :maxExpand, :maxSize, :minRuleThickness, :strict, :throwOnError]
        wild_props = Symbol[]
        return Component("dashlatex", "DashLatex", "dash_latex", available_props, wild_props; kwargs...)
end

dashlatex(children::Any; kwargs...) = dashlatex(;kwargs..., children = children)
dashlatex(children_maker::Function; kwargs...) = dashlatex(children_maker(); kwargs...)

max_stars_repo_name: ndrezn/dash-latex

max_issues_repo_name: ndrezn/dash-latex

max_forks_repo_name: ndrezn/dash-latex