text:
using Base.Meta
using JuMP: _valid_model, _error_if_cannot_register, object_dictionary, variable_type
using JuMP.Containers

# Parse raw input to define the upper bound for an interval set
function _parse_one_operator_parameter(
    _error::Function, infoexpr::_ParameterInfoExpr, ::Union{Val{:<=}, Val{:≤}},
    upper)
    JuMP._set_upper_bound_or_error(_error, infoexpr, upper)
end

# Parse raw input to define the lower bound for an interval set
function _parse_one_operator_parameter(
    _error::Function, infoexpr::_ParameterInfoExpr, ::Union{Val{:>=}, Val{:≥}},
    lower)
    JuMP._set_lower_bound_or_error(_error, infoexpr, lower)
end

# Parse raw input to define the distribution for a dist set
function _parse_one_operator_parameter(
    _error::Function, infoexpr::_ParameterInfoExpr, ::Union{Val{:in}, Val{:∈}},
    value)
    # check if interval set
    if isexpr(value.args[1], :vect)
        JuMP._set_lower_bound_or_error(_error, infoexpr,
                                  JuMP._esc_non_constant(value.args[1].args[1]))
        JuMP._set_upper_bound_or_error(_error, infoexpr,
                                  JuMP._esc_non_constant(value.args[1].args[2]))
    else
        _set_or_error(_error, infoexpr, value)
    end
end

# Default function for unrecognized operators
function _parse_one_operator_parameter(
    _error::Function, infoexpr::_ParameterInfoExpr, ::Val{S}, value) where S
    _error("Unknown sense $S.")
end

# In that case we assume the variable is the lhs.
function _parse_parameter(_error::Function, infoexpr::_ParameterInfoExpr,
                          sense::Symbol, var, value)
    _parse_one_operator_parameter(_error, infoexpr, Val(sense),
                                JuMP._esc_non_constant(value))
    return var
end

# If the lhs is a number and not the rhs, we can deduce that the rhs is
# the variable.
function _parse_parameter(_error::Function, infoexpr::_ParameterInfoExpr,
                          sense::Symbol, value::Number, var)
    _parse_one_operator_parameter(_error, infoexpr,
                                  JuMP.reverse_sense(Val(sense)),
                                  JuMP._esc_non_constant(value))
    return var
end

# Parse raw input to define the upper and lower bounds for an interval set
function _parse_ternary_parameter(_error::Function, infoexpr::_ParameterInfoExpr,
                                  ::Union{Val{:<=}, Val{:≤}}, lower,
                                  ::Union{Val{:<=}, Val{:≤}}, upper)
    JuMP._set_lower_bound_or_error(_error, infoexpr, lower)
    JuMP._set_upper_bound_or_error(_error, infoexpr, upper)
end

# Parse raw input to define the upper and lower bounds for an interval set
function _parse_ternary_parameter(_error::Function, infoexpr::_ParameterInfoExpr,
                                  ::Union{Val{:>=}, Val{:≥}}, upper,
                                  ::Union{Val{:>=}, Val{:≥}}, lower)
    _parse_ternary_parameter(_error, infoexpr, Val(:≤), lower, Val(:≤), upper)
end

# Default for unrecognized operators
function _parse_ternary_parameter(_error::Function, infoexpr::_ParameterInfoExpr,
                                  ::Val, lvalue, ::Val, rvalue)
    _error("Use the form lb <= ... <= ub.")
end

# Interpret a param_expr from macro input
function _parse_parameter(_error::Function, infoexpr::_ParameterInfoExpr, lvalue,
                          lsign::Symbol, var, rsign::Symbol, rvalue)
    # lvalue lsign var rsign rvalue
    _parse_ternary_parameter(_error, infoexpr, Val(lsign),
                             JuMP._esc_non_constant(lvalue), Val(rsign),
                             JuMP._esc_non_constant(rvalue))
    return var
end

# Extend to increase the param group id
function _assert_valid_model_call(m, macrocode)
    # assumes m is already escaped
    quote
        JuMP._valid_model($m, $(quot(m.args[1])))
        $(m).next_param_id += 1
        $macrocode
    end
end

"""
    @infinite_parameter(model, kw_args...)

Add an *anonymous* infinite parameter to the model `model` described by the
keyword arguments `kw_args` and returns the parameter reference.

```julia
@infinite_parameter(model, expr, kw_args...)
```

Add a parameter to the model `model` described by the expression `expr`, the
positional arguments `args` and the keyword arguments `kw_args`. (note that in
the following the symbol `<=` can be used instead of `≤`, the symbol `>=`can
be used instead of `≥`, and the symbo `in` can be used instead of `∈`) The
expression `expr` can be of the form:
- `paramexpr` creating parameters described by `paramexpr`.
- `lb ≤ paramexpr ≤ ub` creating parameters described by `paramexpr` characterized
   by a continuous interval set with lower bound `lb` and upper bound `ub`.
- `paramexpr ∈ [lb, ub]` creating parameters described by `paramexpr` characterized
   by a continuous interval set with lower bound `lb` and upper bound `ub`.
- `paramexpr ∈ dist` creating parameters described by `paramexpr` characterized
   by the `Distributions.jl` distribution object `dist`.

The expression `paramexpr` can be of the form:
- `paramname` creating a scalar parameter of name `paramname`
- `paramname[...]` or `[...]` creating a container of parameters

The recognized keyword arguments in `kw_args` are the following:
- `base_name`: Sets the name prefix used to generate parameter names. It
  corresponds to the parameter name for scalar parameter, otherwise, the
  parameter names are set to `base_name[...]` for each index `...` of the axes
  `axes`.
- `lower_bound`: Sets the value of the parameter lower bound for an interval set.
- `upper_bound`: Sets the value of the parameter upper bound for an interval set.
- `set`: The `InfiniteSet` characterizing the parameters see [`IntervalSet`](@ref)
   and [`DistributionSet`](@ref).
- `distribution`: Sets the `Distributions.jl` distribution object that characterizes
  the parameters.
- `supports`: Sets the support points for the parameters.
- `num_supports`: Specifies the number of supports to be automatically generated.
                  Note that `supports` takes precedence. Defaults to 50.
- `sig_figs`: Specifies the number of significant digits that should be used
              in automatic support generation. Defaults to 5.
- `independent`: Specifies if the each parameter is independent from each other
  or not. Defaults to false.
- `container`: Specify the container type. Defaults to automatic

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP, Distributions; m = InfiniteModel())
julia> @infinite_parameter(m, 0 <= x <= 1)
x

julia> @infinite_parameter(m, y[i = 1:2] in Normal(), num_supports = 10)
2-element Array{ParameterRef,1}:
 y[1]
 y[2]

julia> z = @infinite_parameter(m, ["a", "b"], distribution = Uniform(), independent = true)
2-dimensional DenseAxisArray{ParameterRef,2,...} with index sets:
    Dimension 1, "a"
    Dimension 2, "b"
And data, a 1×1 Array{ParameterRef,2}:
 noname
```
"""
macro infinite_parameter(model, args...)
    esc_model = esc(model)

    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)

    # check to see if error_args are provided define error function
    error_args = filter(kw -> kw.args[1] == :error_args, kw_args)
    if length(error_args) != 0
        new_args = error_args[1].args[2]
        macro_name = :finite_parameter
    else
        new_args = args
        macro_name = :infinite_parameter
    end
    _error(str...) = JuMP._macro_error(macro_name, (model, new_args...),
                                       str...)

    # if there is only a single non-keyword argument, this is an anonymous
    # variable spec and the one non-kwarg is the model
    if length(extra) == 0
        x = gensym()
        anon_singleton = true
    elseif length(extra) == 1
        x = popfirst!(extra)
        anon_singleton = false
    else
        x = popfirst!(extra)
        arg = popfirst!(extra)
        _error("Unrecognized argument $arg provided.")
    end

    info_kw_args = filter(_is_set_keyword, kw_args)
    extra_kw_args = filter(kw -> kw.args[1] != :base_name && !InfiniteOpt._is_set_keyword(kw) && kw.args[1] != :error_args, kw_args)
    base_name_kw_args = filter(kw -> kw.args[1] == :base_name, kw_args)
    infoexpr = InfiniteOpt._ParameterInfoExpr(; JuMP._keywordify.(info_kw_args)...)

    # There are five cases to consider:
    # x                                         | type of x | x.head
    # ------------------------------------------+-----------+------------
    # param                                       | Symbol    | NA
    # param[1:2]                                  | Expr      | :ref
    # param <= ub or param[1:2] <= ub             | Expr      | :call
    # param in [lb, ub]                           | Expr      | :call
    # lb <= param <= ub or lb <= param[1:2] <= ub | Expr      | :comparison
    # In the three last cases, we call parse_variable
    explicit_comparison = isexpr(x, :comparison) || isexpr(x, :call)
    if explicit_comparison
        param = InfiniteOpt._parse_parameter(_error, infoexpr, x.args...)
    else
        param = x
    end

    anonvar = isexpr(param, :vect) || isexpr(param, :vcat) || anon_singleton
    anonvar && explicit_comparison && _error("Cannot use explicit bounds via " *
                                             ">=, <= with an anonymous parameter")
    parameter = gensym()
    name = JuMPC._get_name(param)
    if isempty(base_name_kw_args)
        base_name = anonvar ? "" : string(name)
    else
        base_name = esc(base_name_kw_args[1].args[2])
    end

    if !isa(name, Symbol) && !anonvar
        _error("Expression $name should not be used as a parameter name. Use " *
               "the \"anonymous\" syntax $name = @infinite_parameter(model, " *
               "...) instead.")
    end

    set = InfiniteOpt._constructor_set(_error, infoexpr)
    if isa(param, Symbol)
        # Easy case - a single variable
        buildcall = :( build_parameter($_error, $set, 1, $(extra...)) )
        JuMP._add_kw_args(buildcall, extra_kw_args)
        parametercall = :( add_parameter($esc_model, $buildcall, $base_name,
                                         macro_call = true) )
        creationcode = :($parameter = $parametercall)
    else
        # isa(param, Expr) || _error("Expected $param to be a parameter name") --> not needed... I think
        # We now build the code to generate the variables (and possibly the
        # SparseAxisArray to contain them)
        idxvars, indices = JuMPC._build_ref_sets(_error, param)
        buildcall = :( build_parameter($_error, $set, length(collect($indices)),
                                       $(extra...)) )
        JuMP._add_kw_args(buildcall, extra_kw_args)
        parametercall = :( add_parameter($esc_model, $buildcall,
                                         $(JuMP._name_call(base_name, idxvars)),
                                         macro_call = true, multi_dim = true) )
        creationcode = JuMPC.container_code(idxvars, indices, parametercall,
                                             requestedcontainer)
    end
    if anonvar
        # Anonymous variable, no need to register it in the model-level
        # dictionary nor to assign it to a variable in the user scope.
        # We simply return the variable
        macro_code = creationcode
    else
        # We register the variable reference to its name and
        # we assign it to a variable in the local scope of this name
        macro_code = JuMP._macro_assign_and_return(creationcode, parameter, name,
                                                   model_for_registering = esc_model)
    end
    return _assert_valid_model_call(esc_model, macro_code)
end

"""
    @finite_parameter(model::InfiniteModel, value)

Define and add an anonymous finite parameter to `model` and return its
parameter reference. Its value is equal to `value`.

```julia
    @finite_parameter(model::InfiniteModel, param_expr, value_expr)
```
Define and add a finite parameter(s) to `model` and return appropriate parameter
reference(s). The parameter(s) has/have value(s) indicated by the `value_expr`.
The expression `param_expr` can be of the form:
- `paramname` creating a scalar parameter of name `paramname`
- `paramname[...]` or `[...]` creating a container of parameters

The expression `value_expr` simply expresses the value of the parameter(s). This
is typically a number but could be an array indexed using an index defined in
`param_expr`.

The recognized keyword arguments in `kw_args` are the following:
- `base_name`: Sets the name prefix used to generate parameter names. It
  corresponds to the parameter name for scalar parameter, otherwise, the
  parameter names are set to `base_name[...]` for each index `...` of the axes
  `axes`.
- `container`: Specify the container type.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> par = @finite_parameter(model, 2)
noname

julia> vals = [3, 2];

julia> pars = @finite_parameter(model, [i = 1:2], vals[i], base_name = "par")
2-element Array{ParameterRef,1}:
 par[1]
 par[2]

julia> @finite_parameter(model, par2, 42)
par2
```
"""
macro finite_parameter(model, args...)
    # define error message function
    _error(str...) = JuMP._macro_error(:finite_parameter, (model, args...),
                                       str...)
    # parse the arguments
    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)

    # check number of arguments
    if length(extra) == 0 || length(extra) > 2
        _error("Incorrect number of arguments. Must be of form " *
               "@finite_parameter(model, name_expr, value).")
    end

    # ensure unsupported keywords are not given
    # if length(kw_args) != 0
    #     _error("Unrecognized keyword arguments given. Only the container " *
    #            "keyword is recognized.")
    # end

    # parse the input information from `extra`
    x = popfirst!(extra)
    # we have a single anonymous expression
    if length(extra) == 0
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an " *
                                               "`InfiniteModel`."
            value = $x
            @assert isa(value, Number) "Value must be a number."
            @infinite_parameter($model, lower_bound = value, upper_bound = value,
                                supports = value, error_args = $args,
                                container = $requestedcontainer, ($(kw_args...)))
        end
    # we have a some sort of name expression or vector expression
    elseif isa(x, Symbol) || isexpr(x, :vect) || isexpr(x, :vcat) || isexpr(x, :ref)
        name_expr = x
        value_expr = extra[1]
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an " *
                                               "`InfiniteModel`."
            @infinite_parameter($model, $name_expr, lower_bound = $value_expr,
                                upper_bound = $value_expr, supports = $value_expr,
                                container = $requestedcontainer, ($(kw_args...)),
                                error_args = $args)
        end
    # we have some other syntax
    else
        _error("Unrecognized input format. Must be of form " *
               "@finite_parameter(model, value) or " *
               "@finite_parameter(model, name[i =..., ...], value_expr).")
    end
    esc(code)
end

## Define helper functions needed to parse infinite variable expressions
# Check rhs to to ensure is not a variable as best we can
function _check_rhs(arg1, arg2)
    if isexpr(arg2, :ref)
        if isexpr(arg2.args[2], :kw)
            temp = arg2
            arg2 = arg1
            arg1 = temp
        end
    # elseif isexpr(arg2, :call) --> infinite loop if both sides are :calls
    #     temp = arg2
    #     arg2 = arg1
    #     arg1 = temp
    end
    return arg1, arg2
end

# Assume variable on lhs, but check if swapped
function _less_than_parse(arg1, arg2)
    # check if swapped
    arg1_orig = arg1
    arg1, arg2 = _check_rhs(arg1, arg2)
    if arg1 != arg1_orig
        return _greater_than_parse(arg1, arg2)
    end
    # check if has parameter tuple
    if isexpr(arg1, :call)
        return Expr(:call, :<=, arg1.args[1], arg2),
                    Expr(:tuple, arg1.args[2:end]...)
    else
        return Expr(:call, :<=, arg1, arg2), nothing
    end
end

# Assume variable on lhs, but check if swapped
function _greater_than_parse(arg1, arg2)
    # check if swapped
    arg1_orig = arg1
    arg1, arg2 = _check_rhs(arg1, arg2)
    if arg1 != arg1_orig
        return _less_than_parse(arg1, arg2)
    end
    # check if has parameter tuple
    if isexpr(arg1, :call)
        return Expr(:call, :>=, arg1.args[1], arg2),
                    Expr(:tuple, arg1.args[2:end]...)
    else
        return Expr(:call, :>=, arg1, arg2), nothing
    end
end

# The variable is on the rhs
function _less_than_parse(arg1::Number, arg2)
    return _greater_than_parse(arg2, arg1)
end

# The variable is on the rhs
function _greater_than_parse(arg1::Number, arg2)
    return _less_than_parse(arg2, arg1)
end

# Assume variable on lhs
function _equal_to_parse(arg1, arg2)
    arg1, arg2 = _check_rhs(arg1, arg2)
    if isexpr(arg1, :call)
        return Expr(:call, :(==), arg1.args[1], arg2),
                    Expr(:tuple, arg1.args[2:end]...)
    else
        return Expr(:call, :(==), arg1, arg2), nothing
    end
end

# The variable is on the rhs
function _equal_to_parse(arg1::Number, arg2)
    return _equal_to_parse(arg2, arg1)
end

# Return the tuple of form (var_expr, param_tuple_expr)
function _parse_parameters(_error::Function, head::Val{:call}, args)
    if args[1] in [:<=, :≤]
        return _less_than_parse(args[2], args[3])
    elseif args[1] in [:>=, :≥]
        return _greater_than_parse(args[2], args[3])
    elseif args[1] == :(==)
        return _equal_to_parse(args[2], args[3])
    elseif !(args[1] in [:in, :∈])
        return args[1], Expr(:tuple, args[2:end]...)
    else
        first = args[1]
        _error("Invalid operator $first.")
    end
end

# Return the tuple of form (var_expr, param_tuple_expr)
function _parse_parameters(_error::Function, head::Val{:comparison}, args)
    if isexpr(args[3], :call)
        return Expr(:comparison, args[1:2]..., args[3].args[1], args[4:5]...),
                    Expr(:tuple, args[3].args[2:end]...)
    else
        return Expr(:comparison, args...), nothing
    end
end

"""
    @infinte_variable(model, kw_args...)

Add an *anonymous* infinite variable to the model `model` described by the
keyword arguments `kw_args` and returns the variable reference. Note that the
`parameter_refs` keyword is required in this case.

```julia
@infinite_variable(model, varexpr, args..., kw_args...)
```

Add an infinite variable to `model` described by the expression `var_expr`, the
positional arguments `args` and the keyword arguments `kw_args`. The expression
`varexpr` can either be (note that in the following the symbol `<=` can be used
instead of `≤` and the symbol `>=`can be used instead of `≥`) of the form:

- `varexpr` creating variables described by `varexpr`
- `varexpr ≤ ub` (resp. `varexpr ≥ lb`) creating variables described by
  `varexpr` with upper bounds given by `ub` (resp. lower bounds given by `lb`)
- `varexpr == value` creating variables described by `varexpr` with fixed values
   given by `value`
- `lb ≤ varexpr ≤ ub` or `ub ≥ varexpr ≥ lb` creating variables described by
  `varexpr` with lower bounds given by `lb` and upper bounds given by `ub`

The expression `varexpr` can be of the form:

- `varname` creating a scalar real variable of name `varname`
- `varname(params...)` creating a scalar real variable of name `varname` with
  infinite parameters `params...` see `parameter_refs` for format.
- `varname[...]` or `[...]` creating a container of variables.
- `varname[...](params...)` or `[...]` creating a container of variables with
  infinite parameters `params...` in the first case.

The recognized positional arguments in `args` are the following:

- `Bin`: Sets the variable to be binary, i.e. either 0 or 1.
- `Int`: Sets the variable to be integer, i.e. one of ..., -2, -1, 0, 1, 2, ...

The recognized keyword arguments in `kw_args` are the following:

- `parameter_refs`: This is mandatory if not specified in `varexpr`. Can be a
  single parameter reference, a single parameter array with parameters defined
  in the same call of [`@infinite_parameter`](@ref) (i.e., have same group ID),
  or a tuple where each element is either of the first two options listed.
- `base_name`: Sets the name prefix used to generate variable names. It
  corresponds to the variable name for scalar variable, otherwise, the
  variable names are set to `base_name[...]` for each index `...` of the axes
  `axes`. Furthermore, the parameter reference tuple is appended on the end of
  the name i.e., `base_name(params...)` or `base_name[...](params...)`.
- `lower_bound`: Sets the value of the variable lower bound.
- `upper_bound`: Sets the value of the variable upper bound.
- `start`: Sets the variable starting value used as initial guess in optimization.
- `binary`: Sets whether the variable is binary or not.
- `integer`: Sets whether the variable is integer or not.
- `container`: Specify the container type.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP, Distributions; model = InfiniteModel())
julia> @infinite_parameter(model, 0 <= t <= 1)
t

julia> @infinite_parameter(model, w[1:2] in Normal())
2-element Array{ParameterRef,1}:
 w[1]
 w[2]

julia> @infinite_variable(model, x(t, w) >= 0)
x(t, w)

julia> x = @infinite_variable(model, parameter_refs = (t, w), base_name = "x",
                              lower_bound = 0)
x(t, w)

julia> lb = [0, 1]; ub = [10, 12];

julia> @infinite_variable(model, lb[i] <= y[i = 1:2](t) <= ub[i], Int)
2-element Array{InfiniteVariableRef,1}:
 y[1](t)
 y[2](t)
```
"""
macro infinite_variable(model, args...)
    _error(str...) = JuMP._macro_error(:infinite_variable, (model, args...),
                                       str...)

    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)
    param_kw_args = filter(kw -> kw.args[1] == :parameter_refs, kw_args)

    # Check for easy case if it is anonymous single variable
    if length(extra) == 0
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an `InfiniteModel`."
            JuMP.@variable($model, ($(args...)), variable_type = Infinite,
                           macro_error = $_error)
        end
    else
        x = popfirst!(extra)

        # There are several cases to consider:
        # x                                                     | type of x | x.head
        # ------------------------------------------------------+-----------+------------
        # var                                                   | Symbol    | NA
        # var[1:2]                                              | Expr      | :ref
        # var(x, y)                                             | Expr      | :call
        # var[1:2](x, y)                                        | Expr      | :call
        # var <= ub or var[1:2] <= ub                           | Expr      | :call
        # var(x, y) <= ub or var[1:2](x, y) <= ub               | Expr      | :call
        # lb <= var <= ub or lb <= var[1:2] <= ub               | Expr      | :comparison
        # lb <= var(x, y) <= ub or lb <= var[1:2](x, y) <= ub   | Expr      | :comparison
        if isexpr(x, :comparison) || isexpr(x, :call)
            inf_expr, params = InfiniteOpt._parse_parameters(_error, Val(x.head),
                                                             x.args)
        else
            inf_expr = x
            params = nothing
        end

        # check for double specification of parameters
        if length(param_kw_args) != 0 && params != nothing
            _error("Cannot specify double specify the infinite parameter " *
                   "references.")
        end

        # easy case where we can don't need to parse extra args
        if length(args) == 1
            code = quote
                @assert isa($model, InfiniteModel) "Model must be an " *
                                                   "`InfiniteModel`."
                JuMP.@variable($model, ($(inf_expr)),
                               variable_type = Infinite,
                               parameter_refs = $params,
                               macro_error = $_error)
            end
        # here we need to parse the extra args and include them in the call
        else
            rest_args = [args[i] for i = 2:length(args)]
            if isa(params, Nothing)
                code = quote
                    @assert isa($model, InfiniteModel) "Model must be an " *
                                                       "`InfiniteModel`."
                    JuMP.@variable($model, ($(inf_expr)), ($(rest_args...)),
                                   variable_type = Infinite,
                                   macro_error = $_error)
                end
            else
                code = quote
                    @assert isa($model, InfiniteModel) "Model must be an " *
                                                       "`InfiniteModel`."
                    JuMP.@variable($model, ($(inf_expr)), ($(rest_args...)),
                                   variable_type = Infinite,
                                   parameter_refs = $params,
                                   macro_error = $_error)
                end
            end
        end
    end
    return esc(code)
end

"""
    @point_variable(model, kw_args...)

Add an *anonymous* point variable to the model `model` described by the
keyword arguments `kw_args` and returns the variable reference. Note that the
`infinite_variable_ref` and `parameter_values` keywords are required in this
case.

```julia
@point_variable(model, infvarexpr, varexpr, args..., kw_args...)
```

Add a point variable to `model` described by the expression `varexpr`, the
positional arguments `args`, and the keyword arguments `kw_args` and the
infinite variable expr `infvarexpr`. The expression `infvarexpr` specifies the
infinite variable this point variable corresponds to and the values at which the
parameters are evaluated and must be of the form: `infvar(param_values...)`
where the parameter values `param_values...` are listed in the same format as
they are in teh definition of `infvar`. The expression `varexpr` is used to
define variable specific bounds and whose name is used as an alias for the point
variable which is simply the infinite variable evaluated at the values indicated.
The expression `varexpr` can either be (note that in the following the symbol
`<=` can be used instead of `≤` and the symbol `>=`can be used instead of `≥`)
of the form:

- `varexpr` creating variables described by `varexpr`
- `varexpr ≤ ub` (resp. `varexpr ≥ lb`) creating variables described by
  `varexpr` with upper bounds given by `ub` (resp. lower bounds given by `lb`)
- `varexpr == value` creating variables described by `varexpr` with fixed values
   given by `value`
- `lb ≤ varexpr ≤ ub` or `ub ≥ varexpr ≥ lb` creating variables described by
  `varexpr` with lower bounds given by `lb` and upper bounds given by `ub`

 Note that be default a point variable inherits all of the same properties as
 the infinite variable it corresponds to, but that these can be overwritten
 by specifying properties such as lower bounds, fix values, etc.

The expression `varexpr` can be of the form:

- `varname` creating a scalar real variable of alias name `varname`
- `varname[...]` or `[...]` creating a container of variables.

The recognized positional arguments in `args` are the following:

- `Bin`: Sets the variable to be binary, i.e. either 0 or 1.
- `Int`: Sets the variable to be integer, i.e. one of ..., -2, -1, 0, 1, 2, ...

The recognized keyword arguments in `kw_args` are the following:

- `infinite_variable_ref`: Sets the infinite variable reference that the point
  variable is associated with.
- `parameter_refs`: Sets the values of the infinite parameters of the infinite
  variable at which this poitn variable is evaluated at. Must be of the same
  format of that specified for the parameters in the definition of the infinite
  variable.
- `base_name`: Sets the name prefix used to generate variable names. It
  corresponds to the variable name for scalar variable, otherwise, the
  variable names are set to `base_name[...]` for each index `...` of the axes
  `axes`. This serves as the alias for `infvarexpr` (the infinite variable
  evaluated at particular parameter values).
- `lower_bound`: Sets the value of the variable lower bound.
- `upper_bound`: Sets the value of the variable upper bound.
- `start`: Sets the variable starting value used as initial guess in optimization.
- `binary`: Sets whether the variable is binary or not.
- `integer`: Sets whether the variable is integer or not.
- `container`: Specify the container type.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP, Distributions; model = InfiniteModel())
julia> @infinite_parameter(model, 0 <= t <= 1)
t

julia> @infinite_parameter(model, w[1:2] in Normal())
2-element Array{ParameterRef,1}:
 w[1]
 w[2]

julia> @infinite_variable(model, x(t, w) >= 0)
x(t, w)

julia> @point_variable(model, x(0, [0, 0]), x0 <= 1)
x0

julia> x0 = @point_variable(model, x(0, [0, 0]), upper_bound = 1, base_name = "x0")
x0

julia> x0 = @point_variable(model, upper_bound = 1, base_name = "x0",
                            infinite_variable_ref = x, parameter_values = (0, [0, 0]))
x0

julia> @point_variable(model, x([0, 1][i], [0, 0]), xf[i = 1:2])
2-element Array{PointVariableRef,1}:
 xf[1]
 xf[2]

julia> lb = [0, 1]; ub = [10, 12];

julia> @infinite_variable(model, lb[i] <= y[i = 1:2](t) <= ub[i], Int)
2-element Array{InfiniteVariableRef,1}:
 y[1](t)
 y[2](t)

julia> @point_variable(model, y[i](0), y0[i = 1:2], Bin)
2-element Array{PointVariableRef,1}:
 y0[1]
 y0[2]

julia> y0 = @point_variable(model, [i = 1:2], binary = true, base_name = "y0",
                             infinite_variable_ref = y[i], parameter_values = 0)
2-element Array{PointVariableRef,1}:
 y0[1]
 y0[2]
```
"""
macro point_variable(model, args...)
    _error(str...) = JuMP._macro_error(:point_variable,
                                       (model, args...), str...)

    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)
    param_kw_args = filter(kw -> kw.args[1] == :infinite_variable_ref || kw.args[1] == :parameter_values, kw_args)

    # Check for easy case if it is anonymous single variable
    if length(extra) == 0
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an " *
                                               "`InfiniteModel`."
            JuMP.@variable($model, ($(args...)), variable_type = Point,
                           macro_error = $_error)
        end
    else
        x = popfirst!(extra)

        if isexpr(x, :call)
            if length(param_kw_args) != 0
                _error("Cannot double specify the infinite variable reference " *
                       "and/or its parameter values.")
            elseif x.args[1] in [:in, :<=, :>=, :(==), :≥, :≤, ∈]
                _error("Invalid input syntax.")
            end
            rest_args = [args[i] for i = 2:length(args)]
            inf_var = x.args[1]
            param_vals = Expr(:tuple, x.args[2:end]...)
            code = quote
                @assert isa($model, InfiniteModel) "Model must be an " *
                                                   "`InfiniteModel`."
                JuMP.@variable($model, ($(rest_args...)), variable_type = Point,
                               infinite_variable_ref = $inf_var,
                               parameter_values = $param_vals,
                               macro_error = $_error)
            end
        elseif isexpr(x, :vect) && length(extra) == 0
            code = quote
                @assert isa($model, InfiniteModel) "Model must be an " *
                                                   "`InfiniteModel`."
                JuMP.@variable($model, ($(args...)), variable_type = Point,
                               macro_error = $_error)
            end
        else
            _error("Invalid input syntax.")
        end
    end
    return esc(code)
end

"""
    @hold_variable(model, kw_args...)

Add an *anonymous* hold variable to the model `model` described by the
keyword arguments `kw_args` and returns the variable reference.

```julia
@hold_variable(model, varexpr, args..., kw_args...)
```

Add a hold variable to `model` described by the expression `varexpr`, the
positional arguments `args` and the keyword arguments `kw_args`. The expression
`varexpr` can either be (note that in the following the symbol `<=` can be used
instead of `≤` and the symbol `>=`can be used instead of `≥`) of the form:

- `varexpr` creating variables described by `varexpr`
- `varexpr ≤ ub` (resp. `varexpr ≥ lb`) creating variables described by
  `varexpr` with upper bounds given by `ub` (resp. lower bounds given by `lb`)
- `varexpr == value` creating variables described by `varexpr` with fixed values
   given by `value`
- `lb ≤ varexpr ≤ ub` or `ub ≥ varexpr ≥ lb` creating variables described by
  `varexpr` with lower bounds given by `lb` and upper bounds given by `ub`

The expression `varexpr` can be of the form:

- `varname` creating a scalar real variable of name `varname`
- `varname[...]` or `[...]` creating a container of variables.

The recognized positional arguments in `args` are the following:

- `Bin`: Sets the variable to be binary, i.e. either 0 or 1.
- `Int`: Sets the variable to be integer, i.e. one of ..., -2, -1, 0, 1, 2, ...

Specifiying a hold variable which applies only to sub-domain of the model's
infinite parameter(s) domain can be done via the `parameter_bounds` keyword
argument. It is specified as a tuple of parameter bound expressions which can
be of the form:

- `(param in [lb, ub], ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub` (note `∈` can be used in place of `in`)
- `(params in [lb, ub], ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(lb <= param <= ub, ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub`
- `(lb <= params <= ub, ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(param == value, ...)` enforcing `param` to be equal to `value`
- `(params == value, ...)` enforcing that all parameter references in `params`
                            each to be equal to `value`
- Any combination of the above forms. Must be inside parentheses and comma
  separated.

The other recognized keyword arguments in `kw_args` are the following:

- `base_name`: Sets the name prefix used to generate variable names. It
  corresponds to the variable name for scalar variable, otherwise, the
  variable names are set to `base_name[...]` for each index `...` of the axes
  `axes`.
- `lower_bound`: Sets the value of the variable lower bound.
- `upper_bound`: Sets the value of the variable upper bound.
- `start`: Sets the variable starting value used as initial guess in optimization.
- `binary`: Sets whether the variable is binary or not.
- `integer`: Sets whether the variable is integer or not.
- `container`: Specify the container type.

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

julia> @hold_variable(model, 0 <= y <= 4, Bin)
y

julia> y = @hold_variable(model, lower_bound = 0, upper_bound = 4,
                            binary = true, base_name = "y")
y

julia> @hold_variable(model, z[2:3] == 0)
1-dimensional DenseAxisArray{HoldVariableRef,1,...} with index sets:
    Dimension 1, 2:3
And data, a 2-element Array{HoldVariableRef,1}:
 z[2]
 z[3]

julia> @hold_variable(model, d, parameter_bounds = (t in [0, 5]))
d
```
"""
macro hold_variable(model, args...)
    _error(str...) = JuMP._macro_error(:hold_variable, (model, args...),
                                       str...)
    # parse the arguments
    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)
    bound_kw_args = filter(kw -> kw.args[1] == :parameter_bounds, kw_args)

    # no bounds given so we don't need to do anything
    if length(bound_kw_args) == 0
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an " *
                                               "`InfiniteModel`."
            JuMP.@variable($model, ($(args...)), variable_type = Hold,
                           macro_error = $_error)
        end
    else
        x = bound_kw_args[1].args[2]
        if isexpr(x, :call) || isexpr(x, :tuple) || isexpr(x, :comparison)
           name_expr, bounds = InfiniteOpt._extract_bounds(_error, x.args,
                                                           Val(x.head))
        else
           _error("Unrecognized input format for parameter bounds. Must be of " *
                  "tuple with elements of form: par(s) in [lb, ub] or " *
                  "par(s) = value.")
        end
        extra_kw_args = filter(kw -> kw.args[1] != :parameter_bounds, kw_args)
        code = quote
            @assert isa($model, InfiniteModel) "Model must be an " *
                                               "`InfiniteModel`."
            JuMP.@variable($model, ($(extra...)), variable_type = Hold,
                           parameter_bounds = ($(bounds)), macro_error = $_error,
                           container = ($(requestedcontainer)),
                           ($(extra_kw_args...)))
        end
    end
    return esc(code)
end

## Helper function for parsing the parameter bounds and the name expression
# Check that name expression is in acceptable form or throw error
function _parse_name_expression(_error::Function, expr)
    if isexpr(expr, :ref) || isexpr(expr, :vect)|| isa(expr, Symbol) || isexpr(expr, :vcat)
        return expr
    else
        _error("Unrecognized input format for name expression. Must be of " *
               "form name[i = ..., ...](param_bounds).")
    end
end

# Return IntervalSet call given a :vect expression or error
function _make_interval_set(_error::Function, expr::Expr)
    if isexpr(expr, :vect) && length(expr.args) == 2
        return Expr(:call, :IntervalSet, expr.args...)
    else
        _error("Unrecognized input format for parameter bounds. Must be of form " *
        "par in [lb, ub], lb <= par <= ub,  or par = value.")
    end
end

# Return ParameterRef(s) as a vector if necessary (to handle JuMP arrays)
function _process_parameter(arg)
    return Expr(:call, :(InfiniteOpt._make_vector), arg)
end

## Return a bound pair of form param => IntervalSet(lb, ub) or error
# Call expressions using :in
function _make_bound_pair(_error::Function, expr::Expr,
                          ::Union{Val{:in}, Val{:∈}})
    set = _make_interval_set(_error, expr.args[3])
    param = _process_parameter(expr.args[2])
    return Expr(:call, :(=>), param, set)
end

# Call expressions using :(==)
function _make_bound_pair(_error::Function, expr::Expr, ::Val{:(==)})
    set = _make_interval_set(_error, Expr(:vect, expr.args[3], expr.args[3]))
    param = _process_parameter(expr.args[2])
    return Expr(:call, :(=>), param, set)
end

# Comparison expressions using something else
function _make_bound_pair(_error::Function, expr::Expr, first)
    _error("Invalid set format operators. Must be of form par in [lb, ub], " *
           "lb <= par <= ub, or par = value.")
end

# Comparison expressions using :(<=)
function _make_bound_pair(_error::Function, expr::Expr, ::Val{:(<=)},
                          ::Val{:(<=)})
    set = _make_interval_set(_error, Expr(:vect, expr.args[1], expr.args[5]))
    param = _process_parameter(expr.args[3])
    return Expr(:call, :(=>), param, set)
end

# Comparison expressions using :(>=)
function _make_bound_pair(_error::Function, expr::Expr, ::Val{:(>=)},
                          ::Val{:(>=)})
    set = _make_interval_set(_error, Expr(:vect, expr.args[5], expr.args[1]))
    param = _process_parameter(expr.args[3])
    return Expr(:call, :(=>), param, set)
end

# Comparison expressions using something else
function _make_bound_pair(_error::Function, expr::Expr, first, second)
    _error("Invalid set format operators. Must be of form par in [lb, ub], " *
           "lb <= par <= ub, or par = value.")
end

# General wrapper
function _make_bound_pair(_error::Function, expr::Expr)
    if isexpr(expr, :call)
        return _make_bound_pair(_error, expr, Val(expr.args[1]))
    elseif isexpr(expr, :comparison)
        return _make_bound_pair(_error, expr, Val(expr.args[2]), Val(expr.args[4]))
    else
        _error("Unrecognized input format for parameter bounds. Must be of form " *
               "par in [lb, ub], lb <= par <= ub, or par = value.")
    end
end

# Return a dictionary of parameter bounds given raw vector of call expressions
function _parse_parameter_bounds(_error::Function, args::Vector)
    dict_args = [_make_bound_pair(_error, arg) for arg in args]
    return Expr(:call, :ParameterBounds, Expr(:call, :Dict, dict_args...))
end

# Only 1 parameter bound is given thus dispatch as a vector to make dictionary
function _parse_parameter_bounds(_error::Function, expr::Expr)
    return _parse_parameter_bounds(_error, [expr])
end

## Return name expression and parameter dictionary expression separatedly
# (:call) In this case we are not sure if a name expression is given
function _extract_bounds(_error::Function, args::Vector, ::Val{:call})
    # check if only one parameter bound was given and nothing else
    if args[1] in [:in, :∈, :(==)]
        bounds = _parse_parameter_bounds(_error, Expr(:call, args...))
        return nothing, bounds
    # otherwise parameter bounds were given with a name expression
    else
        name_expr = _parse_name_expression(_error, args[1])
        bounds = _parse_parameter_bounds(_error, args[2:end])
        return name_expr, bounds
    end
end

# (:tuple) In this case we know mulitple bounds are given with no name expression
function _extract_bounds(_error::Function, args::Vector, ::Val{:tuple})
    bounds = _parse_parameter_bounds(_error, args)
    return nothing, bounds
end

# (:comparison) In this case we have a single interval set and nothing else
function _extract_bounds(_error::Function, args::Vector, ::Val{:comparison})
    bounds = _parse_parameter_bounds(_error, Expr(:comparison, args...))
    return nothing, bounds
end

"""
    @BDconstraint(model::InfiniteModel, [i = ..., ...](bound_expr), constr_expr;
                  [kw_args...])

Add an anonymous bounded constraint to `model` and return an appropriate
container of constraint reference(s).

```julia
    @BDconstraint(model::InfiniteModel, name[i = ..., ...](bound_expr),
                  constr_expr; [kw_args...])
```

Add a named bounded constraint to `model` and return an appropriate
container of constraint reference(s). This defines the constraint as expressed
in `constr_expr` over some sub-domain of parameters as indicated by `bound_expr`.
The format of `constr_expr` must follow the same syntax as that specified in
[`JuMP.@constraint`](@ref). For example, if we want to express the constraint
`2T(t, x) + 3y = 42`, the constraint expression `constr_expr` would be
`2 * T + 3 * y == 42`.

By default, `JuMP.@constraint` would express the above example constraint over
the whole domain of `t` and `x`. However, we can use `@BDconstraint` to express
this constraint over some sub-domain(s) via `bound_expr`. Here `(bound_expr)`
can be of the form:

- `(param in [lb, ub], ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub` (note `∈` can be used in place of `in`)
- `(params in [lb, ub], ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(lb <= param <= ub, ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub`
- `(lb <= params <= ub, ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(param == value, ...)` enforcing `param` to be equal to `value`
- `(params == value, ...)` enforcing that all parameter references in `params`
                            each to be equal to `value`
- Any combination of the above forms. Must be inside parentheses and comma
  separated.

Like typical constraints, the `container` keyword argument can be used to specify
the `JuMP` container type used to store the constraint references. Note this
macro errors if `bound_expr` is ommited or if some unrecognized syntax is used.

**Examples**
```julia-repl
julia> @BDconstraint(model, c1(t in [0, 1]), T^2 + z <= 1)
c1 : T(x, t)² + z ≤ 1.0, ∀ t ∈ [0, 1]

julia> @BDconstraint(model, c2[i = 1:3](x[i] in [0, 1]), T^2 + z + x[i] <= 1)
3-element Array{GeneralConstraintRef,1}:
 c2[1] : T(x, t)² + z + x[1] ≤ 1.0, ∀ x[1] ∈ [0, 1]
 c2[2] : T(x, t)² + z + x[2] ≤ 1.0, ∀ x[2] ∈ [0, 1]
 c2[3] : T(x, t)² + z + x[3] ≤ 1.0, ∀ x[3] ∈ [0, 1]

julia> @BDconstraint(model, (x == 0, t == 0), T^2 + z <= 1)
T(x, t)² + z ≤ 1.0, ∀ x[2] = 0, x[3] = 0, t = 0, x[1] = 0

julia> @BDconstraint(model, [i = 1:3](x[i] == 0), T^2 + z <= 1,
                     container = SparseAxisArray)
JuMP.Containers.SparseAxisArray{GeneralConstraintRef,1,Tuple{Any}} with 3 entries:
  [3]  =  T(x, t)² + z ≤ 1.0, ∀ x[3] = 0
  [2]  =  T(x, t)² + z ≤ 1.0, ∀ x[2] = 0
  [1]  =  T(x, t)² + z ≤ 1.0, ∀ x[1] = 0
```
"""
macro BDconstraint(model, args...)
    # define appropriate error message function
    _error(str...) = JuMP._macro_error(:BDconstraint, (model, args...), str...)

    # parse the arguments
    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)
    bound_kw_args = filter(kw -> kw.args[1] == :parameter_bounds, kw_args)

    # check for double specification of parameter bounds
    if length(bound_kw_args) != 0
        _error("Cannot double specify the parameter bounds. Use form " *
               "@BDconstraint(model, name[i=..., ...](par in [lb, ub], " *
               "par2 = value, ...), expr).")
    end

    # check that enough arguments are given
    if length(extra) != 2
        _error("Incorrect amount of arguments. Must be of form " *
               "@BDconstraint(model, name[i=..., ...](par in [lb, ub], " *
               "par2 = value, ...), expr).")
    # otherwise we have a 3 argument form
    else
        # process the 2nd argument for the parameter bounds if provided
        x = popfirst!(extra)
        if isexpr(x, :ref) || isexpr(x, :vect) || isa(x, Symbol) || isexpr(x, :vcat)
            _error("Must specify at least one parameter bound of the form " *
                   "@BDconstraint(model, name[i=..., ...](par in [lb, ub], " *
                   "par2 = value, ...), expr).")
        elseif isexpr(x, :call) || isexpr(x, :tuple) || isexpr(x, :comparison)
            name_expr, bounds = InfiniteOpt._extract_bounds(_error, x.args,
                                                            Val(x.head))
        else
            _error("Unrecognized input format for name expression. Must be of " *
                   "form @BDconstraint(model, name[i=..., ...](par in [lb, ub], " *
                   "par2 = value, ...), expr).")
        end
        # TODO check the expression type somewhere...
        # we have a single anonymous constraint
        if isa(name_expr, Nothing)
            code = quote
                @assert isa($model, InfiniteModel) "Model must be an " *
                                                   "`InfiniteModel`."
                JuMP.@constraint($model, ($(extra[1])),
                                 parameter_bounds = ($(bounds)), ($(kw_args...)),
                                 container = ($(requestedcontainer)))
            end
        # we have some sort of name expression in our constraint
        else
            code = quote
                @assert isa($model, InfiniteModel) "Model must be an " *
                                                   "`InfiniteModel`."
                JuMP.@constraint($model, ($(name_expr)), ($(extra[1])),
                                 parameter_bounds = ($(bounds)), ($(kw_args...)),
                                 container = ($(requestedcontainer)))
            end
        end
    end
    return esc(code)
end

# TODO add looping capability for multi-dimensional refs
"""
    @set_parameter_bounds(ref, bound_expr; [force = false])

Specify new parameter bounds for a constraint reference or hold variable
reference `ref`. These bounds correspond to bounding a constraint in an
equivalent way to using [`@BDconstraint`](@ref) or to limiting the scope of a
hold variable in an equivalent way to using the `parameter_bounds` keyword
argument in [`@hold_variable`](@ref). Here `(bound_expr)` can be of the form:

- `(param in [lb, ub], ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub` (note `∈` can be used in place of `in`)
- `(params in [lb, ub], ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(lb <= param <= ub, ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub`
- `(lb <= params <= ub, ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(param == value, ...)` enforcing `param` to be equal to `value`
- `(params == value, ...)` enforcing that all parameter references in `params`
                            each to be equal to `value`
- Any combination of the above forms. Must be inside parentheses and comma
  separated.

Errors if the constraint or variable corresponding to `ref` already has bounds.
However, using `force = true` can be used ignore the current bounds and
overwrite them with new ones. Also, note that bounds on dependent constraints
of hold variables will be updated to account for changes in hold variable bounds.

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

julia> @infinite_variable(model, x(t))
x(t)

julia> @hold_variable(model, y)
y

julia> @constraint(model, con, x + y == 0)
con : x(t) + y = 0.0

julia> @set_parameter_bounds(y, t in [0, 5])

julia> con
con : x(t) + y = 0.0, ∀ t ∈ [0, 5]

julia> @set_parameter_bounds(con, t == 0, force = true)

julia> con
con : x(t) + y = 0.0, ∀ t = 0
```
"""
macro set_parameter_bounds(ref, bound_expr, args...)
    # define appropriate error message function
    _error(str...) = JuMP._macro_error(:set_parameter_bounds,
                                       (ref, bound_expr, args...), str...)

    # parse the arguments
    extra, kw_args, requestedcontainer = JuMPC._extract_kw_args(args)
    extra_kw_args = filter(kw -> kw.args[1] != :force, kw_args)
    if length(extra) != 0
        _error("Too many positional arguments.")
    elseif length(extra_kw_args) != 0
        _error("Unrecognized keyword arguments.")
    end

    # parse the bounds
    x = bound_expr
    if isexpr(x, :call) || isexpr(x, :tuple) || isexpr(x, :comparison)
       name_expr, bounds = InfiniteOpt._extract_bounds(_error, x.args,
                                                       Val(x.head))
    else
       _error("Unrecognized input format for parameter bounds. Must be of " *
              "tuple with elements of form: par(s) in [lb, ub] or " *
              "par(s) = value.")
    end
    # prepare the code
    code = quote
        @assert(isa($ref, Union{GeneralConstraintRef, HoldVariableRef}),
                "Reference must correspond to a constraint or hold variable.")
        set_parameter_bounds($ref, $bounds; ($(args...)), _error = $_error)
    end
    return esc(code)
end

"""
    @add_parameter_bounds(ref, bound_expr)

Add new parameter bounds for a constraint reference or hold variable
reference `ref`. These bounds correspond to bounding a constraint in an
equivalent way to using [`@BDconstraint`](@ref) or to limiting the scope of a
hold variable in an equivalent way to using the `parameter_bounds` keyword
argument in [`@hold_variable`](@ref). Here `(bound_expr)` can be of the form:

- `(param in [lb, ub], ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub` (note `∈` can be used in place of `in`)
- `(params in [lb, ub], ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(lb <= param <= ub, ...)` enforcing `param` to be in a sub-domain from `lb`
                             to `ub`
- `(lb <= params <= ub, ...)` enforcing that all parameter references in `params`
                              each be a in sub-domain from `lb` to `ub`
- `(param == value, ...)` enforcing `param` to be equal to `value`
- `(params == value, ...)` enforcing that all parameter references in `params`
                            each to be equal to `value`
- Any combination of the above forms. Must be inside parentheses and comma
  separated.

Errors if the new bounds cause irreconcilable differences with existing measures
and constraints. For example, this occurs when adding hold variable bounds
that are outside the domain of a bounded constraint that uses that hold variable.
Also, note that bounds on dependent constraints of hold variables will be
updated to account for changes in hold variable bounds.

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

julia> @infinite_parameter(model, q in [-2, 2])
q

julia> @infinite_variable(model, x(t, q))
x(t, q)

julia> @hold_variable(model, y)
y

julia> @constraint(model, con, x + y == 0)
con : x(t, q) + y = 0.0

julia> @add_parameter_bounds(y, t in [0, 5])

julia> con
con : x(t, q) + y = 0.0, ∀ t ∈ [0, 5]

julia> @add_parameter_bounds(con, q == 0)

julia> con
con : x(t, q) + y = 0.0, ∀ t ∈ [0, 5], q = 0
```
"""
macro add_parameter_bounds(ref, bound_expr)
    # define appropriate error message function
    _error(str...) = JuMP._macro_error(:add_parameter_bounds,
                                       (ref, bound_expr), str...)

    # parse the bounds
    x = bound_expr
    if isexpr(x, :call) || isexpr(x, :tuple) || isexpr(x, :comparison)
       name_expr, bounds = InfiniteOpt._extract_bounds(_error, x.args,
                                                       Val(x.head))
    else
       _error("Unrecognized input format for parameter bounds. Must be of " *
              "tuple with elements of form: par(s) in [lb, ub] or " *
              "par(s) = value.")
    end
    # prepare the code
    code = quote
        @assert(isa($ref, Union{GeneralConstraintRef, HoldVariableRef}),
                "Reference must correspond to a constraint or hold variable.")
        for (pref, set) in ($(bounds)).intervals
            add_parameter_bound($ref, pref, set.lower_bound, set.upper_bound,
                                 _error = $_error)
        end
    end
    return esc(code)
end


text:
function calculate_residuals(ws::COSMO.Workspace,  IGNORESCALING_FLAG::Bool = false)

	if (ws.settings.scaling != 0 && !IGNORESCALING_FLAG)
		r_prim = norm(ws.sm.Einv * (ws.p.A * ws.vars.x + ws.vars.s - ws.p.b), Inf)
		# ∇f0 + ∑ νi ∇hi(x*) == 0 condition
		r_dual = norm(ws.sm.cinv[] * ws.sm.Dinv * (ws.p.P * ws.vars.x + ws.p.q - ws.p.A' * ws.vars.μ), Inf)
	end
	if (ws.settings.scaling == 0 || IGNORESCALING_FLAG )
		r_prim = norm(ws.p.A * ws.vars.x + ws.vars.s - ws.p.b, Inf)
		r_dual = norm(ws.p.P * ws.vars.x + ws.p.q - ws.p.A' * ws.vars.μ, Inf)
	end
	# FIXME: Why is it -A'μ ?
	return r_prim,r_dual
end

function max_res_component_norm(ws::COSMO.Workspace, IGNORESCALING_FLAG::Bool = false)
	if (ws.settings.scaling != 0 && !IGNORESCALING_FLAG)
		max_norm_prim = max.(norm(ws.sm.Einv * (ws.p.A * ws.vars.x), Inf), norm(ws.sm.Einv * ws.vars.s, Inf), norm(ws.sm.Einv * ws.p.b, Inf))
		max_norm_dual = max.(norm(ws.sm.cinv .* (ws.sm.Dinv * (ws.p.P * ws.vars.x)), Inf), norm(ws.sm.cinv .* (ws.sm.Dinv * ws.p.q), Inf), norm(ws.sm.cinv .* (ws.sm.Dinv * (ws.p.A' * ws.vars.μ)), Inf))
	end
	if (ws.settings.scaling == 0 || IGNORESCALING_FLAG)
		max_norm_prim = max.(norm(ws.p.A * ws.vars.x, Inf), norm(ws.vars.s, Inf), norm(ws.p.b, Inf))
		max_norm_dual = max.(norm(ws.p.P * ws.vars.x, Inf), norm(ws.p.q, Inf), norm(ws.p.A' * ws.vars.μ, Inf))
	end
	return max_norm_prim, max_norm_dual
end

function has_converged(ws::COSMO.Workspace, r_prim::Float64, r_dual::Float64)
	max_norm_prim, max_norm_dual = max_res_component_norm(ws)
	settings = ws.settings
	ϵ_prim = settings.eps_abs + settings.eps_rel * max_norm_prim
	ϵ_dual = settings.eps_abs + settings.eps_rel * max_norm_dual

	# if an optimal objective value was specified for the problem check if current solution is within specified accuracy
	obj_true_FLAG = true
	if !isnan(settings.obj_true)
		current_cost = ws.sm.cinv[] * (1/2 * ws.vars.x' * ws.p.P * ws.vars.x + ws.p.q' * ws.vars.x)[1]
		obj_true_FLAG = abs(settings.obj_true - current_cost) <= settings.obj_true_tol
	end

	return (r_prim < ϵ_prim  && r_dual < ϵ_dual && obj_true_FLAG)
end

ext: jl


text:
## In this notebook, we'll look at least squares on biometric data
A = randn(10,5)
b = randn(10)
x = A\b # built in least-squares
## Solve via Normal equations
y = (A'*A) \ (A'*b)
@show norm(x-y)
## Solve via QR
Q,R = qr(A)
c = Q'*b
c1 = c[1:size(A,2)]
z = R\c1
@show norm(x-z)


text:
using Distributed
@time include("required.jl")

## Conditional arguments
# outcome = "raw" # desired outcome (required)
# outcome = "sdm" # desired outcome (required)
# save_figures = true # should figures be overwritten (optional)

# Make sure "outcome" is defined
if !(@isdefined outcome)
    @warn "'outcome' not defined, must be either 'raw' or 'sdm'"
elseif (outcome != "raw" && outcome != "sdm")
    @warn "'outcome' invalid, must be either 'raw' or 'sdm'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Load distributions for all species
@load joinpath("data", "jld2", "$(outcome)-distributions.jld2") distributions

## Create custom functions
# Function for Pielou's evenness index
function pielou(a::Vector)
    A = filter(!isnothing, a)
    length(A) == 0 && return nothing
    iszero(sum(A)) && return nothing
    p = A ./ sum(A)
    return abs(sum(filter(!isnothing, p .* log.(p))) / length(p))
end
function pielou_allspecies(a::Vector)
    A = filter(!isnothing, a)
    length(A) == 0 && return nothing
    iszero(sum(A)) && return nothing
    p = A ./ length(A)
    return abs(sum(filter(!isnothing, p .* log.(p))) / length(p))
end

# Function for Shannon's diversity index
function shannon(a::Vector)
    A = filter(!isnothing, a)
    length(A) == 0 && return nothing
    iszero(sum(A)) && return nothing
    p = A ./ sum(A)
    return abs(sum(p .* log.(p)))
end

## Calculate diversity/evenness scores
# Empty array for diversity scores
output = zeros(Float32, size(distributions[1])) |> Array{Union{Nothing,Float32}}
output2 = zeros(Float32, size(distributions[1])) |> Array{Union{Nothing,Float32}}
# Loop for each pixel/grid element
@time for i in 1:size(output, 1), j in 1:size(output, 2)
    # Group distributions for all species in pixel [i,j]
    x = map(x -> x.grid[i, j], distributions)
    # Calculate Pielou's evenness index for pixel [i,j]
    output[i, j] = pielou(x)
    output2[i, j] = pielou_allspecies(x)
end
# Create SimpleSDMLayer with diversity/evenness scores
diversity = SimpleSDMResponse(
    output,
    distributions[1].left,
    distributions[1].right,
    distributions[1].bottom,
    distributions[1].top,
)
diversity2 = SimpleSDMResponse(
    output2,
    distributions[1].left,
    distributions[1].right,
    distributions[1].bottom,
    distributions[1].top,
)

## Plot result
diversity_plot = plot_layer(
    diversity;
    c=:BuPu,
    title="Community evenness ($(outcome) distributions)",
    colorbar_title="Pielou's evenness index (site richness)",
)
diversity_plot2 = plot_layer(
    diversity2;
    c=:BuPu,
    title="Community evenness ($(outcome) distributions)",
    colorbar_title="Pielou's evenness index (total richness)",
)

## Save result
# save_figures = true # should figures be overwritten (optional)
if (@isdefined save_figures) && save_figures == true
    @info "Figures saved ($(outcome) evenness)"
    savefig(diversity_plot, joinpath("fig", outcome, "04_$(outcome)_diversity-pielou.png"))
    savefig(
        diversity_plot2, joinpath("fig", outcome, "04_$(outcome)_diversity-pielou2.png")
    )
else
    @info "Figures not saved ($(outcome) evenness)"
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"""
    AffineFEOperator

Reprepresent a fully assembled affine (linear) finite element problem.
See also [FEOperator](@ref)
"""
struct AffineFEOperator <: FEOperator
  trial::FESpace
  test::FESpace
  op::AffineOperator
end

"""
"""
function AffineFEOperator(trial::FESpace,test::FESpace,matrix::AbstractMatrix,vector::AbstractVector)
  @assert num_free_dofs(trial) == size(matrix,2) "Incompatible trial space and matrix"
  @assert num_free_dofs(test) == size(matrix,1) "Incompatible test space and matrix"
  op = AffineOperator(matrix,vector)
  AffineFEOperator(trial,test,op)
end

"""
   AffineFEOperator(test::FESpace,trial::FESpace,assem::Assembler,terms::AffineFETerm...)
   AffineFEOperator(test::FESpace,trial::FESpace,terms::AffineFETerm...)
"""
function AffineFEOperator(trial::FESpace,test::FESpace,assem::Assembler,terms::AffineFETerm...)

  u = get_cell_basis(trial)
  v = get_cell_basis(test)
  @assert is_trial(u)

  uhd = zero(trial)

  data = collect_cell_matrix_and_vector(uhd,u,v,terms)
  A,b = assemble_matrix_and_vector(assem,data)

  #matdata = collect_cell_matrix(u,v,terms)
  #vecdata = collect_cell_vector(uhd,v,terms)
  #A = assemble_matrix(assem,matdata)
  #b = assemble_vector(assem,vecdata)

  AffineFEOperator(trial,test,A,b)
end

function AffineFEOperator(trial::FESpace,test::FESpace,terms::AffineFETerm...)
  assem = SparseMatrixAssembler(trial,test)
  AffineFEOperator(trial,test,assem,terms...)
end

function AffineFEOperator(
  mat::Type{<:AbstractSparseMatrix},trial::FESpace,test::FESpace,terms::AffineFETerm...)
  assem = SparseMatrixAssembler(mat,trial,test)
  AffineFEOperator(trial,test,assem,terms...)
end

function AffineFEOperator(
  mat::Type{<:AbstractSparseMatrix},vec::Type{<:AbstractVector},trial::FESpace,test::FESpace,terms::AffineFETerm...)
  assem = SparseMatrixAssembler(mat,vec,trial,test)
  AffineFEOperator(trial,test,assem,terms...)
end

get_test(feop::AffineFEOperator) = feop.test

get_trial(feop::AffineFEOperator) = feop.trial

"""
"""
get_matrix(feop::AffineFEOperator) = get_matrix(feop.op)

"""
"""
get_vector(feop::AffineFEOperator) = get_vector(feop.op)

get_algebraic_operator(feop::AffineFEOperator) = feop.op

function allocate_residual(feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  allocate_residual(feop.op,x)
end

function residual!(b::AbstractVector,feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  residual!(b,feop.op,x)
end

function residual(feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  residual(feop.op,x)
end

function allocate_jacobian(feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  allocate_jacobian(feop.op,x)
end

function jacobian!(A::AbstractMatrix,feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  jacobian!(A,feop.op,x)
end

function jacobian(feop::AffineFEOperator,u)
  @assert is_a_fe_function(u)
  x = get_free_values(u)
  jacobian(feop.op,x)
end


text:
module OneDimensionalCosineAdvectionDiffusion

using Printf, Statistics

using Oceananigans, Oceananigans.OutputWriters, Oceananigans.Grids

include("analysis.jl")

# Advection and diffusion of a cosine.
c(x, y, z, t, U, κ) = exp(-κ * t) * cos(x - U * t)

function run_test(; Nx, Δt, stop_iteration, U = 1, κ = 1e-4,
                  architecture = CPU(), topo = (Periodic, Periodic, Bounded))

    #####
    ##### Test c and v-advection
    #####

    grid = RegularCartesianGrid(size=(Nx, 1, 1), x=(0, 2π), y=(0, 1), z=(0, 1), topology=topo)

    model = IncompressibleModel(architecture = architecture,
                                        grid = grid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     closure = ConstantIsotropicDiffusivity(ν=κ, κ=κ))

    set!(model, u = U,
                v = (x, y, z) -> c(x, y, z, 0, U, κ),
                c = (x, y, z) -> c(x, y, z, 0, U, κ))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress_frequency=stop_iteration)

    println("Running 1D in x cosine advection diffusion test for v and c with Nx = $Nx and Δt = $Δt...")
    run!(simulation)

    x = xnodes(model.tracers.c)
    c_analytical = c.(x, 0, 0, model.clock.time, U, κ)

    # Calculate errors
    cx_simulation = model.tracers.c
    cx_simulation = interior(cx_simulation)[:, 1, 1]
    cx_errors = compute_error(cx_simulation, c_analytical)

    v_simulation = model.velocities.v
    v_simulation = interior(v_simulation)[:, 1, 1]
    v_errors = compute_error(v_simulation, c_analytical)

    #####
    ##### Test u-advection
    #####

    ygrid = RegularCartesianGrid(size=(1, Nx, 1), x=(0, 1), y=(0, 2π), z=(0, 1), topology=topo)

    model = IncompressibleModel(architecture = architecture,
                                        grid = ygrid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     closure = ConstantIsotropicDiffusivity(ν=κ, κ=κ))

    set!(model, v = U,
                u = (x, y, z) -> c(y, x, z, 0, U, κ),
                c = (x, y, z) -> c(y, x, z, 0, U, κ))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress_frequency=stop_iteration)

    println("Running 1D in y cosine advection diffusion test for u and c with Ny = $Nx and Δt = $Δt...")
    run!(simulation)

    # Calculate errors
    u_simulation = model.velocities.u
    u_simulation = interior(u_simulation)[1, :, 1]
    u_errors = compute_error(u_simulation, c_analytical)

    # Calculate errors
    cy_simulation = model.tracers.c
    cy_simulation = interior(cy_simulation)[1, :, 1]
    cy_errors = compute_error(cy_simulation, c_analytical)

    return (

            cx = (simulation = cx_simulation,
                  analytical = c_analytical,
                          L₁ = cx_errors.L₁,
                          L∞ = cx_errors.L∞),

            cy = (simulation = cy_simulation,
                  analytical = c_analytical,
                          L₁ = cy_errors.L₁,
                          L∞ = cy_errors.L∞),

            v = (simulation = v_simulation,
                 analytical = c_analytical, # same solution as c.
                         L₁ = v_errors.L₁,
                         L∞ = v_errors.L∞),

            u = (simulation = u_simulation,
                 analytical = c_analytical, # same solution as c.
                         L₁ = u_errors.L₁,
                         L∞ = u_errors.L∞),

            grid = grid

            )
end

end # module


text:
import RestrictedBoltzmannMachines as RBMs
import ConvolutionalRBMs as ConvRBMs

using Test: @test, @testset, @inferred, @test_throws
using Random: bitrand, randn!
using LinearAlgebra: dot
using RestrictedBoltzmannMachines: visible, hidden, weights, Binary
using RestrictedBoltzmannMachines: energy, free_energy, interaction_energy
using RestrictedBoltzmannMachines: inputs_v_to_h, inputs_h_to_v
using RestrictedBoltzmannMachines: sample_v_from_v, sample_h_from_h, sample_v_from_h, sample_h_from_v
using ConvolutionalRBMs: vsizes, hsizes, BinaryConvRBM, output_size, out2in

@testset "pooling ConvRBM" begin
    rbm = BinaryConvRBM(3, 2, 3; pool=true)
    randn!(weights(rbm))
    randn!(visible(rbm).θ)
    randn!(hidden(rbm).θ)

    v = bitrand(3,5)
    @test @inferred(free_energy(rbm, v)) isa Number
    v = bitrand(3,5,2)
    @test size(@inferred free_energy(rbm, v)) == (2,)
    v = bitrand(3,5)
    @test @inferred(output_size(rbm, v)) == (3,)

    v = bitrand(3,5,2)
    @test size(@inferred sample_h_from_v(rbm, v)) == (2,3,2)
    @test @inferred(output_size(rbm, v)) == (3,)

    h1 = @inferred sample_h_from_v(rbm, v)
    h2 = @inferred sample_h_from_h(rbm, h1)
    @test size(h1) == size(h2) == (2,3,2)
    # pooling: at most one hidden unit active per group
    @test all(sum(h1; dims=2) .≤ 1)
    @test all(sum(h2; dims=2) .≤ 1)
end


text:
# Julia wrapper for header: /usr/include/scip/pub_disp.h
# Automatically generated using Clang.jl wrap_c

function SCIPdispGetData(disp)
    ccall((:SCIPdispGetData, libscip), Ptr{SCIP_DISPDATA}, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispSetData(disp, dispdata)
    ccall((:SCIPdispSetData, libscip), Cvoid, (Ptr{SCIP_DISP}, Ptr{SCIP_DISPDATA}), disp, dispdata)
end

function SCIPdispGetName(disp)
    ccall((:SCIPdispGetName, libscip), Cstring, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetDesc(disp)
    ccall((:SCIPdispGetDesc, libscip), Cstring, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetHeader(disp)
    ccall((:SCIPdispGetHeader, libscip), Cstring, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetWidth(disp)
    ccall((:SCIPdispGetWidth, libscip), Cint, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetPriority(disp)
    ccall((:SCIPdispGetPriority, libscip), Cint, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetPosition(disp)
    ccall((:SCIPdispGetPosition, libscip), Cint, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispGetStatus(disp)
    ccall((:SCIPdispGetStatus, libscip), SCIP_DISPSTATUS, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispIsInitialized(disp)
    ccall((:SCIPdispIsInitialized, libscip), UInt32, (Ptr{SCIP_DISP},), disp)
end

function SCIPdispLongint(messagehdlr, file, val, width)
    ccall((:SCIPdispLongint, libscip), Cvoid, (Ptr{SCIP_MESSAGEHDLR}, Ptr{FILE}, Clonglong, Cint), messagehdlr, file, val, width)
end

function SCIPdispInt(messagehdlr, file, val, width)
    ccall((:SCIPdispInt, libscip), Cvoid, (Ptr{SCIP_MESSAGEHDLR}, Ptr{FILE}, Cint, Cint), messagehdlr, file, val, width)
end

function SCIPdispTime(messagehdlr, file, val, width)
    ccall((:SCIPdispTime, libscip), Cvoid, (Ptr{SCIP_MESSAGEHDLR}, Ptr{FILE}, Cdouble, Cint), messagehdlr, file, val, width)
end


text:
module UniformMatroid

using Kombinator

using JuMP

import Base: copy
import Kombinator: solve, dimension, value, values, formulation, fastest_exact, make_solution

export UniformMatroidInstance,
    UniformMatroidSolution,
    MinBudgetedUniformMatroidSolution,
    items,
    items_all_budgets

include("um.jl")

include("um_dp.jl")
include("um_greedy.jl")
include("um_lp.jl")

include("um_budgeted_dp.jl")
include("um_budgeted_lp.jl")

end


text:
"""
    _cross(parent1::Vector{Vector{T}}, parent2::Vector{Vector{T}}) where T <: Signed
Simulate a mating cross between two parents, generating one offspring with the same
ploidy as `parent1`. This variant of `cross` is used internally for `simulatekin`.
"""
function _cross(parent1::Vector{Vector{T}}, parent2::Vector{Vector{T}}) where T <: Signed
    ploidy = length(first(parent1))
    ploidy == 1 && error("Haploid crosses are not yet supported. Please file and issue or pull request")
    if ploidy == 2
        p1_contrib = rand.(parent1)
        p2_contrib = rand.(parent2)
        geno_out = sort.(zip(p1_contrib, p2_contrib))
    elseif iseven(ploidy)
        n_allele = ploidy ÷ 2
        p1_contrib = sample.(parent1, n_allele, replace = false)
        p2_contrib = sample.(parent2, n_allele, replace = false)
        geno_out = Tuple.(sort!.(append!.(p1_contrib, p2_contrib)))
    else
        # special method to provide a 50% chance of one parent giving more alleles than the other
        rng = rand()
        contrib_1 = ploidy ÷ 2
        contrib_2 = ploidy - contrib_1
        p1_contrib = rng > 0.5 ? sample.(parent1, contrib_1, replace = false) : sample.(parent1, contrib_2, replace = false)
        p2_contrib = rng > 0.5 ? sample.(parent2, contrib_2, replace = false) : sample.(parent2, contrib_1, replace = false)
        geno_out = Tuple.(sort!.(append!.(p1_contrib, p2_contrib)))
    end
    return geno_out
end

function _parentoffspring(alleles::Dict, loc::Vector{String}, n::Int, ploidy::Signed, padding::Int)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim" * lpad(i, padding, '0')
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        insertcols!(out_df, Symbol(prefix * "_parent") =>  Tuple.(p1))
        insertcols!(out_df, Symbol(prefix * "_offspring") => _cross(p1, p2))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "parent_offspring")
    return out_df
end

function _fullsib(alleles::Dict, loc::Vector{String}, n::Int, ploidy::Signed, padding::Int)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim" * lpad(i, padding, '0')
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        [insertcols!(out_df, Symbol(prefix * "_fullsib_$j") => _cross(p1, p2)) for j in 1:2]
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "fullsib")
    return out_df
end

function _halfsib(alleles::Dict, loc::Vector{String}, n::Int, ploidy::Signed, padding::Int)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim" * lpad(i, padding, '0')
        p1,p2,p3 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:3]
        insertcols!(out_df, Symbol(prefix * "_halfsib_1") => _cross(p1, p2))
        insertcols!(out_df, Symbol(prefix * "_halfsib_2") => _cross(p1, p3))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "halfsib")
    return out_df
end

function _unrelated(alleles::Dict, loc::Vector{String}, n::Int, ploidy::Signed, padding::Int)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim" * lpad(i, padding, '0')
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        insertcols!(out_df, Symbol(prefix * "_unrelated_1") => Tuple.(p1))
        insertcols!(out_df, Symbol(prefix * "_unrelated_2") => Tuple.(p2))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "unrelated")
    return out_df
end

"""
    simulatekin(data::PopData; fullsib::Int, halfsib::Int, unrelated::Int, parentoffspring::Int, ploidy::Signed)
Simulate mating crosses to generate sample pairs with any combination of the specified relationships,
returning a `PopData` object. The simulations will first generate parents of a given
`ploidy` (inferred or specified) by drawing alleles from a global allele pool derived
from the given `data` (i.e. weighted by their frequencies).

#### Relationship
Simulated parents will be crossed to generate offspring depending on the relationship:
- `fullsib` : 2 parents generate 2 full-sibling offspring, return 2 offspring
- `halfsib` : 3 parents generate 2 half-sibling offspring, returns 2 offspring
- `unrelated` : returns 2 randomly generated individuals from the global allele pools
- `parentoffspring` : 2 parents generate 1 offspring, returns 1 offspring and 1 parent

#### Identifying pairs
The relationship between the newly generated samples can be identified by:
- Sample `name`s will specify their simulation number, relationship, and whether parent or offspring
    - Naming convention: [simulation #]_[relationship]_[offspring #]
    - example: sim005_fullsib_1 = [simulation 005]_[full sibling]_[offspring 1]
- Their `population` name will be that of their relationship (e.g. "fullsib")

#### Ploidy
If the samples in your `PopData` are of a single ploidy, then `ploidy = 0` (the default) will infer the ploidy
and generate parents and offspring according to the ploidy of your data. If you have mixed-ploidy data or wish
to generate parents and offspring of a ploidy different than the source `PopData` you can specify the ploidy
with which to simulate parents and offspring. For example, if your `PopData` is diploid, but you wish to generate
triploid or octoploid parents and offspring, you would specify `ploidy = 3` or `ploidy = 8` repectively.
#### Odd ploidy
If trying to create offspring with an odd ploidy (3,5, etc.), each parent has a 50% chance of
contributing (½ × ploidy) + 1 alleles for all loci to the offspring. In other words, if ploidy = 3,
there's a 50% chance parent_1 will give 2 alleles for every locus for that simulated cross.

**Example**
```
julia> cats = @nanycats ;

julia> cat_sims = simulatekin(cats, fullsib = 10, halfsib = 50)
PopData{Diploid, 9 Microsatellite loci}
  Samples: 120
  Populations: 2

julia> cat_sims.sampleinfo
120×3 DataFrame
 Row │ name             population  ploidy
     │ String           String      Int64
─────┼─────────────────────────────────────
   1 │ sim01_fullsib_1  fullsib          2
   2 │ sim01_fullsib_2  fullsib          2
   3 │ sim02_fullsib_1  fullsib          2
   4 │ sim02_fullsib_2  fullsib          2
   5 │ sim03_fullsib_1  fullsib          2
   6 │ sim03_fullsib_2  fullsib          2
  ⋮  │        ⋮             ⋮         ⋮
 115 │ sim48_halfsib_1  halfsib          2
 116 │ sim48_halfsib_2  halfsib          2
 117 │ sim49_halfsib_1  halfsib          2
 118 │ sim49_halfsib_2  halfsib          2
 119 │ sim50_halfsib_1  halfsib          2
 120 │ sim50_halfsib_2  halfsib          2
                           108 rows omitted
```
"""
function simulatekin(data::PopData; fullsib::Int = 0, halfsib::Int = 0, unrelated::Int = 0, parentoffspring::Int = 0, ploidy::Signed = 0)
    if iszero(sum([fullsib, halfsib, unrelated, parentoffspring]))
        throw(ArgumentError("Please specify at least one of: \n- \"fullsib\" \n- \"halfsib\" \n- \"unrelated\"\n- \"parentoffspring\""))
    end
    # automatic ploidy finding
    if ploidy == 0
        ploids = data.metadata.ploidy
        if ploids isa AbstractVector
            error("For mixed ploidy data, please specify a single ploidy with which to generate parents and offspring")
        else
            ploidy += ploids
        end
    end
    loc, alleles = allele_pool(data)
    # how many digits to pad the offspring names
    padding = length(string(maximum([fullsib, halfsib, unrelated, parentoffspring])))
    # perform the simulation if the integer > 0, otherwise return an empty boolean vector
    # the empty vector is just to skip over with Base.Iterators.filter
    fs = fullsib > 0 ? _fullsib(alleles, loc, fullsib, ploidy, padding) : Bool[]
    hs = halfsib > 0 ? _halfsib(alleles, loc, halfsib, ploidy, padding) : Bool[]
    unrl = unrelated > 0 ? _unrelated(alleles, loc, unrelated, ploidy, padding) : Bool[]
    poff = parentoffspring > 0 ? _parentoffspring(alleles, loc, parentoffspring, ploidy, padding) : Bool[]
    # combine the results together into a single df
    geno_df = reduce(vcat, Base.Iterators.filter(!isempty, (fs, hs, unrl, poff)))
    geno_df.name = PooledArray(geno_df.name, compress = true)
    geno_df.population = PooledArray(geno_df.population, compress = true)
    geno_df.locus = PooledArray(geno_df.locus, compress = true)
    #meta_df = select(unique(geno_df, :name), 1, 2)
    #insertcols!(meta_df, :ploidy => ploidy)
    return PopData(geno_df)
end


text:
# ------------ Union of constraints ------------ #

struct ConstrainedUnionSelection <: ConstrainedSelection
    query::Vector{ConstrainedSelection}
end

function has_top(cus::ConstrainedUnionSelection, addr)
    for q in cus.query
        has_top(q, addr) && return true
    end
    return false
end
function dump_queries(cus::ConstrainedUnionSelection)
    arr = Address[]
    for q in cus.query
        append!(arr, collect(dump_queries(q)))
    end
    return arr
end
function get_top(cus::ConstrainedUnionSelection, addr)
    for q in cus.query
        has_top(q, addr) && return get_top(q, addr)
    end
    error("ConstrainedUnionSelection (get_top): query not defined for $addr.")
end

function get_sub(cus::ConstrainedUnionSelection, addr)
    return ConstrainedUnionSelection(map(cus.query) do q
                                         get_sub(q, addr)
                                     end)
end

isempty(cus::ConstrainedUnionSelection) = foldr(x -> isempty(x), cus.query)

# ------------ Unconstrained union selection ------------ #

struct UnconstrainedUnionSelection <: UnconstrainedSelection
    query::Vector{UnconstrainedSelection}
end

function has_top(uus::UnconstrainedUnionSelection, addr)
    for q in uus.query
        has_top(q, addr) && return true
    end
    return false
end

function dump_queries(uus::UnconstrainedUnionSelection)
    arr = Address[]
    for q in uus.query
        append!(arr, collect(dump_queries(q)))
    end
    return arr
end

function get_sub(uus::UnconstrainedUnionSelection, addr)
    return UnconstrainedUnionSelection(map(uus.query) do q
                                           get_sub(q, addr)
                                       end)
end

isempty(uus::UnconstrainedUnionSelection) = isempty(uus.query)


text:
## To control debugging information, set these constants and recompile

# debug level, larger means more verbose. -1 is off.
const MXDEBUGLEVEL = -1

# if true, print filenames as they are included.
const VERBOSEINCLUDE = false

# Print symbol names ("functions") as the dispatch code is written via @doap
const VERBOSE_DOAP = false

const PYDEBUGLEVEL = -1

const SJDEBUGLEVEL = -1

###

macro mdebug(level, a...)
    if level <= MXDEBUGLEVEL
        esc(:((println($(a...)));println()))
    else
        nothing
    end
end

macro inc(filename)
    if VERBOSEINCLUDE
        :(println("--- loading ", $filename); include($filename))
    else
        :(include($filename))
    end
end


text:
#=
These tests require a serial device to echo bytes written from the host computer.

For example, these hardware configurations should work:
 - A standalone USB-to-UART adapter (e.g. FTDI FT232) in loopback mode (TX and RX pins jumped).
 - A USB-to-serial equipped, arduino-compatible microcontroller running examples/serial_example.ino.
 - The USB-to-serial chip on some microcontroller boards can be used directly, bypassing the micro.
   For example, an Arduino UNO R3 can be connected via USB with RESET grounded and TX wired to RX.

Run the tests locally using
$ julia test/runtests.jl <port address> <baudrate>

LibSerialPort.list_ports() may help with finding the correct port address.
=#

using LibSerialPort
using LibSerialPort.Lib
using Test

if haskey(ENV, "CI")
    @testset "LibSerialPort CI port listing" begin
        list_ports()
    end
else
    if length(ARGS) == 0
        port = "/dev/ttyS0"
    else
        port = ARGS[1]
    end

    baudrate = length(ARGS) >= 2 ? ARGS[2] : "115200"

    @testset "LibSerialPort" begin
        @testset "Low level API" begin
            include("test-low-level-api.jl")
            @test test_low_level_api() == nothing
            @test test_low_level_api(port, baudrate) == nothing
        end

        @testset "High level API" begin
            include("test-high-level-api.jl")
            @test test_high_level_api() == nothing
            @test test_high_level_api(port, baudrate) == nothing
        end

    end
end


text:
using DataStructures, Test, Suppressor

graph = Dict()
graph["you"] = ["alice","bob","claire"]
graph["bob"] = ["anuj","peggy"]
graph["claire"] = ["thom","jonny"]
graph["alice"] = ["peggy"]
graph["anuj"] = []
graph["peggy"] = []
graph["thom"] = []
graph["jonny"] = []

function person_is_seller(name)
	return name[end] == 'm'
end

function search(name)
	search_queue = Deque{String}()
	map(n -> push!(search_queue,n),graph[name])
	searched = []
	while isempty(search_queue) == false
		person = popfirst!(search_queue)
		if (person in searched) == false
			if person_is_seller(person)
				println(person , " is a mango seller")
				return true
			else
				map(n -> push!(search_queue,n),graph[person])
				push!(searched,person)
			end
		end
	end
	return false
end

result = @capture_out search("you")
@test result == "thom is a mango seller
"


text:
mutable struct Frame{E}
    entities::Vector{E} # NOTE: I tried StaticArrays; was not faster
    n::Int
end
function Frame{E}(arr::AbstractVector{E}, N::Int=length(arr))
    N ≥ length(arr) || error("capacity cannot be less than entitiy count! (N ≥ length(arr))")
    entities = convert(Vector{E}, arr)
    return Frame{E}(entities, N)
end
function Frame{E}(::Type{E}, N::Int=100)
    entities = Array{E}(N)
    return Frame{E}(entities, 0)
end

Base.show{E}(io::IO, frame::Frame{E}) = @printf(io, "Frame{%s}(%d entities)", string(E), length(frame))

capacity(frame::Frame) = length(frame.entities)
Base.length(frame::Frame) = frame.n
Base.getindex(frame::Frame, i::Int) = frame.entities[i]

Base.endof(frame::Frame) = frame.n
function Base.setindex!{E}(frame::Frame{E}, entity::E, i::Int)
    frame.entities[i] = entity
    return frame
end
function Base.empty!(frame::Frame)
    frame.n = 0
    return frame
end
function Base.deleteat!(frame::Frame, entity_index::Int)
    entity_index > 0 || throw(DomainError())

    for i in entity_index : frame.n - 1
        frame.entities[i] = frame.entities[i+1]
    end
    frame.n -= 1
    frame
end

Base.start(frame::Frame) = 1
Base.done(frame::Frame, i::Int) = i > length(frame)
Base.next(frame::Frame, i::Int) = (frame.entities[i], i+1)

function Base.copy!{E}(dest::Frame{E}, src::Frame{E})
    for i in 1 : src.n
        dest.entities[i] = src.entities[i]
    end
    dest.n = src.n
    return dest
end
Base.copy{E}(frame::Frame{E}) = copy!(Frame(E, capacity(frame)), frame)

function Base.push!{E}(frame::Frame{E}, entity::E)
    frame.n += 1
    frame.entities[frame.n] = entity
    return frame
end

####

const EntityFrame{S,D,I} = Frame{Entity{S,D,I}}
EntityFrame{S,D,I}(::Type{S},::Type{D},::Type{I}) = Frame(Entity{S,D,I})
EntityFrame{S,D,I}(::Type{S},::Type{D},::Type{I},N::Int) = Frame(Entity{S,D,I}, N)

Base.in{S,D,I}(frame::EntityFrame{S,D,I}, id::I) = findfirst(frame, id) != 0
function Base.findfirst{S,D,I}(frame::EntityFrame{S,D,I}, id::I)
    for entity_index in 1 : frame.n
        entity = frame.entities[entity_index]
        if entity.id == id
            return entity_index
        end
    end
    return 0
end
function id2index{S,D,I}(frame::EntityFrame{S,D,I}, id::I)
    entity_index = findfirst(frame, id)
    if entity_index == 0
        throw(BoundsError(frame, id))
    end
    return entity_index
end
get_by_id{S,D,I}(frame::EntityFrame{S,D,I}, id::I) = frame[id2index(frame, id)]
function get_first_available_id{S,D,I}(frame::EntityFrame{S,D,I})
    ids = Set{I}(entity.id for entity in frame)
    id_one = one(I)
    id = id_one
    while id ∈ ids
        id += id_one
    end
    return id
end
function Base.push!{S,D,I}(frame::EntityFrame{S,D,I}, s::S)
    id = get_first_available_id(frame)
    entity = Entity{S,D,I}(s, D(), id)
    push!(frame, entity)
end

Base.delete!{S,D,I}(frame::EntityFrame{S,D,I}, entity::Entity{S,D,I}) = deleteat!(frame, findfirst(frame, entity.id))
function Base.delete!{S,D,I}(frame::EntityFrame{S,D,I}, id::I)
    entity_index = findfirst(frame, id)
    if entity_index != 0
        deleteat!(frame, entity_index)
    end
    return frame
end

###

function Base.write{S,D,I}(io::IO, mime::MIME"text/plain", frames::Vector{EntityFrame{S,D,I}})
    println(io, length(frames))
    for frame in frames
        println(io, length(frame))
        for entity in frame
            write(io, mime, entity.state)
            print(io, "\n")
            write(io, mime, entity.def)
            print(io, "\n")
            write(io, mime, entity.id)
            print(io, "\n")
        end
    end
end
function Base.read{S,D,I}(io::IO, mime::MIME"text/plain", ::Type{Vector{EntityFrame{S,D,I}}})

    n = parse(Int, readline(io))
    frames = Array{EntityFrame{S,D,I}}(n)

    for i in 1 : n
        m = parse(Int, readline(io))
        frame = Frame(Entity{S,D,I}, m)
        for j in 1 : m
            state = read(io, mime, S)
            def = read(io, mime, D)
            id = read(io, mime, I)
            push!(frame, Entity(state,def,id))
        end
        frames[i] = frame
    end

    return frames
end


text:
# the logic is that we generate script from the example so we can utilize #src and #!src things
using Literate

const is_ci = haskey(ENV, "GITHUB_ACTIONS")

# generate files using literate.jl
examples_dir = @__DIR__
examples_generated_dir = joinpath(@__DIR__, "generated")
!ispath(examples_generated_dir) && mkpath(examples_generated_dir)
example_files = [joinpath(examples_dir, f) for f in [
    "recipes.jl",
    "schema_examination.jl",
    "schema_visualization.jl",
]]

for f in example_files
    Literate.script(f, examples_generated_dir)
end


text:
include("inputs_struct.jl")
include("constraints/constraints.jl")
include("neighbourhood_struct.jl")
include("output_struct.jl")
include("models/models.jl")


text:
# Creates sparse spd band matrix w/ only upper triangle stored
function spgenband( m::Int64, b::Int64 )
  nnonzeros = div((2m-b)*(b+1),2)
  nzval = randn(nnonzeros)
  rowval = Array(Int64,nnonzeros)
  colptr = Array(Int64,m+1)
  k = 1
  for i = 1:m
    colptr[i] = k
    for j = max(i-b,1):i
      rowval[k] = j
      k = k+1
    end
  end
  colptr[m+1] = k
  return SparseMatrixCSC(m,m,colptr,rowval,nzval) + m*speye(m,m)
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    return vec(readdlm(fname, ',', Int))
end

function Median(data::Array)
    if isodd(length(data))
        x0 = sort(data)[div(length(data)-1,2)]
        x1 = sort(data)[div(length(data)+1,2)]

        return x0, x1
    else
        x0 = sort(data)[div(length(data),2)]
        return x0
    end
end

# the optimal position will be the median
function Part1(data::Array)
    if isodd(length(data))
        x0, x1 = Median(data)

        c0 = sum(abs(x0-x) for x in data)
        c1 = sum(abs(x1-x) for x in data)

        return min(c0, c1)
    else
        x0 = Median(data)
        return sum(abs(x0-x) for x in data)
    end
end

# the optimal position will be somewhere between the average and the median
function Part2(data::Array)
    xAvg = round(Int, sum(data)/length(data))
    xMed = isodd(length(data)) ? min(Median(data)) : Median(data)

    # search between xAvg and xMed for the minimum
    x0, x1 = min(xAvg, xMed), max(xAvg, xMed)
    cMin = sum(data)*(sum(data)+1)
    cost = x -> div(sum(abs(y-x)*(1+abs(y-x)) for y in data),2)
    for x in x0:x1
        c = cost(x)
        cMin = c < cMin ? c : cMin
    end

    return cMin
end

#data = [16,1,2,0,4,2,7,1,2,14]
data = ReadData("../data/day7.txt")

println(Part1(data))
println(Part2(data))


text:
mutable struct DenseLinearSolver{T<:AbstractFloat} <: AbstractLinearSolver{T}
    Abuf::Vector{Matrix{T}}
    Fbuf::Vector{CholeskyPivoted{T, Matrix{T}}}
end

struct RecursiveLinearSolver{T<:AbstractFloat} <: AbstractLinearSolver{T}
    Abuf::Vector{Matrix{T}}
    Fbuf::Vector{LU{T, Matrix{T}}}
    threshold::Int
end

struct SparseLinearSolver{T<:AbstractFloat} <: AbstractLinearSolver{T}
    A::SparseMatrixCSC{T, Int}
    Fsym::SuiteSparse.CHOLMOD.Factor{T}
    Fbuf::Vector{SuiteSparse.CHOLMOD.Factor{T}}
    perm::Union{Nothing, Vector{Int}}
end

getjac(solver::AbstractLinearSolver, n) = nothing
getjac(solver::DenseLinearSolver, n) = @inbounds solver.Abuf[n]
getjac(solver::RecursiveLinearSolver, n) = @inbounds solver.Abuf[n]
getjac(solver::SparseLinearSolver, n) = solver.A

setjac!(solver::AbstractLinearSolver, n) = nothing
setjac!(solver::SparseLinearSolver, n) = !isnothing(solver.perm) && @inbounds solver.Fbuf[n] = copy(solver.Fsym)

getperm(prob, ::Val{S}) where {S} = nothing
getperm(prob, ::Val{:none}) = 1:size(prob.C, 2)
getperm(prob, ::Val{:amd}) = AMD.symamd(prob.C' * prob.C)
getperm(prob, ::Val{:metis}) = Int.(Metis.permutation(prob.C' * prob.C)[1])

function DenseLinearSolver{T}(prob::NonlinearProblem{T}) where {T}
    Nt = size(prob.C, 2)
    Nd = length(prob.pviews)
    Abuf = map(_ -> diagm(ones(T, Nt)), 1:Nd)
    Fbuf = map(A -> cholesky!(A, Val(true)), Abuf)
    DenseLinearSolver{T}(Abuf, Fbuf)
end

function RecursiveLinearSolver{T}(prob::NonlinearProblem{T}; threshold=100) where {T}
    Nt = size(prob.C, 2)
    Nd = length(prob.pviews)
    Abuf = map(_ -> diagm(ones(T, Nt)), 1:Nd)
    Fbuf = map(A -> RecursiveFactorization.lu!(A, Val(true); threshold), Abuf)
    RecursiveLinearSolver{T}(Abuf, Fbuf, threshold)
end

function SparseLinearSolver{T}(prob::NonlinearProblem{T}; reorder::Symbol=:none) where {T}
    Nt = size(prob.C, 2)
    Nd = length(prob.pviews)
    A = spdiagm(ones(T, Nt))
    perm = getperm(prob, Val(reorder)) # matrix reordering scheme
    Fsym = cholesky(A; perm) # symbolic factorization
    Fbuf = fill(Fsym, Nd)
    SparseLinearSolver{T}(A, Fsym, Fbuf, perm)
end

function solve!(x, b, n, solver::DenseLinearSolver)
    @unpack Abuf, Fbuf = solver
    @inbounds Fbuf[n] = cholesky!(Abuf[n], Val(true))
    # TODO: Somehow faster than 'ldiv!'
    x .= @inbounds Fbuf[n] \ b
    return x
end

function solve!(x, b, n, solver::RecursiveLinearSolver)
    @unpack Abuf, Fbuf, threshold = solver
    @inbounds Fbuf[n] = RecursiveFactorization.lu!(Abuf[n], Val(true); threshold)
    # TODO: Somehow faster than 'ldiv!'
    x .= @inbounds Fbuf[n] \ b
    return x
end

function solve!(x, b, n, solver::SparseLinearSolver)
    @unpack A, Fsym, Fbuf, perm = solver
    if isnothing(perm)
        @inbounds Fbuf[n] = cholesky(A)
        x .= @inbounds Fbuf[n] \ b
    else
        # FIXED: If a non-empty 'perm' was given, then the reordering method is already contained in 'Fsym'
        cholesky!(Fsym, A)
        x .= Fsym \ b
    end
    return x
end


text:
# This file was generated, do not modify it. # hide
rpt2.best_report.controls |> collect

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials