text:
# Attributes

"""
    AbstractOptimizerAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of the optimizer.

### Note

The difference between `AbstractOptimizerAttribute` and `AbstractModelAttribute`
lies in the behavior of `is_empty`, `empty!` and `copy_to`. Typically optimizer
attributes only affect how the model is solved.
"""
abstract type AbstractOptimizerAttribute end

"""
    AbstractModelAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of the model.
"""
abstract type AbstractModelAttribute end

"""
    AbstractVariableAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of variables in the model.
"""
abstract type AbstractVariableAttribute end

"""
    AbstractConstraintAttribute

Abstract supertype for attribute objects that can be used to set or get attributes (properties) of constraints in the model.
"""
abstract type AbstractConstraintAttribute end

const AnyAttribute = Union{AbstractOptimizerAttribute, AbstractModelAttribute, AbstractVariableAttribute, AbstractConstraintAttribute}

# This allows to use attributes in broadcast calls without the need to
# embed it in a `Ref`
Base.broadcastable(attribute::AnyAttribute) = Ref(attribute)

"""
    struct UnsupportedAttribute{AttrType} <: UnsupportedError
        attr::AttrType
        message::String
    end

An error indicating that the attribute `attr` is not supported by the model,
i.e. that [`supports`](@ref) returns `false`.
"""
struct UnsupportedAttribute{AttrType<:AnyAttribute} <: UnsupportedError
    attr::AttrType
    message::String
end
UnsupportedAttribute(attr::AnyAttribute) = UnsupportedAttribute(attr, "")
element_name(err::UnsupportedAttribute) = "Attribute $(err.attr)"

"""
    struct SetAttributeNotAllowed{AttrType} <: NotAllowedError
        attr::AttrType
        message::String # Human-friendly explanation why the attribute cannot be set
    end

An error indicating that the attribute `attr` is supported (see
[`supports`](@ref)) but cannot be set for some reason (see the error string).
"""
struct SetAttributeNotAllowed{AttrType<:AnyAttribute} <: NotAllowedError
    attr::AttrType
	message::String # Human-friendly explanation why the attribute cannot be set
end
SetAttributeNotAllowed(attr::AnyAttribute) = SetAttributeNotAllowed(attr, "")

operation_name(err::SetAttributeNotAllowed) = "Setting attribute $(err.attr)"
message(err::SetAttributeNotAllowed) = err.message

"""
    supports(model::ModelLike, attr::AbstractOptimizerAttribute)::Bool

Return a `Bool` indicating whether `model` supports the optimizer attribute
`attr`. That is, it returns `false` if `copy_to(model, src)` shows a warning in
case `attr` is in the [`ListOfOptimizerAttributesSet`](@ref) of `src`; see
[`copy_to`](@ref) for more details on how unsupported optimizer attributes are
handled in copy.

    supports(model::ModelLike, attr::AbstractModelAttribute)::Bool

Return a `Bool` indicating whether `model` supports the model attribute `attr`.
That is, it returns `false` if `copy_to(model, src)` cannot be performed in case
`attr` is in the [`ListOfModelAttributesSet`](@ref) of `src`.

    supports(model::ModelLike, attr::AbstractVariableAttribute, ::Type{VariableIndex})::Bool

Return a `Bool` indicating whether `model` supports the variable attribute
`attr`. That is, it returns `false` if `copy_to(model, src)` cannot be performed
in case `attr` is in the [`ListOfVariableAttributesSet`](@ref) of `src`.

    supports(model::ModelLike, attr::AbstractConstraintAttribute, ::Type{ConstraintIndex{F,S}})::Bool where {F,S}

Return a `Bool` indicating whether `model` supports the constraint attribute
`attr` applied to an `F`-in-`S` constraint. That is, it returns `false` if
`copy_to(model, src)` cannot be performed in case `attr` is in the
[`ListOfConstraintAttributesSet`](@ref) of `src`.

For all four methods, if the attribute is only not supported in specific
circumstances, it should still return `true`.

Note that `supports` is only defined for attributes for which
[`is_copyable`](@ref) returns `true` as other attributes do not appear in the
list of attributes set obtained by `ListOf...AttributesSet`.
"""
function supports end
function supports(::ModelLike, attr::Union{AbstractModelAttribute,
                                           AbstractOptimizerAttribute})
    if !is_copyable(attr)
        throw(ArgumentError("`supports` is not defined for $attr, it is only" *
                            " defined for attributes such that `is_copyable`" *
                            " returns `true`."))
    end
    return false
end
function supports(::ModelLike, attr::Union{AbstractVariableAttribute,
                                           AbstractConstraintAttribute},
                  ::Type{<:Index})
    if !is_copyable(attr)
        throw(ArgumentError("`supports` is not defined for $attr, it is only" *
                            " defined for attributes such that `is_copyable`" *
                            " returns `true`."))
    end
    return false
end

"""
    get(optimizer::AbstractOptimizer, attr::AbstractOptimizerAttribute)

Return an attribute `attr` of the optimizer `optimizer`.

    get(model::ModelLike, attr::AbstractModelAttribute)

Return an attribute `attr` of the model `model`.

    get(model::ModelLike, attr::AbstractVariableAttribute, v::VariableIndex)

If the attribute `attr` is set for the variable `v` in the model `model`, return
its value, return `nothing` otherwise. If the attribute `attr` is not supported
by `model` then an error should be thrown instead of returning `nothing`.

    get(model::ModelLike, attr::AbstractVariableAttribute, v::Vector{VariableIndex})

Return a vector of attributes corresponding to each variable in the collection `v` in the model `model`.

    get(model::ModelLike, attr::AbstractConstraintAttribute, c::ConstraintIndex)

If the attribute `attr` is set for the constraint `c` in the model `model`,
return its value, return `nothing` otherwise. If the attribute `attr` is not
supported by `model` then an error should be thrown instead of returning
`nothing`.

    get(model::ModelLike, attr::AbstractConstraintAttribute, c::Vector{ConstraintIndex{F,S}})

Return a vector of attributes corresponding to each constraint in the collection `c` in the model `model`.

    get(model::ModelLike, ::Type{VariableIndex}, name::String)

If a variable with name `name` exists in the model `model`, return the
corresponding index, otherwise return `nothing`. Errors if two variables
have the same name.

    get(model::ModelLike, ::Type{ConstraintIndex{F,S}}, name::String) where {F<:AbstractFunction,S<:AbstractSet}

If an `F`-in-`S` constraint with name `name` exists in the model `model`, return
the corresponding index, otherwise return `nothing`. Errors if two constraints
have the same name.

    get(model::ModelLike, ::Type{ConstraintIndex}, name::String)

If *any* constraint with name `name` exists in the model `model`, return the
corresponding index, otherwise return `nothing`. This version is available for
convenience but may incur a performance penalty because it is not type stable.
Errors if two constraints have the same name.

### Examples

```julia
get(model, ObjectiveValue())
get(model, VariablePrimal(), ref)
get(model, VariablePrimal(5), [ref1, ref2])
get(model, OtherAttribute("something specific to cplex"))
get(model, VariableIndex, "var1")
get(model, ConstraintIndex{ScalarAffineFunction{Float64},LessThan{Float64}}, "con1")
get(model, ConstraintIndex, "con1")
```
"""
function get end
# We want to avoid being too specific in the type arguments to avoid method ambiguity.
# For model, get(::ModelLike, ::AbstractVariableAttribute, ::Vector{VariableIndex}) would not allow
# to define get(::SomeModel, ::AnyAttribute, ::Vector)
get(model::ModelLike, attr::AnyAttribute, idxs::Vector) = get.(model, attr, idxs)

function get(model::ModelLike, attr::AnyAttribute, args...)
    throw(ArgumentError("ModelLike of type $(typeof(model)) does not support accessing the attribute $attr"))
end

"""
    get!(output, model::ModelLike, args...)

An in-place version of `get`.
The signature matches that of `get` except that the the result is placed in the vector `output`.
"""
function get! end
function get!(output, model::ModelLike, attr::AnyAttribute, args...)
    throw(ArgumentError("ModelLike of type $(typeof(model)) does not support accessing the attribute $attr"))
end

"""
    set(optimizer::AbstractOptimizer, attr::AbstractOptimizerAttribute, value)

Assign `value` to the attribute `attr` of the optimizer `optimizer`.

    set(model::ModelLike, attr::AbstractModelAttribute, value)

Assign `value` to the attribute `attr` of the model `model`.

    set(model::ModelLike, attr::AbstractVariableAttribute, v::VariableIndex, value)

Assign `value` to the attribute `attr` of variable `v` in model `model`.

    set(model::ModelLike, attr::AbstractVariableAttribute, v::Vector{VariableIndex}, vector_of_values)

Assign a value respectively to the attribute `attr` of each variable in the collection `v` in model `model`.

    set(model::ModelLike, attr::AbstractConstraintAttribute, c::ConstraintIndex, value)

Assign a value to the attribute `attr` of constraint `c` in model `model`.

    set(model::ModelLike, attr::AbstractConstraintAttribute, c::Vector{ConstraintIndex{F,S}}, vector_of_values)

Assign a value respectively to the attribute `attr` of each constraint in the collection `c` in model `model`.

An [`UnsupportedAttribute`](@ref) error is thrown if `model` does not support
the attribute `attr` (see [`supports`](@ref)) and a [`SetAttributeNotAllowed`](@ref)
error is thrown if it supports the attribute `attr` but it cannot be set.

### Replace set in a constraint

    set(model::ModelLike, ::ConstraintSet, c::ConstraintIndex{F,S}, set::S)

Change the set of constraint `c` to the new set `set` which should be of the
same type as the original set.

#### Examples

If `c` is a `ConstraintIndex{F,Interval}`

```julia
set(model, ConstraintSet(), c, Interval(0, 5))
set(model, ConstraintSet(), c, GreaterThan(0.0))  # Error
```

### Replace function in a constraint

    set(model::ModelLike, ::ConstraintFunction, c::ConstraintIndex{F,S}, func::F)

Replace the function in constraint `c` with `func`. `F` must match the original
function type used to define the constraint.

#### Note

Setting the constraint function is not allowed if `F` is
[`SingleVariable`](@ref), it throws a
[`SettingSingleVariableFunctionNotAllowed`](@ref) error. Indeed, it would
require changing the index `c` as the index of `SingleVariable` constraints
should be the same as the index of the variable.

#### Examples

If `c` is a `ConstraintIndex{ScalarAffineFunction,S}` and `v1` and `v2` are
`VariableIndex` objects,

```julia
set(model, ConstraintFunction(), c,
    ScalarAffineFunction(ScalarAffineTerm.([1.0, 2.0], [v1, v2]), 5.0))
set(model, ConstraintFunction(), c, SingleVariable(v1)) # Error
```
"""
function set end
# See note with get
function set(model::ModelLike,
             attr::Union{AbstractVariableAttribute,
                         AbstractConstraintAttribute},
             idxs::Vector, vector_of_values::Vector)
    if length(idxs) != length(vector_of_values)
        throw(DimensionMismatch("Number of indices ($(length(idxs))) does " *
                                "not match the number of values " *
                                "($(length(vector_of_values))) set to `$attr`."))
    end
    return set.(model, attr, idxs, vector_of_values)
end

function set(model::ModelLike, attr::AnyAttribute, args...)
    throw_set_error_fallback(model, attr, args...)
end
# throw_set_error_fallback is included so that we can return type-specific error
# messages without needing to overload set and cause ambiguity errors. For
# examples, see ConstraintSet and ConstraintFunction. throw_set_error_fallback should
# not be overloaded by users of MOI.
function throw_set_error_fallback(model::ModelLike,
                                  attr::Union{AbstractModelAttribute,
                                              AbstractOptimizerAttribute},
                                  value;
                                  error_if_supported = SetAttributeNotAllowed(attr))
    if supports(model, attr)
        throw(error_if_supported)
    else
        throw(UnsupportedAttribute(attr))
    end
end
function throw_set_error_fallback(model::ModelLike,
                                  attr::Union{AbstractVariableAttribute,
                                              AbstractConstraintAttribute},
                                  index::Index, value;
                                  error_if_supported = SetAttributeNotAllowed(attr))
    if supports(model, attr, typeof(index))
        throw(error_if_supported)
    else
        throw(UnsupportedAttribute(attr))
    end
end

"""
    ListOfOptimizerAttributesSet()

Error type that should be thrown when the user [`set`](@ref) the
[`ConstraintFunction`](@ref) of a [`SingleVariable`](@ref) constraint.
"""
struct SettingSingleVariableFunctionNotAllowed <: Exception end

## Optimizer attributes

"""
    ListOfOptimizerAttributesSet()

An optimizer attribute for the `Vector{AbstractOptimizerAttribute}` of all optimizer attributes that were set.
"""
struct ListOfOptimizerAttributesSet <: AbstractOptimizerAttribute end

"""
    SolverName()

An optimizer attribute for the string identifying the solver/optimizer.
"""
struct SolverName <: AbstractOptimizerAttribute end

"""
    Silent

An optimizer attribute for silencing the output of an optimizer. When `set`
to `true`, it takes precedence over any other attribute controlling verbosity
and requires the solver to produce no output. The default value is `false`
which has no effect. In this case the verbosity is controlled by other
attributes.

## Note

Every optimizer should have verbosity on by default. For instance, if a solver
has a solver-specific log level attribute, the MOI implementation should set it
to `1` by default. If the user sets `Silent` to `true`, then the log level
should be set to `0`, even if the user specifically sets a value of log level.
If the value of `Silent` is `false` then the log level set to the solver is the
value given by the user for this solver-specific parameter or `1` if none is
given.
"""
struct Silent <: AbstractOptimizerAttribute end

"""
    RawParameter(name)

An optimizer attribute for the solver-specific parameter identified by `name`
which is typically an `Enum` or a `String`.
"""
struct RawParameter <: AbstractOptimizerAttribute
    name::Any
end

## Model attributes

"""
    ListOfModelAttributesSet()

A model attribute for the `Vector{AbstractModelAttribute}` of all model
attributes `attr` such that 1) `is_copyable(attr)` returns `true` and 2) the
attribute was set to the model.
"""
struct ListOfModelAttributesSet <: AbstractModelAttribute end

"""
    Name()

A model attribute for the string identifying the model. It has a default value
of `""` if not set`.
"""
struct Name <: AbstractModelAttribute end

"""
    ObjectiveSense()

A model attribute for the `OptimizationSense` of the objective function, which
can be `MIN_SENSE`, `MAX_SENSE`, or `FeasiblitySense`. The default is
`FEASIBILITY_SENSE`.
"""
struct ObjectiveSense <: AbstractModelAttribute end

@enum OptimizationSense MIN_SENSE MAX_SENSE FEASIBILITY_SENSE

"""
    NumberOfVariables()

A model attribute for the number of variables in the model.
"""
struct NumberOfVariables <: AbstractModelAttribute end

"""
    ListOfVariableIndices()

A model attribute for the `Vector{VariableIndex}` of all variable indices present in the model
(i.e., of length equal to the value of `NumberOfVariables()`) in the order in
which they were added.
"""
struct ListOfVariableIndices <: AbstractModelAttribute end

"""
    ListOfConstraintIndices{F,S}()

A model attribute for the `Vector{ConstraintIndex{F,S}}` of all constraint indices of type
`F`-in-`S` in the model (i.e., of length equal to the value of
`NumberOfConstraints{F,S}()`) in the order in which they were added.
"""
struct ListOfConstraintIndices{F,S} <: AbstractModelAttribute end

"""
    NumberOfConstraints{F,S}()

A model attribute for the number of constraints of the type `F`-in-`S` present in the model.
"""
struct NumberOfConstraints{F,S} <: AbstractModelAttribute end

"""
    ListOfConstraints()

A model attribute for the list of tuples of the form `(F,S)`, where `F` is a function type
and `S` is a set type indicating that the attribute `NumberOfConstraints{F,S}()`
has value greater than zero.
"""
struct ListOfConstraints <: AbstractModelAttribute end

"""
    ObjectiveFunction{F<:AbstractScalarFunction}()

A model attribute for the objective function which has a type `F<:AbstractScalarFunction`.
`F` should be guaranteed to be equivalent but not necessarily identical to the function type provided by the user.
Throws an `InexactError` if the objective function cannot be converted to `F`,
e.g. the objective function is quadratic and `F` is `ScalarAffineFunction{Float64}` or
it has non-integer coefficient and `F` is `ScalarAffineFunction{Int}`.
"""
struct ObjectiveFunction{F<:AbstractScalarFunction} <: AbstractModelAttribute end

"""
    ObjectiveFunctionType()

A model attribute for the type `F` of the objective function set using the
`ObjectiveFunction{F}` attribute.

## Examples

In the following code, `attr` should be equal to `MOI.SingleVariable`:
```julia
x = MOI.add_variable(model)
MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
         MOI.SingleVariable(x))
attr = MOI.get(model, MOI.ObjectiveFunctionType())
```
"""
struct ObjectiveFunctionType <: AbstractModelAttribute end

## Optimizer attributes

"""
    ObjectiveValue(resultidx::Int=1)

A model attribute for the objective value of the `result_index`th primal result.
"""
struct ObjectiveValue <: AbstractModelAttribute
    result_index::Int
    (::Type{ObjectiveValue})(result_index=1) = new(result_index)
end

"""
    DualObjectiveValue(result_index::Int=1)

A model attribute for the value of the objective function of the dual problem
for the `result_index`th dual result.
"""
struct DualObjectiveValue <: AbstractModelAttribute
    result_index::Int
    (::Type{DualObjectiveValue})(result_index=1) = new(result_index)
end

"""
    ObjectiveBound()

A model attribute for the best known bound on the optimal objective value.
"""
struct ObjectiveBound <: AbstractModelAttribute end

"""
    RelativeGap()

A model attribute for the final relative optimality gap, defined as ``\\frac{|b-f|}{|f|}``, where ``b`` is the best bound and ``f`` is the best feasible objective value.
"""
struct RelativeGap <: AbstractModelAttribute  end

"""
    SolveTime()

A model attribute for the total elapsed solution time (in seconds) as reported by the optimizer.
"""
struct SolveTime <: AbstractModelAttribute end

"""
    SimplexIterations()

A model attribute for the cumulative number of simplex iterations during the optimization process.
In particular, for a mixed-integer program (MIP), the total simplex iterations for all nodes.
"""
struct SimplexIterations <: AbstractModelAttribute end

"""
    BarrierIterations()

A model attribute for the cumulative number of barrier iterations while solving a problem.
"""
struct BarrierIterations <: AbstractModelAttribute end

"""
    NodeCount()

A model attribute for the total number of branch-and-bound nodes explored while solving a mixed-integer program (MIP).
"""
struct NodeCount <: AbstractModelAttribute end

"""
    RawSolver()

A model attribute for the object that may be used to access a solver-specific API for this optimizer.
"""
struct RawSolver <: AbstractModelAttribute end

"""
    ResultCount()

A model attribute for the number of results available.
"""
struct ResultCount <: AbstractModelAttribute end

## Variable attributes

"""
    ListOfVariableAttributesSet()

A model attribute for the `Vector{AbstractVariableAttribute}` of all variable
attributes `attr` such that 1) `is_copyable(attr)` returns `true` and 2) the
attribute was set to variables.
"""
struct ListOfVariableAttributesSet <: AbstractModelAttribute end

"""
    VariableName()

A variable attribute for a string identifying the variable. It is *valid* for
two variables to have the same name; however, variables with duplicate names
cannot be looked up using [`get`](@ref). It has a default value of `""` if not
set`.
"""
struct VariableName <: AbstractVariableAttribute end

"""
    VariablePrimalStart()

A variable attribute for the initial assignment to some primal variable's value that the optimizer may use to warm-start the solve.
"""
struct VariablePrimalStart <: AbstractVariableAttribute end

"""
    VariablePrimal(N)
    VariablePrimal()

A variable attribute for the assignment to some primal variable's value in result `N`.
If `N` is omitted, it is 1 by default.
"""
struct VariablePrimal <: AbstractVariableAttribute
    N::Int
end
VariablePrimal() = VariablePrimal(1)

"""
    BasisStatusCode

An Enum of possible values for the `ConstraintBasisStatus` attribute.
This explains the status of a given element with respect to an optimal solution basis.
Possible values are:
* `BASIC`: element is in the basis
* `NONBASIC`: element is not in the basis
* `NONBASIC_AT_LOWER`: element is not in the basis and is at its lower bound
* `NONBASIC_AT_UPPER`: element is not in the basis and is at its upper bound
* `SUPER_BASIC`: element is not in the basis but is also not at one of its
  bounds

Note: `NONBASIC_AT_LOWER` and `NONBASIC_AT_UPPER` should be used only for
constraints with the `Interval`. In this case cases they are necessary to
distinguish which side of the constraint. One-sided constraints (e.g.,
`LessThan` and `GreaterThan`) should use `NONBASIC` instead of the
`NONBASIC_AT_*` values.
"""
@enum(BasisStatusCode, BASIC, NONBASIC, NONBASIC_AT_LOWER, NONBASIC_AT_UPPER,
      SUPER_BASIC)

## Constraint attributes

"""
    ListOfConstraintAttributesSet{F, S}()

A model attribute for the `Vector{AbstractConstraintAttribute}` of all
constraint attributes `attr` such that 1) `is_copyable(attr)` returns `true` and
2) the attribute was set to `F`-in-`S` constraints.

## Note

The attributes [`ConstraintFunction`](@ref) and [`ConstraintSet`](@ref) should
not be included in the list even if then have been set with [`set`](@ref).
"""
struct ListOfConstraintAttributesSet{F,S} <: AbstractModelAttribute end

"""
    ConstraintName()

A constraint attribute for a string identifying the constraint. It is *valid*
for constraints variables to have the same name; however, constraints with
duplicate names cannot be looked up using [`get`](@ref) regardless of if they
have the same `F`-in-`S` type. It has a default value of `""` if not set.
"""
struct ConstraintName <: AbstractConstraintAttribute end

"""
    ConstraintPrimalStart()

A constraint attribute for the initial assignment to some constraint's primal value(s) that the optimizer may use to warm-start the solve.
"""
struct ConstraintPrimalStart <: AbstractConstraintAttribute end

"""
    ConstraintDualStart()

A constraint attribute for the initial assignment to some constraint's dual value(s) that the optimizer may use to warm-start the solve.
"""
struct ConstraintDualStart <: AbstractConstraintAttribute end

"""
    ConstraintPrimal(N)
    ConstraintPrimal()

A constraint attribute for the assignment to some constraint's primal value(s) in result `N`.
If `N` is omitted, it is 1 by default.

Given a constraint `function-in-set`, the `ConstraintPrimal` is the value of the
function evaluated at the primal solution of the variables. For example, given
the constraint `ScalarAffineFunction([x,y], [1, 2], 3)`-in-`Interval(0, 20)` and
a primal solution of `(x,y) = (4,5)`, the `ConstraintPrimal` solution of the
constraint is `1 * 4 + 2 * 5 + 3 = 17`.
"""
struct ConstraintPrimal <: AbstractConstraintAttribute
    N::Int
end
ConstraintPrimal() = ConstraintPrimal(1)

"""
    ConstraintDual(N)
    ConstraintDual()

A constraint attribute for the assignment to some constraint's dual value(s) in result `N`.
If `N` is omitted, it is 1 by default.
"""
struct ConstraintDual <: AbstractConstraintAttribute
    N::Int
end
ConstraintDual() = ConstraintDual(1)

"""
    ConstraintBasisStatus()

A constraint attribute for the `BasisStatusCode` of some constraint, with
respect to an available optimal solution basis.

**For the basis status of a variable, query the corresponding `SingleVariable`
constraint that enforces the variable's bounds.**
"""
struct ConstraintBasisStatus <: AbstractConstraintAttribute end

"""
    ConstraintFunction()

A constraint attribute for the `AbstractFunction` object used to define the constraint.
It is guaranteed to be equivalent but not necessarily identical to the function provided by the user.
"""
struct ConstraintFunction <: AbstractConstraintAttribute end

function throw_set_error_fallback(::ModelLike, attr::ConstraintFunction,
                                  ::ConstraintIndex{F, S}, ::F;
                                  error_if_supported = SetAttributeNotAllowed(attr)) where {F <: AbstractFunction, S}
    throw(error_if_supported)
end
func_type(c::ConstraintIndex{F, S}) where {F, S} = F
function throw_set_error_fallback(::ModelLike, ::ConstraintFunction,
                                  constraint_index::ConstraintIndex,
                                  func::AbstractFunction;
                                  kwargs...)
    throw(ArgumentError("""Cannot modify functions of different types.
    Constraint type is $(func_type(constraint_index)) while the replacement
    function is of type $(typeof(func))."""))
end

"""
    ConstraintSet()

A constraint attribute for the `AbstractSet` object used to define the constraint.
"""
struct ConstraintSet <: AbstractConstraintAttribute end

function throw_set_error_fallback(::ModelLike, attr::ConstraintSet,
                                  ::ConstraintIndex{F, S}, ::S;
                                  error_if_supported = SetAttributeNotAllowed(attr)) where {F, S <: AbstractSet}
    throw(error_if_supported)
end
set_type(::ConstraintIndex{F, S}) where {F, S} = S
function throw_set_error_fallback(::ModelLike, ::ConstraintSet,
                                  constraint_index::ConstraintIndex,
                                  set::AbstractSet;
                                  kwargs...)
    throw(ArgumentError("""Cannot modify sets of different types. Constraint
    type is $(set_type(constraint_index)) while the replacement set is of
    type $(typeof(set)). Use `transform` instead."""))
end

## Termination status
"""
    TerminationStatus()

A model attribute for the `TerminationStatusCode` explaining why the optimizer stopped.
"""
struct TerminationStatus <: AbstractModelAttribute end

"""
    TerminationStatusCode

An Enum of possible values for the `TerminationStatus` attribute. This attribute
is meant to explain the reason why the optimizer stopped executing in the most
recent call to [`optimize!`](@ref).

If no call has been made to [`optimize!`](@ref), then the `TerminationStatus`
is:

* `OPTIMIZE_NOT_CALLED`: The algorithm has not started.

## OK

These are generally OK statuses, i.e., the algorithm ran to completion normally.

* `OPTIMAL`: The algorithm found a globally optimal solution.
* `INFEASIBLE`: The algorithm concluded that no feasible solution exists.
* `DUAL_INFEASIBLE`: The algorithm concluded that no dual bound exists for the
  problem. If, additionally, a feasible (primal) solution is known to
  exist, this status typically implies that the problem is unbounded, with some
  technical exceptions.
* `LOCALLY_SOLVED`: The algorithm converged to a stationary point, local
  optimal solution, could not find directions for improvement, or otherwise
  completed its search without global guarantees.
* `LOCALLY_INFEASIBLE`: The algorithm converged to an infeasible point or
  otherwise completed its search without finding a feasible solution, without
  guarantees that no feasible solution exists.
* `INFEASIBLE_OR_UNBOUNDED`: The algorithm stopped because it decided that the
  problem is infeasible or unbounded; this occasionally happens during MIP
  presolve.

## Solved to relaxed tolerances

* `ALMOST_OPTIMAL`: The algorithm found a globally optimal solution to relaxed
  tolerances.
* `ALMOST_INFEASIBLE`: The algorithm concluded that no feasible solution exists
  within relaxed tolerances.
* `ALMOST_DUAL_INFEASIBLE`: The algorithm concluded that no dual bound exists for
  the problem within relaxed tolerances.
* `ALMOST_LOCALLY_SOLVED`: The algorithm converged to a stationary point, local
  optimal solution, or could not find directions for improvement within relaxed
  tolerances.

## Limits

The optimizer stopped because of some user-defined limit.

* `ITERATION_LIMIT`: An iterative algorithm stopped after conducting the maximum
  number of iterations.
* `TIME_LIMIT`: The algorithm stopped after a user-specified computation time.
* `NODE_LIMIT`: A branch-and-bound algorithm stopped because it explored a
  maximum number of nodes in the branch-and-bound tree.
* `SOLUTION_LIMIT`: The algorithm stopped because it found the required number of
  solutions. This is often used in MIPs to get the solver to return the first
  feasible solution it encounters.
* `MEMORY_LIMIT`: The algorithm stopped because it ran out of memory.
* `OBJECTIVE_LIMIT`: The algorthm stopped because it found a solution better than
  a minimum limit set by the user.
* `NORM_LIMIT`: The algorithm stopped because the norm of an iterate became too
  large.
* `OTHER_LIMIT`: The algorithm stopped due to a limit not covered by one of the
  above.

## Problematic

This group of statuses means that something unexpected or problematic happened.

* `SLOW_PROGRESS`: The algorithm stopped because it was unable to continue making
  progress towards the solution.
* `NUMERICAL_ERROR`: The algorithm stopped because it encountered unrecoverable
  numerical error.
* `INVALID_MODEL`: The algorithm stopped because the model is invalid.
* `INVALID_OPTION`: The algorithm stopped because it was provided an invalid
  option.
* `INTERRUPTED`: The algorithm stopped because of an interrupt signal.
* `OTHER_ERROR`: The algorithm stopped because of an error not covered by one of
  the statuses defined above.
"""
@enum(TerminationStatusCode,
    OPTIMIZE_NOT_CALLED,
    # OK
    OPTIMAL, INFEASIBLE, DUAL_INFEASIBLE, LOCALLY_SOLVED, LOCALLY_INFEASIBLE,
    INFEASIBLE_OR_UNBOUNDED,
    # Solved to relaxed tolerances
    ALMOST_OPTIMAL, ALMOST_INFEASIBLE, ALMOST_DUAL_INFEASIBLE,
    ALMOST_LOCALLY_SOLVED,
    # Limits
    ITERATION_LIMIT, TIME_LIMIT,  NODE_LIMIT, SOLUTION_LIMIT, MEMORY_LIMIT,
    OBJECTIVE_LIMIT, NORM_LIMIT, OTHER_LIMIT,
    # Problematic
    SLOW_PROGRESS, NUMERICAL_ERROR, INVALID_MODEL, INVALID_OPTION, INTERRUPTED,
    OTHER_ERROR
)

"""
    RawStatusString()

A model attribute for a solver specific string explaining why the optimizer
stopped.
"""
struct RawStatusString <: AbstractModelAttribute end

## Result status

"""
    ResultStatusCode

An Enum of possible values for the `PrimalStatus` and `DualStatus` attributes.
The values indicate how to interpret the result vector.

* `NO_SOLUTION`: the result vector is empty.
* `FEASIBLE_POINT`: the result vector is a feasible point.
* `NEARLY_FEASIBLE_POINT`: the result vector is feasible if some constraint
  tolerances are relaxed.
* `INFEASIBLE_POINT`: the result vector is an infeasible point.
* `INFEASIBILITY_CERTIFICATE`: the result vector is an infeasibility certificate.
  If the `PrimalStatus` is `INFEASIBILITY_CERTIFICATE`, then the primal result
  vector is a certificate of dual infeasibility. If the `DualStatus` is
  `INFEASIBILITY_CERTIFICATE`, then the dual result vector is a proof of primal
  infeasibility.
* `NEARLY_INFEASIBILITY_CERTIFICATE`: the result satisfies a relaxed criterion for
  a certificate of infeasibility.
* `REDUCTION_CERTIFICATE`: the result vector is an ill-posed certificate; see
  [this article](https://arxiv.org/abs/1408.4685) for details. If the
  `PrimalStatus` is `REDUCTION_CERTIFICATE`, then the primal result vector is a
  proof that the dual problem is ill-posed. If the `DualStatus` is
  `REDUCTION_CERTIFICATE`, then the dual result vector is a proof that the primal
  is ill-posed.
* `NEARLY_REDUCTION_CERTIFICATE`: the result satisfies a relaxed criterion for
  an ill-posed certificate.
* `UNKNOWN_RESULT_STATUS`: the result vector contains a solution with an unknown
  interpretation.
* `OTHER_RESULT_STATUS`: the result vector contains a solution with an
  interpretation not covered by one of the statuses defined above.
"""
@enum(ResultStatusCode, NO_SOLUTION,
      FEASIBLE_POINT, NEARLY_FEASIBLE_POINT, INFEASIBLE_POINT,
      INFEASIBILITY_CERTIFICATE, NEARLY_INFEASIBILITY_CERTIFICATE,
      REDUCTION_CERTIFICATE, NEARLY_REDUCTION_CERTIFICATE,
      UNKNOWN_RESULT_STATUS, OTHER_RESULT_STATUS)

"""
    PrimalStatus(N)
    PrimalStatus()

A model attribute for the `ResultStatusCode` of the primal result `N`.
If `N` is omitted, it defaults to 1.
"""
struct PrimalStatus <: AbstractModelAttribute
    N::Int
end
PrimalStatus() = PrimalStatus(1)

"""
    DualStatus(N)
    DualStatus()

A model attribute for the `ResultStatusCode` of the dual result `N`.
If `N` is omitted, it defaults to 1.
"""
struct DualStatus <: AbstractModelAttribute
    N::Int
end
DualStatus() = DualStatus(1)

"""
    is_set_by_optimize(::AnyAttribute)

Return a `Bool` indicating whether the value of the attribute is modified
during an [`optimize!`](@ref) call, that is, the attribute is used to query
the result of the optimization.

## Important note when defining new attributes

This function returns `false` by default so it should be implemented for
attributes that are modified by [`optimize!`](@ref).
"""
is_set_by_optimize(::AnyAttribute) = false
function is_set_by_optimize(::Union{ObjectiveValue,
                                    ObjectiveBound,
                                    RelativeGap,
                                    SolveTime,
                                    SimplexIterations,
                                    BarrierIterations,
                                    NodeCount,
                                    RawSolver,
                                    ResultCount,
                                    TerminationStatus,
                                    RawStatusString,
                                    PrimalStatus,
                                    DualStatus,
                                    VariablePrimal,
                                    ConstraintPrimal,
                                    ConstraintDual,
                                    ConstraintBasisStatus})
    return true
end

"""
    is_copyable(::AnyAttribute)

Return a `Bool` indicating whether the value of the attribute may be copied
during [`copy_to`](@ref) using [`set`](@ref).

## Important note when defining new attributes

By default `is_copyable(attr)` returns `!is_set_by_optimize(attr)`. A specific
method should be defined for attributes which are copied indirectly during
[`copy_to`](@ref). For instance, both `is_copyable` and
[`is_set_by_optimize`](@ref) return `false` for the following attributes:

* [`ListOfOptimizerAttributesSet`](@ref), [`ListOfModelAttributesSet`](@ref),
  [`ListOfConstraintAttributesSet`](@ref) and
  [`ListOfVariableAttributesSet`](@ref).
* [`SolverName`](@ref) and [`RawSolver`](@ref): these attributes cannot be set.
* [`NumberOfVariables`](@ref) and [`ListOfVariableIndices`](@ref): these
  attributes are set indirectly by [`add_variable`](@ref) and
  [`add_variables`](@ref).
* [`ObjectiveFunctionType`](@ref): this attribute is set indirectly when setting
  the [`ObjectiveFunction`](@ref) attribute.
* [`NumberOfConstraints`](@ref), [`ListOfConstraintIndices`](@ref),
  [`ListOfConstraints`](@ref), [`ConstraintFunction`](@ref) and
  [`ConstraintSet`](@ref): these attributes are set indirectly by
  [`add_constraint`](@ref) and [`add_constraints`](@ref).
"""
function is_copyable(attr::AnyAttribute)
    return !is_set_by_optimize(attr)
end
function is_copyable(::Union{ListOfOptimizerAttributesSet,
                             ListOfModelAttributesSet,
                             ListOfConstraintAttributesSet,
                             ListOfVariableAttributesSet,
                             SolverName,
                             RawSolver,
                             NumberOfVariables,
                             ListOfVariableIndices,
                             NumberOfConstraints,
                             ObjectiveFunctionType,
                             ListOfConstraintIndices,
                             ListOfConstraints,
                             ConstraintFunction,
                             ConstraintSet})
    return false
end


text:
using Distributions
using Test
using DataFrames
using CSVFiles
using DataStructures
# using Plots

# using VegaLite

import GlobalSensitivityAnalysis: ishigami
import StatsBase: ordinalrank

################################################################################
## SETUP
################################################################################

# define the (uncertain) parameters of the problem and their distributions
data = DeltaData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = 1000
)

N = data.N
D = length(data.params)

@testset "Ishigami Function - Sampling" begin

    py_samples = load("data/delta/py_ishigami/py_samples.csv", header_exists=false, colnames = ["x1", "x2", "x3"]) |> DataFrame
    julia_samples = GlobalSensitivityAnalysis.sample(data)

    quants = [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]
    sig_level = 0.05
    output_dir = joinpath(@__DIR__, "../output/LHS_QuantileTesting/ishigami")
    mkpath(output_dir)

    for i = 1:D
        sampleA = py_samples[:,i]
        sampleB = julia_samples[:,i]

        # Make dataframe
        # df = DataFrame(A = sampleA, B = sampleB)
        # save(joinpath(output_dir, "samples$(N)_p$i.csv"), df)

        # # Quick plot comparison
        # p1 = df |> @vlplot(:bar, x = {:A, bin = {step=0.1}}, y="count()")
        # save(joinpath(output_dir, "sampleA_$(N)_p$i.png"), p1)
        # p2 = df |> @vlplot(:bar, x = {:B, bin = {step=0.1}}, y="count()", background=:white)
        # save(joinpath(output_dir, "sampleB_$(N)_p$i.png"), p2)

        # Run WRS quantile matching
        results = WRS.pb2gen(sampleA, sampleB, quantiles = quants) |> DataFrame
        # save(joinpath(output_dir, "LHS Sampling Quantile Comparison N$(N)_p$i.csv"), results)

        @test sum(results[!,:signif]) == 0

    end

end

@testset "Ishigami Function - Analysis" begin

    # Get the samples - here we will use the Python samples to keep things consistent,
    # could also use Julia samples run through both the SALib and GSA functions
    samples = load("data/delta/py_ishigami/py_samples.csv", header_exists=false, colnames = ["x1", "x2", "x3"]) |> DataFrame
    samples = Matrix(samples)

    # check ishigami
    python_Y = load("data/delta/py_ishigami/py_ishigami.csv", header_exists=false) |> DataFrame
    julia_Y = ishigami(samples)
    @test python_Y[:,1] ≈ julia_Y atol = ATOL_sample

    # julia results
    julia_results = analyze(data, samples, julia_Y; num_resamples = 1_000)

    # python results
    py_firstorder = load("data/delta/py_ishigami/py_firstorder.csv", header_exists=false) |> DataFrame
    py_delta = load("data/delta/py_ishigami/py_delta.csv", header_exists=false) |> DataFrame
    py_firstorder_conf = load("data/delta/py_ishigami/py_firstorder_conf.csv", header_exists=false) |> DataFrame
    py_delta_conf = load("data/delta/py_ishigami/py_delta_conf.csv", header_exists=false) |> DataFrame

    # test indices - check values and orderings
    @test julia_results[:firstorder] ≈ Matrix(py_firstorder) atol = ATOL_delta
    @test ordinalrank(julia_results[:firstorder]) == ordinalrank(py_firstorder[!,:Column1])

    @test julia_results[:delta] ≈ Matrix(py_delta) atol = 0.05 # TODO - this seems too high?
    @test ordinalrank(julia_results[:delta]) == ordinalrank(py_delta[!, :Column1])

    # test confidence intervals
    @test julia_results[:firstorder_conf] ≈ Matrix(py_firstorder_conf) atol = ATOL_CI
    @test julia_results[:delta_conf] ≈ Matrix(py_delta_conf) atol = ATOL_CI

end

ext: jl

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
using KFEstimate
using LinearAlgebra, Plots, Zygote, Statistics, Revise
using Flux, Flux.Optimise
using ProgressBars
pathof(KFEstimate)

dt = 0.001
A = [1.0 dt 1/2*dt^2; 0.0 1.0 dt; 0.0 0.0 1.0]
B = [0.0; 0.0; 1.0]
B = reshape(B, length(B), 1)
Q = 1*Matrix{Float64}(I, 3, 3)

# observation model, assume we can noisily measure position
H = randn(3, 3)
# H = [1 0 0; 0 1 0; 0 0 1]
R = 10.0*Matrix{Float64}(I, 3, 3)

kf = KalmanFilter(A, B, Q, H, R)

# run simulation
time_step = 0:1000
x0 = [0.0; 0.0; 0.0]
action_sequence = [[1] for t in time_step]
sim_states, sim_measurements = run_simulation(kf, x0, action_sequence)
# run kalman filter
P0 = Matrix{Float64}(I, 3, 3)
s0 = State(x0, P0) # initial state belief
filtered_states = run_filter(kf, s0, action_sequence, sim_measurements)

# unpack sim and filtered states
μ, Σ = unpack(filtered_states)
p = [x[1] for x in sim_states]
v = [x[2] for x in sim_states]
a = [x[3] for x in sim_states]
x = hcat(p, v, a)
# plot
plot(time_step, [x[2:end, 1] x[2:end, 2] x[2:end, 3]], label = ["simulated p" "simulated v" "simulated a"])
plot!(time_step, [μ[2:end, 1] μ[2:end, 2] μ[2:end, 3]], label = ["filtred p" "filtered v" "filtered a"])
xlabel!("time step (t)")

##

# parametrized matrix estimates
Ahat(θ) = [θ[1] dt 1/2*dt^2; 0.0 θ[1] dt; 0.0 0.0 θ[1]]
Bhat(θ) = B
Qhat(θ) = Q
Hhat(θ) = H
Rhat(θ) = R
# define a parametrized kalman filter
param_kf = ParamKalmanFilter(Ahat, Bhat, Qhat, Hhat, Rhat)

θ0 = [1.1]
pkf = ParamKalmanFilter(Ahat, Bhat, Qhat, Hhat, Rhat)

loss = []
θ_range = 0.9:0.001:1.1
for i in θ_range
    θ_i = [i]
    states = run_param_kf(θ_i, pkf, s0, action_sequence, sim_measurements)
    l = kf_likelihood(θ_i, param_kf, states, action_sequence, sim_measurements)
    push!(loss, l)
end

plot(θ_range, loss)

##

# parametrized matrix estimates
Ahat(θ) = [θ[1] θ[2] θ[3]; 0 θ[1] θ[2]; 0 0 θ[1]]
# Ahat(θ) = [θ[1, 1] θ[1, 2] θ[1, 3]; 0.0 θ[2, 2] θ[2, 3]; 0.0 0.0 θ[3, 3]]
Bhat(θ) = B
Qhat(θ) = Q
Hhat(θ) = H
Rhat(θ) = R
# define a parametrized kalman filter
param_kf = ParamKalmanFilter(Ahat, Bhat, Qhat, Hhat, Rhat)

θ0 = [1.0 dt 1/2*dt^2] + 0.1*randn(1, 3)
opt = ADAM(0.01)
epochs = 200
newθ, loss = run_kf_gradient(θ0, param_kf, s0, action_sequence, sim_measurements, opt, epochs)

grad_states = run_param_filter(newθ, param_kf, s0, action_sequence, sim_measurements)
μgrad, Σgrad = unpack(grad_states)

l = @layout [a{0.7h};grid(1, 3)]
p1 = plot(time_step, [x[2:end, 1] x[2:end, 2] x[2:end, 3]], label = ["simulated p" "simulated v" "simulated a"], xlabel="time step (t)")
p1 = plot!(time_step, [μ[2:end, 1] μ[2:end, 2] μ[2:end, 3]], label = ["filtered p" "filtered v" "filtered a"], xlabel="time step (t)")
p1 = plot!(time_step, [μgrad[2:end, 1] μgrad[2:end, 2] μgrad[2:end, 3]], label = ["learned p" "learned v" "learned a"],  xlabel="time step (t)")
p2 = plot(1:epochs, loss, title="loss", xlabel="number of epochs")
p3 = plot(time_step, (x[2:end, :]-μ[2:end, :]).^2, title="KF vs. sim error", xlabel="time step (t)")
p4 = plot(time_step, (x[2:end, :]-μgrad[2:end, :]).^2, title="grad vs. sim error", xlabel="time step (t)")
plot(p1, p2, p3, p4, layout=l, titlefont = font(12), size=(1000, 700))


text:
using Toxcore

tox_options_default()


text:
using Test
using DataDrivenControl
using LinearAlgebra

@testset "utils" begin
    @testset "convert" begin
        # convert_to_matrix
        n = 2
        w = rand(Int(n*(n+1)/2))
        P = DataDrivenControl.convert_to_matrix(w)
        @test P == [w[1] 0.5*w[2]; 0.5*w[2] w[3]]
        # convert_quadratic_to_linear_basis
        x = rand(n)
        w' * DataDrivenControl.convert_quadratic_to_linear_basis(x) == x' * P * x
    end
    @testset "cost" begin
        x = ones(3)
        u = ones(3)
        Q = Matrix(I, 3, 3)
        R = 2*Matrix(I, 3, 3)

        quadratic_cost = DataDrivenControl.QuadraticCost(Q, R)
        r = quadratic_cost(x, u)
        @test r == x'*Q*x + u'*R*u  # analytically obtained
        Q_func(x) = x'*Q*x
        R_func(x) = R
        quadratic_in_input_cost = DataDrivenControl.QuadraticInInputCost(Q_func, R_func)
        @test r == quadratic_in_input_cost(x, u)
    end
end


text:
@testset "Indexing" begin
    function run_tests(datype)
        data = rand(10, 10)
        na = bitrand(10, 10)
        A = datype(data, na)

        # Scalar getindex
        for i = 1:100
            if na[i]
                @test ismissing(A[i])
            else
                @test A[i] == data[i]
            end
        end
        for i = 1:10, j = 1:10
            if na[i, j]
                @test ismissing(A[i, j])
            else
                @test A[i, j] == data[i, j]
            end
        end

        # getindex with AbstractVectors
        rg = 2:9
        v = A[rg]
        for i = 1:length(rg)
            if na[rg[i]]
                @test ismissing(v[i])
            else
                @test v[i] == data[rg[i]]
            end
        end

        v = A[rg, 9]
        for i = 1:length(rg)
            if na[rg[i], 9]
                @test ismissing(v[i])
            else
                @test v[i] == data[rg[i], 9]
            end
        end

        rg2 = 5:7
        v = A[rg, rg2]
        for j = 1:length(rg2), i = 1:length(rg)
            if na[rg[i], rg2[j]]
                @test ismissing(v[i, j])
            else
                @test v[i, j] == data[rg[i], rg2[j]]
            end
        end

        # getindex with AbstractVector{Bool}
        b = bitrand(10, 10)
        rg = findall(b)
        v = A[b]
        for i = 1:length(rg)
            if na[rg[i]]
                @test ismissing(v[i])
            else
                @test v[i] == data[rg[i]]
            end
        end

        # getindex with DataVectors with missingness throws
        @test_throws MissingException A[@data([1, 2, 3, missing])]

        # setindex! with scalar indices
        data = rand(10, 10)
        for i = 1:100
            A[i] = data[i]
        end
        @test A == data

        data = rand(10, 10)
        for i = 1:10, j = 1:10
            A[i, j] = data[i, j]
        end
        @test A == data

        na = bitrand(10, 10)
        for i = 1:100
            na[i] && (A[i] = missing)
        end

        # setindex! with scalar and vector indices
        rg = 2:9
        data[rg] = 1.0
        A[rg] = 1.0
        for i = 1:length(rg)
            @test A[rg[i]] == 1.0
        end

        # setindex! with missing and vector indices
        rg = 5:13
        na[rg] = true
        A[rg] = missing
        for i = 1:length(rg)
            @test ismissing(A[rg[i]])
        end

        # setindex! with vector and vector indices
        rg = 12:67
        data[rg] = rand(length(rg))
        A[rg] = data[rg]
        for i = 1:length(rg)
            @test A[rg[i]] == data[rg[i]]
        end

        # setindex! with DataArray/PooledDataArray and vector indices
        for datype2 in (DataArray, PooledDataArray)
            newdata = rand(3, 3)
            newdata[1:2:9] = data[1:2:9]
            newna = bitrand(3, 3)
            rg1, rg2 = 1:3, 5:7
            data[rg1, rg2] = newdata
            na[rg1, rg2] = newna
            A[rg1, rg2] = datype2(newdata, newna)
            for j = rg2, i = rg1
                if na[i, j]
                    @test ismissing(A[i, j])
                else
                    @test A[i, j] == data[i, j]
                end
            end
        end
    end

    run_tests(DataArray)
    run_tests(PooledDataArray)
end


text:
# %% Entropy comparison
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

he = Systems.henon()
N, q = 10000, 2
X = trajectory(he, N, Ttr = 10)
ε = 10 .^ (-5.0:0.1:1)
H = [genentropy(X, e; q) for e in ε]
C = correlationsum(X, ε)

fig, axs = subplots(2, 1; figsize = (figx/2, figy), sharex = true)

axs[1].plot(log.(ε), -H, COLORS[1])

x = -5:0
D = 1.22

axs[2].plot(log.(ε), log.(C), COLORS[1])

for (j, (ax, x)) in enumerate(zip(axs, (-5:0, -10:0)))
    Dl = D .* x
    Dl = Dl .- 2
    ax.axvspan(x[1], x[end], color = COLORS[4], alpha = 0.25)
    ax.grid("on")
    ax.plot(x, Dl, c = COLORS[3])
    ax.text(x[2] + 1.5 + 0.5*(j-1), Dl[2], "\$D_2\$ = $D", color = COLORS[3], size = 30)
end

xtick = -11:3:1
xtickl = string.(xtick)
xtickl[3] = ""
axs[2].set_xticks(xtick)
axs[2].set_xticklabels(xtickl)
axs[2].set_xlabel("\$\\log(\\varepsilon)\$"; labelpad = -20)
# axs[1].set_yticks()
axs[1].set_ylabel("\$-H_$q\$")
axs[2].set_ylabel("\$\\log(C)\$")
axs[1].set_yticks(0:-3:-10)
axs[2].set_yticks(0:-4:-12)
fig.tight_layout(pad=0.3)
wsave(plotsdir("5", "dimension"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using MarketData
export simulate

# run model without closing price updates

# convert TimeArray to Timestamps array
clt = Timestamp[]

for i in 1:length(cl)
    push!(clt, Timestamp(cl[i].timestamp[1], cl[i].values[1]))
end

cr = vcat(res, clt) # combine both arrays with txn array first
sd = sort(Date[c.timestamp for c in cr]) # get dates and sort
ud = unique(sd) # remove duplicates

mixed_array = cr[ud[1]] # initialize array to be pushed

for i in 2:length(ud)
    mixed_array = vcat(mixed_array, cr[ud[i]])
end

# this takes mixed_array and updates the following
# txn & mark_price explicity
# A, Pval, Aval by association

function simulate(tsv::Vector{Timestamp})
    for t in 1:length(tsv)
        if typeof(tsv[t].value) <: FinancialAsset
            push!(txn, tsv[t].value)
        elseif typeof(tsv[t].value) == Float64
            push!(mark_price, tsv[t].value)
        else
            error("rogue timestamp in your midst!")
        end
    end
end


text:
abstract type BitMask{T<:Unsigned} end

Base.broadcastable(x::BitMask) = Ref(x)

function generate_bitmask_flags(type, decl)
    identifier, value = decl.args
    :(const $identifier = $type($value))
end

"""
    @bitmask_flag BitFlags::UInt32 begin
        FLAG_A = 1
        FLAG_B = 2
        FLAG_C = 4
    end

Enumeration of bitmask flags that can be combined with `&`, `|` and `xor`, forbidding the combination of flags from different bitmasks.
"""
macro bitmask_flag(typedecl, expr)
    if Base.is_expr(typedecl, :(::), 2)
        type, eltype = typedecl.args
    else
        error("First argument to @bitmask_flag must be of the form 'type::eltype'")
    end
    decls = filter(x -> typeof(x) ≠ LineNumberNode, expr.args)

    exs = [
        :(
            Base.@__doc__ struct $type <: BitMask{$eltype}
                val::$eltype
            end
        );
        generate_bitmask_flags.(type, decls)
    ]

    Expr(:block, esc.(exs)...)
end

(&)(a::BitMask, b::BitMask) = error("Bitwise operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")
(|)(a::BitMask, b::BitMask) = error("Bitwise operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")
xor(a::BitMask, b::BitMask) = error("Bitwise operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")
isless(a::BitMask, b::BitMask) = error("Bitwise operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")
(==)(a::BitMask, b::BitMask) = error("Operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")
in(a::BitMask, b::BitMask) = error("Operation not allowed between incompatible bitmasks '$(typeof(a))', '$(typeof(b))'")

(&)(a::T, b::T) where {T <: BitMask} = T(a.val & b.val)
(|)(a::T, b::T) where {T <: BitMask} = T(a.val | b.val)
xor(a::T, b::T) where {T <: BitMask} = T(xor(a.val, b.val))
isless(a::T, b::T) where {T <: BitMask} = isless(a.val, b.val)
(==)(a::T, b::T) where {T <: BitMask} = a.val == b.val
in(a::T, b::T) where {T <: BitMask} = a & b == a

(&)(a::T, b::Integer) where {T <: BitMask} = T(a.val & b)
(|)(a::T, b::Integer) where {T <: BitMask} = T(a.val | b)
xor(a::T, b::Integer) where {T <: BitMask} = T(xor(a.val, b))
isless(a::T, b::Integer) where {T <: BitMask} = isless(a.val, b)
in(a::T, b::Integer) where {T <: BitMask} = a & b == a

(&)(a::Integer, b::T) where {T <: BitMask} = b & a
(|)(a::Integer, b::T) where {T <: BitMask} = b | a
xor(a::Integer, b::T) where {T <: BitMask} = xor(b, a)
isless(a::Integer, b::T) where {T <: BitMask} = isless(a, b.val) # need b.val to prevent stackoverflow
in(a::Integer, b::T) where {T <: BitMask} = a | b == b

(::Type{T})(bm::BitMask) where {T <: Integer} = T(bm.val)

convert(T::Type{<:Integer}, bm::BitMask) = T(bm.val)
convert(T::Type{<:BitMask}, val::Integer) = T(val)

typemax(T::Type{<:BitMask{_T}}) where {_T} = T(typemax(_T))


text:
"""
    Bridges{T}

A state type for the depth-first search that finds bridges in a graph.
"""
mutable struct Bridges{T<:Integer}
    low::Vector{T}
    depth::Vector{T}
    bridges::Vector{Edge{T}}
    id::T
end

function Bridges(g::AbstractGraph)
    n = nv(g)
    T = eltype(g)
    return Bridges(zeros(T, n), zeros(T, n), Edge{T}[], zero(T))
end

"""
    visit!(state, g, u, v)

Perform a depth first search storing the depth (in `depth`) and low-points
(in `low`) of each vertex.
"""
function visit!(state::Bridges, g::AbstractGraph, u::Integer, v::Integer)
    state.id += 1
    state.depth[v] = state.id
    state.low[v] = state.depth[v]

    @inbounds for w in outneighbors(g, v)
        if state.depth[w] == 0
            visit!(state, g, v, w)

            state.low[v] = min(state.low[v], state.low[w])
            if state.low[w] > state.depth[v]
                edge = v < w ? Edge(v, w) : Edge(w, v)
                push!(state.bridges, edge)
            end

        elseif w != u
            state.low[v] = min(state.low[v], state.depth[w])
        end
    end
end

"""
    bridges(g)

Compute the [bridges](https://en.m.wikipedia.org/wiki/Bridge_(graph_theory))
of a connected graph `g` and return an array containing all bridges, i.e edges
whose deletion increases the number of connected components of the graph.
# Examples
```jldoctest
julia> using LightGraphs

julia> bridges(StarGraph(5))
8-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 1 => 4
 Edge 1 => 5

julia> bridges(PathGraph(5))
8-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 4 => 5
 Edge 3 => 4
 Edge 2 => 3
 Edge 1 => 2
```
"""
function bridges end

@traitfn function bridges(g::AG::(!IsDirected)) where {T<:Integer, AG<:AbstractGraph{T}}
    state = Bridges(g)
    for u in vertices(g)
        if state.depth[u] == 0
            visit!(state, g, u, u)
        end
    end
    return state.bridges
end


text:
module Registry

import ..Pkg, ..Types, ..API
using ..Pkg: depots1
using ..Types: RegistrySpec, Context, Context!
using UUIDs

"""
    Pkg.Registry.add(url::String)
    Pkg.Registry.add(registry::RegistrySpec)

Add new package registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.add("General")
Pkg.Registry.add(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
Pkg.Registry.add(RegistrySpec(url = "https://github.com/JuliaRegistries/General.git"))
```
"""
add(reg::Union{String,RegistrySpec}; kwargs...) = add([reg]; kwargs...)
add(regs::Vector{String}; kwargs...) = add([RegistrySpec(name = name) for name in regs]; kwargs...)
add(regs::Vector{RegistrySpec}; kwargs...) = add(Context(), regs; kwargs...)
function add(ctx::Context, regs::Vector{RegistrySpec}; kwargs...)
    Context!(ctx; kwargs...)
    if isempty(regs)
        Types.clone_default_registries(ctx, only_if_empty = false)
    else
        Types.clone_or_cp_registries(ctx, regs)
    end
end

"""
    Pkg.Registry.rm(registry::String)
    Pkg.Registry.rm(registry::RegistrySpec)

Remove registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.rm("General")
Pkg.Registry.rm(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
rm(reg::Union{String,RegistrySpec}; kwargs...) = rm([reg]; kwargs...)
rm(regs::Vector{String}; kwargs...) = rm([RegistrySpec(name = name) for name in regs]; kwargs...)
rm(regs::Vector{RegistrySpec}; kwargs...) = rm(Context(), regs; kwargs...)
function rm(ctx::Context, regs::Vector{RegistrySpec}; kwargs...)
    Context!(ctx; kwargs...)
    Types.remove_registries(ctx, regs)
end

"""
    Pkg.Registry.update()
    Pkg.Registry.update(registry::RegistrySpec)
    Pkg.Registry.update(registry::Vector{RegistrySpec})

Update registries. If no registries are given, update
all available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.update()
Pkg.Registry.update("General")
Pkg.Registry.update(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
update(reg::Union{String,RegistrySpec}; kwargs...) = update([reg]; kwargs...)
update(regs::Vector{String}; kwargs...) = update([RegistrySpec(name = name) for name in regs]; kwargs...)
update(regs::Vector{RegistrySpec} = RegistrySpec[]; kwargs...) =
    update(Context(), regs; kwargs...)
function update(ctx::Context,
                regs::Vector{RegistrySpec} = RegistrySpec[];
                kwargs...)
    isempty(regs) && (regs = Types.collect_registries(depots1()))
    Context!(ctx; kwargs...)
    Types.update_registries(ctx, regs; force=true)
end

"""
    Pkg.Registry.status()

Display information about available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.status()
```
"""
status(; kwargs...) = status(Context(); kwargs...)
function status(ctx::Context; io::IO=stdout, as_api=false, kwargs...) # TODO split as_api into own function
    Context!(ctx; io=io, kwargs...)
    regs = Types.collect_registries()
    regs = unique(r -> r.uuid, regs; seen=Set{Union{UUID,Nothing}}())
    as_api && return regs
    Types.printpkgstyle(ctx, Symbol("Registry Status"), "")
    if isempty(regs)
        println(ctx.io, "  (no registries found)")
    else
        for reg in regs
            printstyled(ctx.io, " [$(string(reg.uuid)[1:8])]"; color = :light_black)
            print(ctx.io, " $(reg.name)")
            reg.url === nothing || print(ctx.io, " ($(reg.url))")
            println(ctx.io)
        end
    end
end

end # module


text:
arr = [1,2,3]
# arr[0] #> causes a program stop with: ERROR: BoundsError()
# sqrt(-3) #>  causes an ERROR: DomainError: sqrt with -3.0:
# sqrt will only return a complex result if called with a complex argument.
sqrt(complex(-3)) #> 0.0 + 1.7320508075688772im

mutable struct CustomException <: Exception
end

# throw:
codes = ["AO", "ZD", "SG", "EZ"]
# code = "AR"
code = "AO"
if code in codes # in(code, codes)
    println("This is an acceptable code")
else
    throw(DomainError())
end

# error:
function sqroot(n::Int)
    n >= 0 || error("n must be non-negative")
    n == 0 && return 0
    sqrt(n)
end
sqroot(4) #> 2.0
sqroot(0) #> 0.0
# sqroot(-6) #> ERROR: n must be non-negative

# try / catch
a = []
try
    pop!(a)
catch ex
    println(typeof(ex))  #> ArgumentError
    backtrace()
    showerror(stdout, ex) #> ArgumentError: array must be non-empty
end

try
# try this code
catch ex
    if isa(ex, DomainError)
# do this
    elseif isa(ex, BoundsError)
# do this
    end
end

ret =   try
            global a = 4 * 2
        catch ex
        end
println(ret) #> 8

# finally
try
    global f = open("file1.txt", "r") #> IOStream(<file file1.txt>)
    # operate on file f, for example:
    # alternative: open("file1.txt", "r") do f
    k = 0
    while(!eof(f))
        a=readline(f)
        println(a)
        k += 1
    end
    println("\nNumber of lines in file: $k")
catch ex
finally
    close(f)
end


text:
using Trixi
using OrdinaryDiffEq
using Plots
using LinearAlgebra

###############################################################################
# Definition of the 2D scalar "KPP equation"
#
# See: Kurganov, A., Petrova, G., and Popov, B. (2007).
# Adaptive Semidiscrete Central-Upwind Schemes for Nonconvex Hyperbolic Conservation Laws,
# SIAM Journal on Scientific Computing, 29(6), 2381--2401
# DOI: https://doi.org/10.1137/040614189

struct KPPEquation2D <: Trixi.AbstractEquations{2, 1} end

# The KPP flux is F(u) = (sin(u), cos(u))
@inline function Trixi.flux(u, orientation::Integer, ::KPPEquation2D)
    if orientation == 1
        return SVector(sin(u[1]))
    else
        return SVector(cos(u[1]))
    end
end

# Since the KPP problem is a scalar equation, the entropy-conservative flux is uniquely determined
@inline function Trixi.flux_ec(u_ll, u_rr, orientation::Integer, ::KPPEquation2D)
    # The tolerance of 1e-12 is based on experience and somewhat arbitrarily chosen
    if abs(u_ll[1] - u_rr[1]) < 1e-12
        return 0.5 * (flux(u_ll, orientation, KPPEquation2D()) + flux(u_rr, orientation, KPPEquation2D()))
    else
        factor = 1.0 / (u_rr[1] - u_ll[1])
        if orientation == 1
            return SVector(factor*(-cos(u_rr[1]) + cos(u_ll[1])))
        else
            return SVector(factor*(sin(u_rr[1]) - sin(u_ll[1])))
        end
    end
end

# Wavespeeds
@inline wavespeed(::KPPEquation2D) = 1.0
@inline Trixi.max_abs_speeds(u, equation::KPPEquation2D) = (wavespeed(equation), wavespeed(equation))
@inline Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equation::KPPEquation2D) = wavespeed(equation)
@inline Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, equation::KPPEquation2D) = wavespeed(equation) * norm(normal_direction)

# Compute entropy: we use the square entropy
@inline Trixi.entropy(u::Real, ::KPPEquation2D) = 0.5 * u^2
@inline Trixi.entropy(u, ::KPPEquation2D) = entropy(u[1], equation)

# Convert between conservative, primitive, and entropy variables. The conserved quantity "u" is also
# considered the "primitive variable". Since we use the square entropy, "u" is also the entropy
# variable.
@inline Trixi.cons2prim(u, ::KPPEquation2D) = u
@inline Trixi.cons2entropy(u, ::KPPEquation2D) = u

Trixi.varnames(::Any, ::KPPEquation2D) = ("u",)

# Standard KPP test problem with discontinuous initial condition
function initial_condition_kpp(x, t, ::KPPEquation2D)
    if x[1]^2 + x[2]^2 < 1
        return SVector(0.25 * 14.0 * pi)
    else
        return SVector(0.25 * pi)
    end
end

###############################################################################
# semidiscretization of the KPP problem
equation = KPPEquation2D()
surface_flux = flux_lax_friedrichs
volume_flux = flux_ec
# Can also compare solution obtained without using entropy-conservative flux. This will not converge
# to the correct (entropy) solution!
# volume_flux = flux_central
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
shock_indicator = IndicatorHennemannGassner(equation, basis,
                                            alpha_max=0.5,
                                            alpha_min=0.001,
                                            alpha_smooth=true,
                                            variable=first)
volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

###############################################################################
# Set up the tree mesh (initially a Cartesian grid of [-2,2]^2)
coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)

mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6,
                periodicity=true,
                n_cells_max=500_000)

###############################################################################
# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition_kpp, solver)

###############################################################################
# Set up adaptive mesh refinement
amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=1.0,
                                          alpha_min=0.0001,
                                          alpha_smooth=false,
                                          variable=first)

max_refinement_level = 8

amr_controller = ControllerThreeLevelCombined(semi, amr_indicator, shock_indicator,
                                              base_level=2,
                                              med_level=0, med_threshold=0.0003,
                                              max_level=max_refinement_level, max_threshold=0.003,
                                              max_threshold_secondary=shock_indicator.alpha_max)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=1,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

###############################################################################
# ODE solvers, callbacks etc.
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_callback = AnalysisCallback(semi, interval=200)

alive_callback = AliveCallback(analysis_interval=200)

save_solution = SaveSolutionCallback(interval=200,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2cons)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution, amr_callback)

###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(), save_everystep=false, callback=callbacks)

summary_callback() # Print the timer summary


text:
###
### An abstract biological sequence type.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

abstract type BioSequence{A<:Alphabet} end

# Aliases and shorthands for describing subsets of the BioSequence type...
const NucleotideSeq = BioSequence{<:NucleicAcidAlphabet}
const AminoAcidSeq = BioSequence{AminoAcidAlphabet}

###
### Required traits and methods
###

# Base.length must be defined for every T<:BioSequence.

# As must the following...

"Get the length of a biological sequence."
@inline function Base.length(seq::BioSequence)
    error(
        string(
            "Base.length has not been defined for BioSequence type: ",
            typeof(seq),
            ". Any compatible concrete BioSequence subtype must have this method implemented."
        )
    )
end

"""
Return the data member of `seq` that stores the encoded sequence data.
"""
@inline function encoded_data(seq::BioSequence)
    error(
        string(
            "encoded_data has not been defined for BioSequence type: ",
            typeof(seq),
            ". Any compatible concrete BioSequence subtype must have this method implemented."
        )
    )
end

###
### Provided traits and methods
###

# These traits and methods are defined automatically for any subtype of BioSequence{A}.
# They may be overloaded for your concrete BioSequence sub-type if it is nessecery.

"Get the vector of bits storing a sequences packed encoded elements."
@inline encoded_data_type(seq::BioSequence) = typeof(encoded_data(seq))

"Get the element type of the vector of bits storing a sequences packed encoded elements."
@inline encoded_data_eltype(seq::BioSequence) = eltype(encoded_data_type(seq))

"""
Return the `Alpahbet` defining the possible biological symbols
and their encoding for a given biological sequence.
"""
@inline function Alphabet(::Type{<:BioSequence{A}}) where {A <: Alphabet}
    return A()
end

"Return the `Alpahbet` type that defines the biological symbols allowed for `seq`."
@inline function Alphabet(seq::BioSequence)
    return Alphabet(typeof(seq))
end

BioSymbols.alphabet(::Type{BioSequence{A}}) where {A<:Alphabet} = alphabet(A)

BitsPerSymbol(seq::BioSequence) = BitsPerSymbol(Alphabet(seq))

"Get the number of bits each symbol packed into a BioSequence uses, as an integer value."
bits_per_symbol(seq::BioSequence) = bits_per_symbol(Alphabet(seq))

# The generic functions for any BioSequence...
include("indexing.jl")
include("conversion.jl")
include("predicates.jl")
include("find.jl")
include("printing.jl")
include("transformations.jl")
include("counting.jl")


text:
function depthfirstsearch!(g::AbstractGraph, v::Vertex)
    #previsit(g, v)
    setmark!(g, v, 1)
    let w = first(g, v)
        while w <= countvertices(g)
            if getmark(g, w) == 0
                #process(g, v, w)
                println(v, w)
                depthfirstsearch!(g, w)
            end
            w = next(g, v, w)
        end
    end
    #postvisit(g, v)
end


text:
module SEP

import Base: collect, broadcast!, -, show
export Background, background, rms

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("SEP not properly installed. Please run Pkg.build(\"SEP\")")
end

# Julia Type -> numeric type code from sep.h
const SEPMaskType = Union{Bool, UInt8, Cint, Cfloat, Cdouble}
const PixType = Union{UInt8, Cint, Cfloat, Cdouble}
sep_typecode(::Type{Bool}) = Cint(11)
sep_typecode(::Type{UInt8}) = Cint(11)
sep_typecode(::Type{Cint}) = Cint(31)
sep_typecode(::Type{Cfloat}) = Cint(42)
sep_typecode(::Type{Cdouble}) = Cint(82)

# definitions from sep.h
const SEP_NOISE_NONE = Cshort(0)
const SEP_NOISE_STDDEV = Cshort(1)
const SEP_NOISE_VAR = Cshort(2)

const SEP_THRESH_REL = Cint(0)
const SEP_THRESH_ABS = Cint(1)

const SEP_FILTER_CONV = Cint(0)
const SEP_FILTER_MATCHED = Cint(1)

function sep_assert_ok(status::Cint)
    if status != 0
        msg = Vector{UInt8}(61)
        ccall((:sep_get_errmsg, libsep), Void, (Int32,Ptr{UInt8}), status, msg)
        msg[end] = 0  # ensure NULL-termination, just in case.
        error(unsafe_string(pointer(msg)))
    end
end

# internal use only: mirrors `sep_image` struct in sep.h
struct sep_image
data::Ptr{Void}
noise::Ptr{Void}
mask::Ptr{Void}
dtype::Cint
ndtype::Cint
mdtype::Cint
w::Cint
h::Cint
noiseval::Cdouble
noise_type::Cshort
gain::Cdouble
maskthresh::Cdouble
end

function sep_image(data::Array{T, 2};
                   noise=nothing, mask=nothing, noise_type=:stddev,
                   gain=0.0, mask_thresh=0.0) where {T<:PixType}
    sz = size(data)

    # data is required
    data_ptr = Ptr{Void}(pointer(data))
    dtype = sep_typecode(T)

    # mask options
    mask_ptr = C_NULL
    mdtype = Cint(0)
    if mask !== nothing
        isa(mask, Matrix) || error("mask must be a 2-d array")
        size(mask) == sz || error("data and mask must be same size")
        mask_ptr = Ptr{Void}(pointer(mask))
        mdtype = sep_typecode(eltype(mask))
    end

    # noise options
    ndtype = Cint(0)
    noise_ptr = C_NULL
    noiseval = 0.0
    noise_typecode = SEP_NOISE_NONE
    if noise !== nothing
        if isa(noise, Matrix)
            size(noise) == sz || error("noise array must be same size as data")
            noise_ptr = sep_typecode(eltype(noise))
            ndtype = sep_typecode(eltype(ndtype))
        elseif isa(noise, Number)
            noiseval = Cdouble(noise)
        else
            error("noise must be array or number")
        end
        noise_typecode = ((noise_type == :stddev) ? SEP_NOISE_STDDEV :
                          (noise_type == :var) ? SEP_NOISE_VAR :
                          error("noise_type must be :stddev or :var"))
    end

    return sep_image(data_ptr, noise_ptr, mask_ptr,
                     dtype, ndtype, mdtype,
                     sz[1], sz[2],
                     noiseval, noise_typecode,
                     Cdouble(gain),
                     Cdouble(maskthresh))
end

# ---------------------------------------------------------------------------
# Background functions

"""
    Background

"""
mutable struct Background
    ptr::Ptr{Void}
    data_size::Tuple{Int, Int}
end

free!(bkg::Background) = ccall((:sep_bkg_free, libsep), Void, (Ptr{Void},),
                               bkg.ptr)

global_mean(bkg::Background) = ccall((:sep_bkg_global, libsep), Cfloat,
                                     (Ptr{Void},), bkg.ptr)

global_rms(bkg::Background) = ccall((:sep_bkg_globalrms, libsep), Cfloat,
                                    (Ptr{Void},), bkg.ptr)

function show(io::IO, bkg::Background)
    print(io, "Background $(bkg.data_size[1])×$(bkg.data_size[2])\n")
    print(io, " - global mean: $(global_mean(bkg))\n")
    print(io, " - global rms : $(global_rms(bkg))\n")
end

function background(data::Array{T, 2}; mask=nothing, boxsize=(64, 64),
                    filtersize=(3, 3), filterthresh=0.0, mask_thresh=0) where {T<:PixType}
    im = sep_image(data; mask=mask, mask_thresh=mask_thresh)
    result = Ref{Ptr{Void}}(C_NULL)
    status = ccall((:sep_background, libsep), Cint,
                   (Ptr{sep_image}, Cint, Cint, Cint, Cint, Cdouble,
                    Ref{Ptr{Void}}),
                   &im, boxsize[1], boxsize[2], filtersize[1], filtersize[2],
                   filterthresh, result)
    sep_assert_ok(status)
    bkg = Background(result[], sz)
    finalizer(bkg, free!)
    return bkg
end

function collect(::Type{T}, bkg::Background) where {T<:PixType}
    result = Array{T}(bkg.data_size)
    status = ccall((:sep_bkg_array, libsep), Cint, (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, result, sep_typecode(T))
    sep_assert_ok(status)
    return result
end
# default collection type is Float32, because that's what's natively stored in
# background.
# TODO: make default the input array type in `background` instead?
collect(bkg::Background) = collect(Float32, bkg)

"""
    rms(bkg)
    rms(T, bkg)

Return an array of the standard deviation of the background. The result is
the size of the original image and is of type `T`, if given.
"""
function rms(::Type{T}, bkg::Background) where {T<:PixType}
    result = Array{T}(bkg.data_size)
    status = ccall((:sep_bkg_rmsarray, libsep), Cint,
                   (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, result, sep_typecode(T))
    sep_assert_ok(status)
    return result
end
rms(bkg::Background) = rms(Float32, bkg)

# In-place background subtraction: A .-= bkg
function broadcast!(-, A::Array{T, 2},  ::Array{T, 2}, bkg::Background) where {T<:PixType}
    if size(A) != bkg.data_size
         throw(DimensionMismatch("dimensions must match"))
    end
    status = ccall((:sep_bkg_subarray, libsep), Cint,
                   (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, A, sep_typecode(T))
    sep_assert_ok(status)
end

function (-)(A::Array{T, 2}, bkg::Background) where {T<:PixType}
    B = copy(A)
    B .-= bkg
    return B
end

# ---------------------------------------------------------------------------
# Source Extraction

# Mirror of C struct
struct sep_catalog
    nobj::Cint                 # number of objects (length of all arrays)
    thresh::Ptr{Cfloat}              # threshold (ADU)
    npix::Ptr{Cint}              # # pixels extracted (size of pix array)
    tnpix::Ptr{Cint}                # # pixels above thresh (unconvolved)
    xmin::Ptr{Cint}
    xmax::Ptr{Cint}
    ymin::Ptr{Cint}
    ymax::Ptr{Cint}
    x::Ptr{Cdouble}              # barycenter (first moments)
    y::Ptr{Cdouble}
    x2::Ptr{Cdouble}     # second moments
    y2::Ptr{Cdouble}
    xy::Ptr{Cdouble}
    errx2::Ptr{Cdouble}  # second moment errors
    erry2::Ptr{Cdouble}
    errxy::Ptr{Cdouble}
    a::Ptr{Cfloat}       # ellipse parameters
    b::Ptr{Cfloat}
    theta::Ptr{Cfloat}
    cxx::Ptr{Cfloat}        # alternative ellipse parameters
    cyy::Ptr{Cfloat}
    cxy::Ptr{Cfloat}
    cflux::Ptr{Cfloat}   # total flux of pixels (convolved)
    flux::Ptr{Cfloat}    # total flux of pixels (unconvolved)
    cpeak::Ptr{Cfloat}   # peak pixel flux (convolved)
    peak::Ptr{Cfloat}    # peak pixel flux (unconvolved)
    xcpeak::Ptr{Cint}    # x, y coords of peak (convolved) pixel
    ycpeak::Ptr{Cint}
    xpeak::Ptr{Cint}     # x, y coords of peak (unconvolved) pixel
    ypeak::Ptr{Cint}
    flag::Ptr{Cshort}    # extraction flags
    pix::Ptr{Ptr{Cint}}  # array giving indicies of object's pixels in
                         # image (linearly indexed). Length is `npix`.
                         # (pointer to within the `objectspix` buffer)
    objectspix::Ptr{Cint}  # buffer holding pixel indicies for all objects
end

function extract(data::Array{T, 2}, thresh::Real;
                 noise=nothing, mask=nothing, noise_type=:stddev,
                 thresh_type=:relative, minarea=5,
                 filter_kernel=Float32[1 2 1; 2 4 2; 1 2 1],
                 filter_type=:matched,
                 deblend_nthresh=32, deblend_cont=0.005,
                 clean=true, clean_param=1.0,
                 gain=0.0, mask_thresh=0) where {T}

    im = sep_image(data; noise=noise, mask=mask, noise_type=noise_type,
                   gain=gain, mask_thresh=mask_thresh)

    thresh_typecode = ((thresh_type == :relative)? SEP_THRESH_REL :
                       (thresh_type == :absolute)? SEP_THRESH_ABS :
                       error("thresh_type must be :relative or :absolute"))

    filter_typecode = ((filter_type == :matched)? SEP_FILTER_MATCHED :
                       (filter_type == :convolution)? SEP_FILTER_CONV :
                       error("filter_type must be :matched or :convolution"))

    # convert filter kernel to Cfloat array
    filter_kernel_cfloat = convert(Array{Cfloat, 2}, filter_kernel)
    filter_size = size(filter_kernel_cfloat)

    result = Ref{Ptr{sep_catalog}}(C_NULL)
    status = ccall((:sep_extract, libsep), Cint,
                   (Ptr{sep_image},
                    Cfloat, Cint, # thresh, thresh_type
                    Cint,  # minarea
                    Ptr{Cfloat}, Cint, Cint,  # conv, convw, convh
                    Cint, # filter_type
                    Cint, Cdouble,  # deblend_nthresh, deblend_cont
                    Cint, Cdouble,  # clean_flag, clean_param
                    Ref{Ptr{sep_catalog}}),
                   &im,
                   thresh, thresh_typecode,
                   minarea,
                   filter_kernel_cfloat, filter_size[1], filter_size[2],
                   filter_typecode,
                   deblend_nthresh, deblend_cont,
                   clean_flag, clean_param,
                   result)
    sep_assert_ok(status)

    # parse result
    catalog_ptr = result[]
    # ...

    # free result allocated in sep_extract
    ccall((:sep_catalog_free, libsep), Void, (Ptr{sep_catalog},), catalog_ptr)

    # todo: return actual result
    nothing
end

end # module


text:
struct TSNode{T}
    head::T
    tail::Union{TSNode{T},Nothing}
end

const TSList{T} = Union{TSNode{T},Nothing}

struct TreiberStack{T,Ref<:Reagents.Ref{TSList{T}}}
    head::Ref
end

TreiberStack{T}() where {T} = TreiberStack(Reagents.Ref{TSList{T}}(nothing))

Base.eltype(::Type{<:TreiberStack{T}}) where {T} = T
Base.isempty(stack::TreiberStack) = stack.head[] === nothing

pushing(stack::TreiberStack{T}) where {T} =
    Update((xs, x) -> (TSNode{T}(x, xs), nothing), stack.head)

trypopping(stack::TreiberStack) =
    Update(stack.head) do xs, _ignored
        if xs === nothing
            return (nothing, nothing)
        else
            return (xs.tail, Some(xs.head))
        end
    end

Base.push!(stack::TreiberStack, x) = pushing(stack)(convert(eltype(stack), x))
maybepop_nowait!(stack::TreiberStack) = trypopping(stack)()

Base.IteratorSize(::Type{<:TreiberStack}) = Base.SizeUnknown()
function Base.iterate(stack::TreiberStack, list = stack.head[])
    node = @something(list, return nothing)
    return (node.head, node.tail)
end


text:
using Literate

function comments2script(content)
    content = replace(content, "#~" => "# ")
    return content
end

function comments2md(content)
    content = replace(content, "#~" => "# ")
    return content
end

function replace_includes(str)
    included = ["lit_tutorial_inputs.jl",]
    path = "docs/src/literate/tutorial_lit/"

    for ex in included
        content = read(path*ex, String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    str = comments2md(str)

    return str
end

function replace_includes_advanced(str)
    included = ["lit_advanced_tutorial_inputs.jl",]
    path = "docs/src/literate/advanced_tutorial_lit/"

    for ex in included
        content = read(path*ex, String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    str = comments2md(str)

    return str
end

function replace_includes_et(str)
    included = ["lit_inputs.jl",]
    path = "docs/src/literate/empty_template_lit/"

    for ex in included
        content = read(path*ex, String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    str = comments2md(str)

    return str
end

#====== Tutorial ======================================================#
Literate.script("docs/src/literate/tutorial_lit/lit_runTutorial.jl",
    "examples/tutorial/", name="runTutorial",
    keep_comments=true, postprocess=comments2script;)

Literate.script("docs/src/literate/tutorial_lit/lit_tutorial_inputs.jl",
    "examples/tutorial/", name="tutorial_inputs",
    keep_comments=true, postprocess=comments2script;)

Literate.markdown("docs/src/literate/tutorial_lit/lit_runTutorial.jl",
    "docs/src/", preprocess = replace_includes,
    name="tutorial", documenter=false)

Literate.notebook("docs/src/literate/tutorial_lit/lit_runTutorial.jl",
"examples/notebooks/", preprocess = replace_includes,
name="Tutorial", documenter=false)

# Plotting
Literate.script("docs/src/literate/tutorial_lit/lit_plotting.jl",
    "examples/tutorial/", name="plotting",
    keep_comments=true, postprocess=comments2script;)

Literate.markdown("docs/src/literate/tutorial_lit/lit_plotting.jl",
    "docs/src/",  name="plotting", preprocess=comments2script, documenter=false)

#====== Advanced Tutorial ======================================================#
Literate.script("docs/src/literate/advanced_tutorial_lit/lit_runAdvancedTutorial.jl",
    "examples/advanced_tutorial/", name="runAdvancedTutorial",
    keep_comments=true, postprocess=comments2script;)

Literate.script("docs/src/literate/advanced_tutorial_lit/lit_advanced_tutorial_inputs.jl",
    "examples/advanced_tutorial/", name="advanced_inputs",
    keep_comments=true, postprocess=comments2script;)

Literate.markdown("docs/src/literate/advanced_tutorial_lit/lit_runAdvancedTutorial.jl",
    "docs/src/", preprocess = replace_includes_advanced,
    name="advanced_tutorial", documenter=false)

Literate.notebook("docs/src/literate/advanced_tutorial_lit/lit_runAdvancedTutorial.jl",
"examples/notebooks/", preprocess = replace_includes_advanced,
name="AdvancedTutorial", documenter=false)

#====== BLUE ======================================================#
Literate.script("docs/src/literate/BLUE_lit/lit_runBLUE.jl",
    "examples/BLUE/", name="runBLUE",
    keep_comments=true, postprocess=comments2script;)

Literate.markdown("docs/src/literate/BLUE_lit/lit_runBLUE.jl",
    "docs/src/", name="BLUE", preprocess=comments2md, documenter=false)

Literate.notebook("docs/src/literate/BLUE_lit/lit_runBLUE.jl",
"examples/notebooks/", name="BLUE", preprocess=comments2md)

#====== Empty Template ======================================================#
Literate.script("docs/src/literate/empty_template_lit/lit_runTemplate.jl",
    "examples/empty_template/", name="runEFTfitter",
    keep_comments=true, postprocess=comments2script;)

Literate.script("docs/src/literate/empty_template_lit/lit_inputs.jl",
    "examples/empty_template/", name="inputs",
    keep_comments=true, postprocess=comments2script;)

Literate.notebook("docs/src/literate/empty_template_lit/lit_runTemplate.jl",
"examples/notebooks/", preprocess = replace_includes_et,
name="EmptyTemplate", documenter=false)


text:
using IntervalArithmetic
using RevisedAffineArithmetic
using StaticArrays
using Test

const RAA = RevisedAffineArithmetic

include("test_affineform.jl")
include("test_arithmetic.jl")
include("test_conversions.jl")
include("test_vars.jl")


text:
import Random: Sampler
using Random
function direct_sample!(rng::AbstractRNG, a::UnitRange, x::AbstractArray)
    s = Sampler(rng, 1:length(a))
    b = a[1] - 1
    if b == 0
        for i = 1:length(x)
            @inbounds x[i] = rand(rng, s)
        end
    else
        for i = 1:length(x)
            @inbounds x[i] = b + rand(rng, s)
        end
    end
    return x
end
function sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray;
                 replace::Bool=true, ordered::Bool=false)
    n = length(a)
    k = length(x)
    k == 0 && return x

    if replace  # with replacement
        if ordered
            sort!(direct_sample!(rng, a, x))
        else
            direct_sample!(rng, a, x)
        end

    else  # without replacement
        k <= n || error("Cannot draw more samples without replacement.")

        if ordered
            if n > 10 * k * k
                seqsample_c!(rng, a, x)
            else
                seqsample_a!(rng, a, x)
            end
        else
            if k == 1
                @inbounds x[1] = sample(rng, a)
            elseif k == 2
                @inbounds (x[1], x[2]) = samplepair(rng, a)
            elseif n < k * 24
                fisher_yates_sample!(rng, a, x)
            else
                self_avoid_sample!(rng, a, x)
            end
        end
    end
    return x
end
direct_sample!(a::UnitRange, x::AbstractArray) = direct_sample!(Random.GLOBAL_RNG, a, x)
function sample(rng::AbstractRNG, a::AbstractArray{T}, n::Integer;
                replace::Bool=true, ordered::Bool=false) where T
    sample!(rng, a, Vector{T}(undef, n); replace=replace, ordered=ordered)
end
sample(a::AbstractArray, n::Integer; replace::Bool=true, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, n; replace=replace, ordered=ordered)

    sample(rng::AbstractRNG, a) = a[rand(rng, 1:length(a))]
    sample(a::AbstractArray) = sample(Random.GLOBAL_RNG, a)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module FieldValues

using ...VersionWeights
importall .....Base.Operators

export FieldValue, Field, validmax, secondmax

# FieldValue is a numeric type which helps dealing with
# infinities. It holds 5 numbers l0,l1,l2,l3,l4. It can
# be interpreted as a polynomial
#  x = a^4 * l0 + a^3 * l1 + a^2 + l2 + a^1 * l3 + l4
# where a -> Inf
# The levels are used as such:
#  l0 : for hard constraints (dependencies and requirements)
#  l1 : for favoring higher versions of the explicitly required
#       packages
#  l2 : for favoring higher versions of all other packages (and
#       favoring uninstallation of non-needed packages)
#  l3 : for favoring dependants over dependencies
#  l4 : for symmetry-breaking random noise
#
immutable FieldValue
    l0::Int
    l1::VersionWeight
    l2::VersionWeight
    l3::Int
    l4::Int128
end
FieldValue(l0::Int,l1::VersionWeight,l2::VersionWeight,l3::Int) = FieldValue(l0, l1, l2, l3, Int128(0))
FieldValue(l0::Int,l1::VersionWeight,l2::VersionWeight) = FieldValue(l0, l1, l2, 0)
FieldValue(l0::Int,l1::VersionWeight) = FieldValue(l0, l1, zero(VersionWeight))
FieldValue(l0::Int) = FieldValue(l0, zero(VersionWeight))
FieldValue() = FieldValue(0)

typealias Field Vector{FieldValue}

Base.zero(::Type{FieldValue}) = FieldValue()

Base.typemin(::Type{FieldValue}) = (x=typemin(Int); y=typemin(VersionWeight); FieldValue(x,y,y,x,typemin(Int128)))

Base.:-(a::FieldValue, b::FieldValue) = FieldValue(a.l0-b.l0, a.l1-b.l1, a.l2-b.l2, a.l3-b.l3, a.l4-b.l4)
Base.:+(a::FieldValue, b::FieldValue) = FieldValue(a.l0+b.l0, a.l1+b.l1, a.l2+b.l2, a.l3+b.l3, a.l4+b.l4)

function Base.isless(a::FieldValue, b::FieldValue)
    a.l0 < b.l0 && return true
    a.l0 > b.l0 && return false
    c = cmp(a.l1, b.l1)
    c < 0 && return true
    c > 0 && return false
    c = cmp(a.l2, b.l2)
    c < 0 && return true
    c > 0 && return false
    a.l3 < b.l3 && return true
    a.l3 > b.l3 && return false
    a.l4 < b.l4 && return true
    return false
end

==(a::FieldValue, b::FieldValue) =
    a.l0 == b.l0 && a.l1 == b.l1 && a.l2 == b.l2 && a.l3 == b.l3 && a.l4 == b.l4

Base.abs(a::FieldValue) = FieldValue(abs(a.l0), abs(a.l1), abs(a.l2), abs(a.l3), abs(a.l4))

# if the maximum field has l0 < 0, it means that
# some hard constraint is being violated
validmax(a::FieldValue) = a.l0 >= 0

# like usual indmax, but favors the highest indices
# in case of a tie
function Base.indmax(f::Field)
    m = typemin(FieldValue)
    mi = 0
    for j = length(f):-1:1
        if f[j] > m
            m = f[j]
            mi = j
        end
    end
    @assert mi != 0
    return mi
end

# secondmax returns the (normalized) value of the second maximum in a
# field. It's used to determine the most polarized field.
function secondmax(f::Field)
    m = typemin(FieldValue)
    m2 = typemin(FieldValue)
    for i = 1:length(f)
        a = f[i]
        if a > m
            m2 = m
            m = a
        elseif a > m2
            m2 = a
        end
    end
    return m2 - m
end

end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia